

In [4]:

    
from scipy.integrate import ode
from models import nmda, LR, LR_ICa, HH, na_current, cable
#from sympy import plot_implicit, symbols, Eq
#import sympy
%pylab inline
%load_ext autoreload
%autoreload 2









    



Populating the interactive namespace from numpy and matplotlib



In [26]:

    
def nmda_impulse(t, y, T=10, t_in=500, t_end=1500):
    rb = 0.005
    dy = nmda(y)
    if t_in <= t < t_end:
        dy[0] += -rb * T * y[0]
        dy[1] += rb * T * y[0] - rb * y[1] * T
        dy[2] += rb * y[1]
    return dy

def ica(O, V):
    Mg = 2.0
    gNMDA = 1.0
    VCa = 100.0
    B = 1.0 / (1 + exp(-0.062 * V)*Mg/3.57)
    ICa = gNMDA * B * O * (V - VCa)
    return ICa

def v_dendrite(t, y):
    dy = zeros_like(y)
    dy[:4] = nmda_impulse(t, y[:4])
    ICa = ica(y[3], y[4])
    dy[4] = -ICa - 0.3*y[4]
    return dy


def nmda_LR(t, y):
    dy = zeros_like(y)
    dy[:4] = nmda_impulse(t, y[:4])
    ICa = ica(y[3], y[4])
    dy[4:] = LR(y[4:], ICa)
    #dy[6] = 0.007*(-0.0275* ICa - 0.02*y[6]) 
    return dy
    
def nmda_HH(t, y):
    dy = zeros_like(y)
    dy[:4] = HH(y[:4])
    dy[4:] = nmda_impulse(t, y[4:])
    Ica = ca_current(y[7], y[0])
    dy[0] -= Ica
    return dy

solver = ['dopri5', 'dop853', 'vode']

Hodgkin and Huxley



In [3]:

    
x, y = symbols('x y')
gna = 120.0
gk = 36.0
gl = 0.3
Vna = 115.0
Vk = -12.0
Vl = 10.6
alpha_m = 0.1 * (25.0 - x) / (sympy.exp((25.0 - x) / 10.0) - 1.0)
beta_m = 4.0 * sympy.exp(-x / 18.0)
m_inf = alpha_m / (alpha_m + beta_m)
INa = gna * m_inf**3 * (0.8 - y) * (x - Vna)
IK = gk * y**4 * (x - Vk)
IL = gl * (x - Vl)


plot_implicit(Eq(INa + IK + IL, 0), (x, -20, 120), (y, 0, 1), adaptive=False);



In [44]:

    
#r = ode(  lambda t, y: HH(y) + array([5, 0, 0, 0]) if 5<= t < 7 else HH(y)).set_integrator(solver[0])
HH_input = lambda t, y: HH(y) + array([140, 0, 0, 0]) if 100<= t < 1100 else HH(y)
#print(HH_input(array(1, 0, 0, 0),0))
y0 = array([0.0, 0.3176, 0.053, 0.596])
r = ode(HH_input).set_integrator(solver[0])
r.set_initial_value(y0, 0)
t1 = 1200.0
dt = .01
v = []
time = []
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
    v.append(r.y[0])
    time.append(r.t)
figure(1)
plot(time, v)
xlabel('timepo (ms)')
ylabel('V (mV)')
grid(True)

Rinzel and Lee

$C_m \frac{dv}{dt}=-I_{Ca} −\left(g_K n^4 +\frac{g_{K, Ca} c}{K_d + c}\right) (V − V K ) − g_L (V − V_L )$

$\tau_n \frac{dn}{dt}=n_{\inf}(V)-n$

$\frac{dc}{dt} = f(-k_1 I_{Ca} -k_c c)$



In [13]:

    
#y0 = array([-50.0, 0.530, 0.7])
#r = ode(lambda t, y: LR(y)).set_integrator(solver[0])
y0 = array([-50.0, 0.530, 0.1])
r = ode(lambda t, y: array([LR(y)[0], LR(y)[1], 0])).set_integrator(solver[0])
r.set_initial_value(y0, 0)
t1 = 1000.0
dt = 1.0
v = []
n = []
ca = []
time = []
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
    v.append(r.y[0])
    n.append(r.y[1])
    ca.append(r.y[2])
    time.append(r.t)
