hodgkin huxley model of a neuron (part 2)

 

In my last article, I talked about the differential equations that govern the working of a neuron.

 

I'll solve the differential equations to plot a V-t graph and analyse the effect of the input current on the rate of spiking of the neuron.


We will be using numpy and matplotlib modules for this project.
Also we will import division so that we dont end up making any mistake in division

from __future__ import division
from matplotlib import pylab as pl
import numpy as np

First, we will define our functions \(\alpha_n , \beta_n, \alpha_m, \beta_m, \alpha_h, \beta_h\) and also \(n_\infty,m_\infty,h_\infty\)

alpha_n=lambda v: 0.01*(-v+10)/(np.exp((-v+10)/10) - 1) if v!=10 else 0.1
beta_n= lambda v: 0.125*np.exp(-v/80)
n_inf = lambda v: alpha_n(-v)/(alpha_n(-v)+beta_n(-v))

alpha_m=lambda v: 0.1*(-v+25)/(np.exp((-v+25)/10) - 1 ) if v!=25 else 1
beta_m= lambda v: 4*np.exp(-v/18)
m_inf = lambda v: alpha_m(-v)/(alpha_m(-v)+beta_m(-v))

alpha_h=lambda v: 0.07*np.exp(-v/20)
beta_h= lambda v: 1/(np.exp((-v+30)/10)+1)
h_inf = lambda v: alpha_h(-v)/(alpha_h(-v)+beta_h(-v))

Since now we have defined our \(\alpha\) and \(\beta\) functions, we can calculate n,m and h at each instant of time using their differential equations. I will use euler method to solve the differential equations.
since now we can solve for \(n(t) , m(t) , h(t)\), we can calculate \(V(t)\) within a single loop.
But for that we need the values of some constants in our equation of \(\frac{dV}{dt}\) (see part 1). These values are taken from the original paper by Hodgkin and Huxley.
\(g_n=120\)
\(g_k=36\)
\(g_l=0.3\)
\(v_n=115\)
\(v_k=-12\)
\(v_l=10.613\)
\(c=1\)
And we will set up these variables in our program along with lists for V and t. We will initialise n,m,h to their rest values. We will also set the input current signal of 3 different strengths at three different periods (maybe looking at the code will make it clearer) . I think now the following code will make relevance.

g_n=120
g_k=36
g_l=0.3
v_n=115
v_k=-12
v_l=10.613
c=1

dt= 0.025
time=np.arange(0,500,dt)
v=np.zeros(len(time))

v_rest=0
m=m_inf(v_rest)
h=h_inf(v_rest)
n=n_inf(v_rest)
v[0]=v_rest

I_s=np.zeros(len(time))
I_s[1000:5000]=[6.6]*4000
I_s[7000:11000]=[25]*4000
I_s[13000:17000]=[40]*4000

Final part of the program is a loop that will calculate V for each instant of time and put the values in an array that is initialised in the code above (v[0]=v_rest)
Now, we just have to plot V vs t and also I_s (input current) vs t.

for i in range(1,len(time)):
    m+=(alpha_m(v[i-1])*(1-m) - beta_m(v[i-1])*m)*dt
    n+=(alpha_n(v[i-1])*(1-n) - beta_n(v[i-1])*n)*dt
    h+=(alpha_h(v[i-1])*(1-h) - beta_h(v[i-1])*h)*dt
    dv=(1./c)*(I_s[i-1] - g_n*m**3*h*(v[i-1]-v_n) - g_k*n**4*(v[i-1]-v_k) - g_l*(v[i-1]-v_l))*dt
    v[i]=v[i-1]+dv

p=pl.plot(time,v,time,I_s-50)
pl.legend(('voltage spiking','input signal strength'),"upper right")
pl.xlabel("time (ms) --->")
pl.ylabel("potential difference,V (mV) ---->")
pl.show()

The output of the graph is:

 

You can see that the increase in input signal strength increases the rate of spiking.

In the next part of the series I will discuss the phase space of the variables.