Previous Contents Next

Chapter 2  The Basic Hodgkin - Huxley Model:

The standard Hodgkin - Huxley model of an excitatory neuron consists of the equation for the total membrane current, IM, obtained from Ohm's law:
IM =
Cm
dV
dt
+ IK + INa + IL  
    (2.1)

where V denotes the membrane voltage, IK is the potassium current, INa is the sodium current and IL is leakage current carried by other ions that move passively through the membrane. This equation is derived by modeling the potassium, sodium and leakage currents using a simple electrical circuit model of the membrane. We think of a gate in the membrane as having an intrinsic resistance and the cell membrane itself as having an intrincis capacitance as shown in Figure 2.1:

Figure 2.1: The Membrane and Gate Circuit Model

Here we show an idealized cell with a small portion of the membrane blown up into an idealized circuit. We see a small piece of the lipid membrane with an inserted gate. We think of the gate as having some intrinsic resistance and capacitance. Now for our simple Hodgkin - Huxley model here, we want to model a sodium and potassium gate as well as the cell capacitance. So we will have a resistance for both the sodium and potassium. In addition, we know that other ions move across the membrane due to pumps, other gates and so forth. We will temporarily model this additional ion current as a leakage current with its own resistance. We also know that each ion has its own equilibrium potential which is determined by applying the Nernst equation. The driving electomotive force or driving emf is the difference between the ion equilibrium potential and the voltage across the membrane itself. Hence, if Ec is the equilibrium potential due to ion c and Vm is the membrane potential, the driving force is Vc - Vm. In Figure 2.2, we see an electric schematic that summarizes what we have just said. We model the membrane as a parellel circuit with a branch for the sodium and potassium ion, a branch for the leakage current and a branch for the membrane capacitance.

Figure 2.2: The Simple Hodgkin - Huxley Membrane Circuit Model

From circuit theory, we know that the charge q across a capacitator is q = C E, where C is the capacitance and E is the voltage across the capicitor. Hence, if the capacitance C is a constant, we see that the current through the capacitor is given by the time rate of change of the charge
dq
dt
=
C
dE
dt

If the voltage E was also space dependent, then we would write E(z,t) to indicate its dependence on both a space variable z and the time t. Then the capacitive current would be
dq
dt
=
C
E
t

From Ohm's law, we know that voltage is current times resistance; hence for each ion c, we can say
Vc = Ic Rc

where we label the voltage, current and resistance due to this ion with the subscript c. This implies
Ic =
1
Rc
Vc
  = gc Ic

where gc is the reciprocal resistance or conductance of ion c. Hence, we can model all of our ionic currents using a conductance equation of the form above. Of course, the potassium and sodium conductances are nonlinear functions of the membrane voltage V and time t. This reflects the fact that the amount of current that flows through the membrane for these ions is dependent on the voltage differential across the membrane which in turn is also time dependent. The general functional form for an ion c is thus
Ic = gc(V,t) ( V(t) - Ec(t) )

where as we mentioned previously, the driving force, V - Ec, is the difference between the voltage across the membrane and the equilibrium value for the ion in question, Ec. Note, the ion battery voltage Ec itself might also change in time (for example, extracellular potassium concentration changes over time ). Hence, the driving force is time dependent. The conductance is modeled as the product of a activation, m, and an inactivation, h, term that are essentially sigmoid nonlinearities. The activation and inactivationa are functions of V and t also. The conductance is assumed to have the form
gc(V,t) = g0 mp(V,t) hq(V,t)

where appropriate powers of p and q are found to match known data for a given ion conductance.

We model the leakage current, IL, as
IL = gL ( V(t) - EL )

where the leakage battery voltage, EL, and the conductance gL are constants that are data driven.

Hence, our full model would be
Im =
Cm
dV
dt
+ gK (V - EK) + gNa (V - ENa) + gL ( V - EL )

2.1  Activation and Inactivation Variables:

We assume that the voltage dependence of our activation and inactivation has been fitted from data. Hodgkin and Huxley modeled the time dependence of these variables using first order kinetics. They assumed a typical variable of this type, say m, satisfies for each value of voltage, V:
dx(V)
dt
= ax(V) (1 - x(V)) + bx(V) x(V)
x(V,0) = x0(V)