figure(1)
plot(time, v)
xlabel('timepo (ms)')
ylabel('V (mV)')
grid(True)

figure(2)
plot(time, n)
xlabel('timepo (ms)')
ylabel('n')
grid(True)

figure(3)
plot(time, ca)
xlabel('timepo (ms)')
ylabel('[Ca] (M)')
grid(True)

show()

NMDA Channel



In [82]:

    
y0 = array([1, 0, 0, 0, 0])
#r = ode(lambda t, y: nmda_impulse(t, y)).set_integrator(solver[0])
r = ode(lambda t, y: v_dendrite(t, y)).set_integrator(solver[0])
r.set_initial_value(y0, 0)
t1 = 4000.0
dt = 1.0
o = []
ICa = []
v = []
time = []
while r.successful() and r.t < t1:
    r.integrate(r.t + dt)
    o.append(r.y[3])
    v.append(r.y[4])
    time.append(r.t)
figure(1)
plot(time, o)
grid(True)
xlabel('timepo (ms)')
ylabel('Fraccion de canales abiertos')

figure(2)
plot(time, v)
grid(True)
xlabel('timepo (ms)')
ylabel('V')


show()



In [87]:

    
y0 = array([1, 0, 0, 0, -50.0, 0.530, 0.4])
r = ode(nmda_LR).set_integrator(solver[0])
r.set_initial_value(y0, 0)
t1 = 4000.0
dt = 1.0
v = []
ca = []
time = []
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
    v.append(r.y[4])
    ca.append(r.y[6])
    time.append(r.t)
figure(1)
plot(time, v)
xlabel('timepo (s)')
ylabel('V (mV)')
grid(True)
figure(2)
plot(time, ca)
grid(True)

show()



In [28]:

    
y0 = array([0.0, 0.3176, 0.053, 0.596])
r = ode(lambda t, y: HH(y) + array([5, 0, 0, 0]) if 5<= t < 7 else HH(y)).set_integrator(solver[0])
r.set_initial_value(y0, 0)
t1 = 20
dt = 0.01
v = []
time = []
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
    v.append(r.y[0])
    time.append(r.t)
figure(1)
plot(time, v)
xlabel('timepo (ms)')
ylabel('V (mV)')
grid(True)



In [27]:

    
y0 = array([-40.0, 0.3176, 0.053, 0.596, 1, 0, 0, 0])
r = ode(nmda_HH).set_integrator(solver[0])
r.set_initial_value(y0, 0)
t1 = 2000.0
dt = 1.0
v = []
o = []
ICa = []
time = []
while r.successful() and r.t < t1: #
    r.integrate(r.t+dt)
    v.append(r.y[0])
    o.append(r.y[7])
    ICa.append(ica(r.y[7], r.y[0]))
    time.append(r.t)
figure(1)
plot(time, v)
xlabel('timepo (ms)')
ylabel('V (mV)')
grid(True)

figure(2)
plot(time, o)
xlabel('timepo (ms)')
ylabel('fraction open')
grid(True)

figure(3)
plot(time, ICa)
xlabel('timepo (ms)')
ylabel('Calcium Current')
grid(True)

Cable



In [16]:

    
y0 = zeros((100,))
y0[40:50] = 10 * ones(10,)
r = ode(lambda t, y: cable(y)).set_integrator(solver[0])
r.set_initial_value(y0, 0)
t1 = 10.0
dt = 1.0
v = [y0]
time = []
while r.successful() and r.t < t1: #
    r.integrate(r.t+dt)
    v.append(r.y)
    time.append(r.t)
figure(1)
for x in v:
    plot(x)
xlim([30, 50])
ylim([0, 0.01])
xlabel('distancia')
ylabel('V (mV)')
grid(True)



In [18]:

    
a = linspace(0, 9, 10)
N = a.shape[0]
print(a[1:N-1].shape)
print(a[:N-2].shape)
print(a[2:].shape)









    



(8,)
(8,)
(8,)