Hodgkin Huxley intro

HodgkinHuxleyIPythonNotebook.py

# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>

# <markdowncell>

# Full [Hodgkin Huxley Model](http://en.wikipedia.org/wiki/Hodgkin-Huxley_model)
# ==========================
# 
# By Stephen Larson ([email protected]) for [OpenWorm](http://openworm.org)
# 
# Adapted from UCSD Computational Lab by Jeff Bush ([email protected])

# <codecell>

import scipy as sp
import pylab as plt
from scipy.integrate import odeint
from scipy import stats
import scipy.linalg as lin

# <markdowncell>

# ## Constants
# ### $C_m$ is membrane capacitance, in $\frac{uF}{cm^2}$

# <codecell>

C_m  =   1.0 

# <markdowncell>

# ### $g_{Na}$ is maximum conductance of sodium, in $\frac{mS}{cm^2}$

# <codecell>

g_Na = 120.0 

# <markdowncell>

# ### $g_{K}$ is maximum conductance of potassium, in $\frac{mS}{cm^2}$

# <codecell>

g_K  =  36.0

# <markdowncell>

# ### $g_{L}$ is maximum conductance of leak, in $\frac{mS}{cm^2}$

# <codecell>

g_L  =   0.3

# <markdowncell>

# ### The Nernst reversal potentials of sodium $E_{Na}$, potassium $E_K$ and leak $E_L$, in millivolts $mV$

# <codecell>

E_Na =  50.0 
E_K  = -77.0
E_L  = -54.387

# <markdowncell>

# # Channel gating kinetics
# ## Activation variables are functions of membrane voltage
# 

# <markdowncell>

# $\alpha_m(V)\ = 0.1\ \frac{V + 40}{ 1 - \exp(-\frac{V + 40}{10})}$
# 
# $$\beta_m(V)\ = 4\ \exp(-\frac{V + 65}{18}) $$
# 
# $$\alpha_h(V)\ = 0.07\ \exp(-\frac{V + 65}{20}) $$
# 
# $$\beta_h(V)\ = \frac{1}{1 + \exp(-\frac{V + 35}{10}) }$$
# 
# $$\alpha_n(V)\ = 0.01 \ \frac{V + 55}{ 1 - \exp(-\frac{V + 55}{10})} $$
# 
# $$\beta_n(V)\ = 0.125\ \exp(-\frac{V + 65}{ 80})$$

# <codecell>

def alpha_m(V): return 0.1*(V+40.0)/(1.0 - sp.exp(-(V+40.0) / 10.0))
def beta_m(V):  return 4.0*sp.exp(-(V+65.0) / 18.0)
def alpha_h(V): return 0.07*sp.exp(-(V+65.0) / 20.0)
def beta_h(V):  return 1.0/(1.0 + sp.exp(-(V+35.0) / 10.0))
def alpha_n(V): return 0.01*(V+55.0)/(1.0 - sp.exp(-(V+55.0) / 10.0))
def beta_n(V):  return 0.125*sp.exp(-(V+65) / 80.0)

# <markdowncell>

# # Membrane currents (in $\frac{uA}{cm^2}$)
# 
# ## Sodium (Na = element name)
# $$I_{Na} = g_{Na}\ m^3 h\ (V - E_{Na})$$

# <codecell>

def I_Na(V,m,h):return g_Na * m**3 * h * (V - E_Na)

# <markdowncell>

# ##  Potassium (K = element name)
# $$I_K = g_K\ n^4\ (V - E_K)$$

# <codecell>

def I_K(V, n):  return g_K  * n**4     * (V - E_K)

# <markdowncell>

# ##Leak
# $$I_L = g_L\ (V - E_L)$$

# <codecell>

def I_L(V):     return g_L             * (V - E_L)

# <markdowncell>

# # External current
# ## step up 10 $\frac{uA}{cm^2}$ after $10ms$

# <codecell>

def I_inj(t): 
    return 10*(t>10) - 10*(t>100)
    #return 10*(t>100) - 10*(t>200) + 35*(t>300)
    #return 10*t

# <markdowncell>

# # The time to integrate over

# <codecell>

t = sp.arange(0.0, 100.0, 0.1)

# <markdowncell>

# # Integrate!
# $C_m\ {dV \over dt} = -I_{Na} -I_K - I_L + I_{ext}$

# <codecell>

def dALLdt(X, t):
    V, m, h, n = X
    
    #calculate membrane potential & activation variables
    dVdt = (I_inj(t) - I_Na(V, m, h) - I_K(V, n) - I_L(V)) / C_m
    dmdt = alpha_m(V)*(1.0-m) - beta_m(V)*m
    dhdt = alpha_h(V)*(1.0-h) - beta_h(V)*h
    dndt = alpha_n(V)*(1.0-n) - beta_n(V)*n
    return dVdt, dmdt, dhdt, dndt
    
X = odeint(dALLdt, [-65, 0.05, 0.6, 0.32], t)

# <markdowncell>

# # Retrieve the values from the output matrix, X

# <codecell>

V = X[:,0]
m = X[:,1]
h = X[:,2]
n = X[:,3]
ina = I_Na(V,m,h)
ik = I_K(V, n)
il = I_L(V)

# <markdowncell>

# # Plot!

# <markdowncell>

# ## Injected current: $I_{ext}$

# <codecell>

plt.title("Injected current: $I_{ext}$")
plt.plot(t, I_inj(t), 'k')
plt.xlabel('t (ms)')
plt.ylabel('$I_{inj}$ ($\\mu{A}/cm^2$)')
plt.ylim(-1, 31)

plt.show()

# <markdowncell>

# Hodgkin-Huxley Neuron membrane potential: $V$

# <codecell>

plt.title('Hodgkin-Huxley Neuron membrane potential: $V$')
plt.plot(t, V, 'k')
plt.ylabel('V (mV)')
plt.xlabel('t (ms)')

plt.show()

# <markdowncell>

# ## Currents: $I_{Na}$, $I_{K}$, $I_{L}$

# <codecell>

plt.title("Currents: $I_{Na}$, $I_{K}$, $I_{L}$")
plt.plot(t, ina, 'c', label='$I_{Na}$')
plt.plot(t, ik, 'y', label='$I_{K}$')
plt.plot(t, il, 'm', label='$I_{L}$')
plt.ylabel('Current')
plt.xlabel('t (ms)')
plt.legend()

plt.show()

# <markdowncell>

# ## Activation variables: m, n, h

# <codecell>

plt.title("Activation variables: m, n, h")
plt.plot(t, m, 'r', label='m')
plt.plot(t, h, 'g', label='h')
plt.plot(t, n, 'b', label='n')
plt.ylabel('Gating Value (unitless)')
plt.xlabel('t (ms)')
plt.legend()

plt.show()