For convenience of exposition, we usually drop the functional dependence of x on V and just write:
dx
dt
= ax (1 - x) + bx x
x(0) = x0

Rewriting, we see
1
ax + bx
dx
dt
=
ax
ax + bx
- x
x(0) = x0

We let
tx =
1
ax + bx
x¥ =
ax
ax + bx

allowing us to rewrite our rate equation as
tx
dx
dt
= x¥ - x
x(0) = x0

This first order equation is easily solved to give
x(t) =
(x0 - x¥) e
-t
tx
 
+ x¥

2.2  The Hodgkin-Huxley Sodium and Potassium Model:

Hodgkin and Huxley modeled the sodium and potassium gates as
gNa(V,t) = g0Na m3(V,t) h(V,t)
gK(V,t) = g0K n4(V,t)

where the two activation variables, m and n, and the one inactivation variable all satisfy first order kinetics as we have discussed. Hence, if the parameters tm and m¥ and so forth were constants, we would know
m(t) =
(m0 - m¥) e
-t
tm
 
) + m¥
h(t) =
(h0 - h¥) e
-t
th
 
) + h¥
n(t) =
(n0 - n¥) e
-t
tn
 
) + n¥

with
tm =
1
am + bm
m¥ =
am
am + bm
th =
1
ah + bh
h¥ =
ah
ah + bh
tn =
1
an + bn
n¥ =
an
an + bn

