Trabajo Práctico de Biofísica 2C 2014

El objetivo de este IPython notebook es reproducir la figura 1d de [1] con un modelo computacional.

1 Jia, H., Rochefort, N. L., Chen, X., & Konnerth, A. (2010). Dendritic organization of sensory input to cortical neurons in vivo. Nature, 464(7293), 1307-1312.


In [1]:
from scipy.integrate import ode
%pylab inline


Populating the interactive namespace from numpy and matplotlib

Canales NMDA

La dinámica de los canales NMDA se puede describir con el siguiente diagrama de reacciones $$ C + T \rightleftarrows C_1 + T \rightleftarrows C_2 \rightleftarrows O $$

$$C_2 \rightleftarrows I$$

In [2]:
def nmda(y):
    ru = .0129
    rd = 0.0084
    rr = 0.0068
    r0 = 0.0465
    rc = 0.0738
    C1 = y[1]
    C2 = y[2]
    O = y[3]
    D = 1 - sum(y)
    dy = np.zeros_like(y)
    dy[0] = ru * C1
    dy[1] = -ru * C1 + ru * C2
    dy[2] = -ru * C2 - rd * C2 + rr * D - r0 * C2 + rc * O
    dy[3] = r0 * C2 - rc * O
    return dy

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

y0 = array([1, 0, 0, 0])
r = ode(lambda t, y: nmda_impulse(t, y)).set_integrator('dopri5')
r.set_initial_value(y0, 0)
t1 = 4000.0
dt = 1.0
o = []
ICa = []
time = []
while r.successful() and r.t < t1:
    r.integrate(r.t + dt)
    o.append(r.y[3])
    time.append(r.t)
figure(1, figsize=[10, 6])
plot(time, o)
fill_between([500, 1500],[0.14, 0.14], alpha=0.3)
ylim([0, 0.14])
grid(True)
xlabel('Tiempo (ms)', fontsize=14)
ylabel('Fracción de canales abiertos', fontsize=14)
show()


Corriente de calcio como corriente aplicada al modelo HH


In [4]:
def ca_current(O, V):
    Mg = 2.0
    gNMDA = 1.0
    VCa = 100.0
    B = 1.0 / (1 + np.exp(-0.062 * V)*Mg/3.57)
    ICa = gNMDA * B * O * (V - VCa)
    return ICa

def HH(y, Vna=115, Vk=-12.0, Vl=10.6):
    gna = 120.0
    gk = 36.0
    gl = 0.3
    Cm = 1.0
    dy = np.zeros_like(y)
    V = y[0]
    n = y[1]
    m = y[2]
    h = y[3]
    alpha_m = 0.1 * (25.0 - V) / (np.exp((25.0 - V) / 10.0) - 1.0)
    beta_m = 4.0 * np.exp(-V / 18.0)
    alpha_h = 0.07 * np.exp(-V / 20.0)
    beta_h = 1.0 / (1.0 + np.exp((30.0 - V) / 10.0))
    alpha_n = 0.01 * (10.0 - V) / (np.exp((10.0 - V) / 10.0) - 1.0)
    beta_n = 0.125 * np.exp(-V / 80.0)
    INa = gna * m**3 * h * (V - Vna)
    IK = gk * n**4 * (V - Vk)
    IL = gl * (V - Vl)
    dy[0] = -(INa + IK + IL) / Cm
    dy[1] = alpha_n * (1 - n) - beta_n * n
    dy[2] = alpha_m * (1 - m) - beta_m * m
    dy[3] = alpha_h * (1 - h) - beta_h * h
    return dy


def nmda_HH(t, y):
    dy = zeros_like(y)
    dy[:4] = HH(y[:4], Vk=-2.0, Vl=-40.0)
    dy[4:] = nmda_impulse(t, y[4:])
    Ica = ca_current(y[7], y[0])
    dy[0] -= 1.7*Ica
    return dy


y0 = array([-40.0, 0.3176, 0.053, 0.596, 1, 0, 0, 0])
r = ode(nmda_HH).set_integrator('dopri5')
r.set_initial_value(y0, -20)
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(ca_current(r.y[7], r.y[0]))
    time.append(r.t)
figure(1, figsize=[10, 6])
plot(time, v)
xlabel('Tiempo (ms)')
ylabel('V (mV)')
grid(True)
fill_between([500, 1500],[120, 120],-40, alpha=0.3)
ylim([-40, 120])
xlim([0, 2000])


Hiperpolarización de la célula

Para simular la hiperpolarización de la célula disminuimos el valor de $V_l$


In [6]:
def nmda_HH(t, y):
    dy = zeros_like(y)
    dy[:4] = HH(y[:4], Vk=-2.0, Vl=-50.0)
    dy[4:] = nmda_impulse(t, y[4:])
    Ica = ca_current(y[7], y[0])
    dy[0] -= 1.7*Ica
    return dy


y0 = array([-50.0, 0.3176, 0.053, 0.596, 1, 0, 0, 0])
r = ode(nmda_HH).set_integrator('dopri5')
r.set_initial_value(y0, -20)
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(ca_current(r.y[7], r.y[0]))
    time.append(r.t)
figure(1, figsize=[10, 6])
plot(time, v)
xlabel('Tiempo (ms)')
ylabel('V (mV)')
grid(True)
fill_between([500, 1500],[120, 120],-60, alpha=0.3)
ylim([-60, 100])
xlim([0, 2000])