However, all of the coefficient functions, a and b for each variable required data fits as functions of voltage. These were determined to be
am =
-0.10
V + 35.0
e -.1 (V +35.0) - 1.0
bm =
4.0 e
-(V+60.0
18.0
 
ah = 0.07 e-.05 (V+60.0)
bh =
1.0
(1.0 + e-.1 (V+30.0))
an =
-
0.01*(V+50.0
(e -.1 (V+50.0) - 1.0)
bn = 0.125 e-0.0125 (V+60.0

Of course these data fits were obtained at a certain temperature and assumed values for all the other constants needed. These other parameters are given in the units below
voltage mV milli volts 10-3 Volts
current na nano amps 10-9 Amps
time ms milli seconds 10-3 Seconds
concentration mM milli moles 10-3 Moles
conductance µ S micro Siemens 10-6 ohms-1
capacitance nF nano Fahrads 10-9 Fahrads

Our model of the membrane dynamics here thus consists of the following differential equations:
tm
dm
dt
= m¥ - m
th
dh
dt
= h¥ - h
tn
dn
dt
= n¥ - n
dV
dt
=
IM - IK - INa - IL
CM

with initial conditions
m(0) = m¥(V0,0)
h(0) = h¥(V0,0)
n(0) = m¥(V0,0)
V(0) = V0

We note that at equilibrium there is no current across the membrane. Hence, the sodium and potassium currents are zero and the activation and inactivation variables should achieve their steady state values which would be m¥, h¥ and n¥ computed at the equilibrium membrane potential which is here denoted by V0.

2.2.1  Encoding The Dynamics:

Now these dynamics are more difficult to solve than you might think. The sequence of steps is this: The basic Hodgkin - Huxley model thus needs four independent variables which we place in a four dimensional vector y which components assigned as follows:
y =
é
ê
ê
ê
ë
y[0] = V
y[1] = m
y[2] = h
y[3] = n
ù
ú
ú
ú
û

We encode the dynamics calculated above into a four dimensional vector f whose components are interpreted as follows:
f = =
é
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ë
f[0]
=
IM - IT
CM
f[1]
=
m¥ - m
tm
f[2]
=
h¥ - h
th
f[3]
=
n¥ - n
tn
ù
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
û

which in terms of our vector y becomes
é
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ë
f[1]
=
m¥ - y[1]
tm
f[2]
=
h¥ - y[2]
th
f[3]
=
n¥ - y[3]
tn
ù
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
û

So to summarize, at each time t and V, we need to calculate the four dimensional dynamics vector f for use in our choice of ODE solver.

2.2.2  Computing the Solution Numerically:

As an example of what we might do to solve this kind of a problem, we note that this system can now be written in vector form as
dy
dt
= f(t,y)
y(0) = y0

where we have found that
y =
é
ê
ê
ê
ë
V
m
h
n
ù
ú
ú
ú
û		 ,

f =
é
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ë
IM - IT
CM
m¥ - m
tm
h¥ - h
th
n¥ - n
tn
ù
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
û

with initial data
y0 =
é
ê
ê
ê
ë
V0
m¥
h¥
n¥
ù
ú
ú
ú
û

The hardest part is to encode all of this in the vector form. The dynamic component f[0] is actually a complicated function of m, h and n as determined by the encoding sequence we have previously discussed. Hence, the calculation of the dynamics vector f is not as simple as defining each component of f as a functional expression.

Estimating The Solution Numerically:

If we want to solve the scalar equation (this means x is not a vector)
dx
dt
= f(t,x)
x(t0) = x0

one way we can attack it is to assume the solution x has enough smoothness to allow us to expand it in a Taylor series expansion about the base point (t0,x0). Using the Taylor Remainder theorem, this gives for a given time point t:
x(t) =
x(t0) +
dx
dt
(t0) (t - t0) +
1
2
d2x
dt2
(t0) (t - t0)2 +
1
6
d3x
dt3
(ct) (t - t0)3

where ct is some point in the interval [t0,t]. We also know by the chain rule that since x' is f that
d2x
dt2
=
f
t
+
f
t
dx
dt
d2x
dt2
=
f
t
+
f
x
f

which implies, switching to a standard subscript notation for partial derivatives (yes, these calculations are indeed yucky2!)
d3x
dt3
= (ftt + ftx f) f + (fxt + fxx f) f + fx (ft + fx f)

The important thing to note is that this third order derivative is made up of algebraic combinations on f and its various partial derivatives. We typically assume that on the interval [t0,t], that all of these functions are continous and bounded. Thus, letting ||g|| represent the maximum value of the continous function g(s,u) on the interval [t0,t] × [x0,x1] for suitable x1, we know there is a constant C so that
|
d3x
dt3
(ct)|
£ (||ftt|| + ||ftx|| ||f||) ||f|| + (||fxt|| + ||fxx|| ||f||) ||f|| + ||fx|| (||ft|| + ||fx|| ||f||)
= C

Now using the standard abbreviations
f0 = f(t0,x))
ft0 =
f
t
(t0,x0)
fx0 =
f
x
(t0,x0)
ftt0 =
2 f
t2
(t0,x(t0))
ftx0 =
2 f
x t
(t0,x(t0))
fxt0 =
2 f
t t
(t0,x(t0))
fxx0 =
2 f
x2
(t0,x(t0))

we see our solution can be written as
x(t) =
x0 + f0 (t - t0) +
1
2
(ft0 + fx0 f0) (t - t0)2 +
1
6
d3x
dt3
(ct) (t - t0)3

Some Detailed Error Estimates:

Now let's use the results we just derived as follows: let's choose a time step h and set t to be t0 + h. Then we have (t - t0) is h and our solution is
x(t0+h) =
x0 + f0 h +
1
2
(ft0 + fx0 f0) h2 +
1
6
d3x
dt3
(ct) h3

Now, note that using a Taylor expansion for the function f, we can write for a base point (t0,x0)
f(t0+
1
2
h,x0 +
1
2
h f0)
=
f0 + ft0
1
2
h + fx0
1
2
h f0 +
1
2
H0(d,u)

where H0(d,u) is what is called a Hessian term and is defined by
H0(d,u) =
é
ê
ê
ë
1
2
h
1
2
hf0
ù
ú
ú
û		 é
ë
ftt(d,u) ftx(d,u)
fxt(d,u) fxx(d,u)
ù
û		 é
ê
ê
ê
ê
ê
ê
ë
1
2
h
1
2
h f0
ù
ú
ú
ú
ú
ú
ú
û

The point (d,u) satisfies d is in the interval [t0,t0+ 1/2 h] and u is in the interval [x0,x0+ 1/2 h f0]. Note that we can estimate the effect of this term by
1
2
||H0||
£
1
2
( ||ftt||
h2
4
+ ||ftx||
h2
4
|f0| + ||fxx||
h2
4
|f0|2
  =
h3
8
( ||ftt|| + ||ftx|| |f0| + ||fxx|| |f0|2 )

Thus for suitable constant D we can say
1
2
||H0||
£
h3
8
D

From the above, we note we can write
x(t0+h) =
x0 + h f(t0+
1
2
h,x0 +
1
2
h f0) + H0(d,u) +
1
6
d3x
dt3
(ct) h3

Finally, we see that if we take our maximum function values on the domain [t0,t0+h] × [x0, x0 + h f0] for all the terms in d3x/dt3, then there is a constant C so that
1
6
||
d3x
dt3
||
£
h3
6
C

and so the error in approximating the true solution x by the term x0 + h f(t0+ 1/2h,x0 + 1/2 h f0) can be written as e(t) where on [t0,t0+h]
|e(t)| £
h3
8
D +
h3
6
C

Thus, we can build an approximate solution f to the true solution x that uses the information we know at t0 with data x0 to construct the approximation
f(t) =
x0 + h f(t0+
1
2
h,x0 +
1
2
h f0)

which has a known error proportional to h3. Since this error is due to the fact that we truncate the Taylor series expansions for f earlier, this is known as truncation error. We can follow the process above to build other sorts of approximations that use more terms of the Taylor Series expansions. We would call this approximation an order 3 method because the local error is proportional to the third power of the step size h. Note our simple approximation If we carried out our Taylor series approximations one more term, that is we use remainder terms due to d4 x/dt4, we would find our local error is proportional to h4 and we would need three function values to build our approximation. In general, to obtain a local order hp method, we need p-1 function evaluations if we use this technique. The example we have worked out here is called the Runga - Kutte Order 2 Method. Now to solve a problem on a long time interval, we do this approximation for each step and consecutively build an approximation solution sequence x0, x(t0+h1), x(t0+h1+h2) and so forth where hi denotes the step size we use at the ith step.

Finally, if we solve a vector system, we basically apply the technique above to each component of the vector solution. So if the problem is four dimensional, we are gettting a constant C and D for each component. This of course complicates our life as one component might have such a fast growing error that the C and D for it is very large compared to the others. To control this, we generally cut the step size h to bring down the error. Hence, the component whose error grows the fastest controls how big a step size we can use in this type of technique.

Runge-Kutta Methods:

As discussed above, these methods are based on multiple function evaluations. Say you want to solve our problem on the interval [0,tf]. If we divide this interval into a collection of discrete points somehow, say { tn } for n ranging from 0 to a largest index N, then at each such time point, we generate a new estimate of the vector solution yn. So we start at the initial data y0 at time 0 and use an estimation formula to generate a new approximate solution value y1 for time t1 and so forth until we hit the final time tN with associated final approximate solution value yN.

If we want a solution on the interval [0,100] for example, if we used a uniform step size of h = 1.0e-3, we would generate 1000 y vectors for each unit of time. Hence, our solution would be encoded in a vector of size 105 and yn would be the numerical solution to our problem for time tn = n h.

Now in the theory of the numerical solution of ordinary differential equations (ODEs) we find that there a several sources of error in this process:
  1. As discussed earlier, usually we approximate the function f using a few terms of the Taylor series expansion of f around the point (tn,yn) at each iteration. This truncation of the Taylor series expansion is not the same as the true function f of course, so there is an error made. Depending on how many terms we use in the Taylor series expansion, the error can be large or small. If we let hn denote the difference tn+1 - tn, a fourth order method is one where this error is of the from C hn5 for some constant C. This means that is you use a step size which is one half hn, the error decreases by a factor of 32. For a Runga - Kutte method to have a local truncation error of order 5, we would need to do 4 function evaluations. Now this error is local to this time step, so if we solve over a very long interval with say N being 100,000 or more, the global error can grow quite large due to the addition of so many small errors. So the numerical solution of an ODE can be very acurrate locally but still have a lot of problems when we try to solve over a long time interval. For example, if we want to track a voltage pulse over 1000 milliseconds, if we use a step size of 10-4 this amounts to 107 individual steps. Even a fifth order method can begin to have problems. This error is called truncation error.
  2. We can't represent numbers with perfect precision on any computer system, so there is an error due to that which is called round - off error which is significant over long computations.
  3. There is usually a modeling error also. The model we use does not perfectly represent the physics, biology and so forth of our problem and so this also introduces a mismatch between our solution and the reality we might measure in a laboratory.
The solution scheme discussed above generates a sequence { yn } starting at y0 using the some sort of recursion equation:
yn+1 = yn + hn × F(tn,yn,hn,f)
y0 = y0

where hn is the step size we use at time point tn. and the function F is some function of the previous approximate solution yn, the step size hn and the dynamics vector f. We usually choose a numerical method which allows us to estimate what a good step size would be at each step n and then alter the step size to that optimal choice. Our order two method discussed earlier used
F(tn,yn,hn,f) =
f(tn +
1
2
hn, xn +
1
2
hn f(tn,xn))

A typical numerical algorithm which will do this is the Runge - Kutta - Fehlberg 45, RKF45 algorithm which works like this: Since this algorithm uses siz function evaluations per time step that is actually used, it is very costly. Indeed, if we alter h, we go back and recompute all of the solutions. There are better methods such as Adams - Bashford and Adams - Moulton. However, for our tutorial here, we will start with this method.

Finally, just to give you the flavor of what needs to be computed, here is a outline of the method: in what follows K is a six dimensional vector which we use to store intermediate results. Now recall that the vector function f we use in our Hodgkin - Huxley simulations is forbiddingly complex and so these six function evaluations really slow us down!!

2.3  The Koch Modified Sodium/ Potassium Model:

We will not use these data fits. Instead we will use a more complicated model given in (Koch, 8) which models not only Sodium and Potassium voltage dependent gates but also Calcium gates, extracellular Potassium concentration outside the cellular membrane, Calcium pumps and even the diffusion of Calcium ions across the cellular interior. You will probably want to look at this article to get a better appreciation of how complex it is! We will start with only Sodium, Potassium gates and leakage current following Equation2.1. This means we will not use many of the other gates and so forth in the Koch model. Our simple model then consists of the transient outward potassium current and one type of potassium current called the transient outward Potassium Current. Since there will eventually be three potassium currents, this particular current will modeled using variables with K1 subscripts while as usual the sodium current will be labeled with variables using NA subscripts. In addition, the Koch model uses a somewhat diferent gate model from that of Hodgkin and Huxley:
gNa(V,t) = g0Na mNA2(V,t) hNA(V,t)
gK1(V,t) = g0K1 mK1(V,t) hK1(V,t)

where you will note that the Potassium gate now has an inactivation component. Our previous discussion for the kinetics of the activation and inactivation still applies though so we know these variables satisfy
mNA(V,t) =
(m0 - mNA¥(V)) e
-t
tNAm(V)
 
+ mNA¥(V)

where the initial constant m0 should be chosen as follows: we assume that the cell is at equilibrium voltage EM. Then the steady state mNA value is obtained by letting t go to infinity: hence we should choose
m0 = mNA¥(EM)

Hence, we have
mNA(V,t) =
(mNA¥(EM) - mNA¥(V)) e
-t
tNAm(V)
 
+ mNA¥(V)

A similar analysis shows
hNA(V,t) =
(hNA(EM) - hNA¥(V)) e
-t
tNAh(V)
 
+ hNA¥(V)
mK1(V,t) =
(mK1(EM) - mK1¥(V)) e
-t
tK1m(V)
 
+ mK1¥
hK1(V,t) =
(hK1¥ - hK1¥(V)) e
-t
tK1h(V)
 
+ hK1¥

The important steady state parameters ( m¥ types) and time constants ( t types) are then typically given in terms of the voltage dependent pairs of constants (a,b) as previously discussed. Now these functions must be provided from curve fit data, spline or neural net models and so forth. Here we will provide nonlinear curve fits as given by (Koch, 8). The t parameters in (Koch, 8) are doubled from what would follow from our discussion of the kinetics; hence there is a factor of two in the t equations that follow:
tNAm =
2
aNAm + bNAm
mNA¥ =
aNAm
aNAm + bNAm
tNAh =
2
aNAh + bNAh
hNA¥ =
aNAh
aNAh + bNAh

The a and b curve fits are given to us as follows:
aNAm =
0.36 (V+33.0)
1.0 - e
-(V+33.0)
3.0
 
bNAm =
-0.40 (V+42.0)
1.0 - e
V+42.0
20.0
 
aNAh =
-0.10 (V+55.0)
1.0 - e
V+55.0
6.0
 
bNAh =
4.5
1.0 + e
-
V
10.0
 

In Figure 2.3 we can see the a and beta curves for the activation variable mNA plotted together.
<