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.

Figure 2.3: Activation a and b Curves for Sodium:

Further, we can the the a and beta curves plotted together for the inactivation variable hNA in Figure 2.4

Figure 2.4: Inactivation a and b Curves for Sodium:

The (a,b) for the activation and inactivation variables for the h parameters of the K1 gate are not provided in (Koch, 8); instead the mK1¥, tK1m, hK1¥ and tK1h are given to us directly.
mK1¥ =
1
1.0 + e
-
V+42.0
13.0
 
tK1m = 1.38
mK1h =
1
1.0 + e
V+110.0
18.0
 
tK1h =
ì
í
î
50 V < -80
150 V ³ -80

If the voltage was a constant, we can actually integrate these kinetics in closed form obtaining
gNA(V,t) = gNA0 mNA2(V,t) hNA(V,t)
gK1(V,t) = gK10 mK1(V,t) hK1(V,t)

Of course, the voltage is not constant, so we treat mNA, hNA and so forth as independent variables in our ODE integration scheme. We need to have five independent variables here so we set up a five dimensional vector y whose components will be interpeted as follows:
y[0] = V
y[1] = mNA
y[2] = hNA
y[3] = mK1
y[4] = hK1

where we use as starting conditions
T = 22.0
gL = 0.02
EL = -49.0
CM = 0.15
gNa0 = 2.0
gK10 = 2.7708

with initial concentrations
[NA]o = 491.0
[NA]i = 50.0
[K]o = 7.859
[K]i = 140.0

The Koch model has many additional gates and so forth with associated conductances. For this simple one gate Sodium and one gate Potassium model we will lump all these other conductances into the K1 term giving the number 2.7708 obtained from
g0K1 = 0.120 + 1.17 + 0.084 + 1.2 + 0.054 + 0.02675 + 0.116
  = 2.7708

Note that we can clearly see that the Sodium conductance activates early and the Potassium conductance lags from the three plots in Figure 2.5 (gNA), Figure 2.6 (gK1) and Figure 2.7 (a surface plot of .10 gNA and gK1 together).

Figure 2.5: Sodium Conductance:


Figure 2.6: Potassium Conductance:


Figure 2.7: Sodium and Potassium Conductance:

2.4  Hodgkin-Huxley Simulation Code:

We will use both Matlab and C++ to develop simulations of a working cell.

2.4.1  C++ Simulation Code: Hodgkin - Huxley

We can use our standard ODE X Windows simulation package (see the Chapter 6 ) to develop a working simulation which we will be able to modify to handle calcium dynamics. We begin with a simple Hodgkin - Huxley giant squid axon simulation. The system to solve is encoded in the following function:

function.c:

The giant squid axon system is encode here for use in a standard Runga - Kutta - Fehlberg Order 5 method. We only use simple sodium and potassium gates here.

function.c 
  1 #include "simulation.h"
  2 
  3 /*-------------------------------------------------------
  4   This function models simple HH dynamics: Na and K gates
  5   only.  
  6   -------------------------------------------------------*/
  7 
  8 void funcp(double t,DOUBLE_VECTOR& y,DOUBLE_VECTOR& f,
  9            DOUBLE_VECTOR& p)
 10 {
 11   double g_NA_bar    = p[0];
 12   double g_K1_bar    = p[1];
 13   double g_leak_bar  = p[2];
 14   double E_NA        = p[3];
 15   double E_K         = p[4];
 16   double E_L         = p[5];
 17     
 18   double sum;
 19   double I_NA,I_K1;
 20   //units
 21   /*===============================================
 22     Methods in Neuronal Modeling: Chapter 4: 
 23       Dynamic Channels by 
 24       Yamada, Koch and Adams
 25     Their Units:
 26     voltage       mV
 27     current       na
 28     time          ms
 29     concentration mM
 30     conductance   micro Siemens
 31     capacitance   nF
 32     ===============================================*/
 33   // y vector assignments
 34   /*===============================================
 35     y[0] = V
 36     y[1] = m_NA
 37     y[2] = h_NA
 38     y[3] = m_K1
 39     y[4] = h_K1
 40     ===============================================*/
 41   /*===============================================
 42     Fast Sodium Current
 43     ===============================================*/
 44 
 45   // activation parameter for I_NA
 46   double alpha_mNA;
 47   double beta_mNA;
 48   double m_NA_infinity;
 49   double t_m_NA;
 50 
 51   sum = y[0]+33.0;
 52   if(sum>0){
 53     alpha_mNA     = 0.36*sum/(1.0 - exp(-sum/3.0));
 54     }
 55   else{
 56     alpha_mNA     = 0.36*exp( sum/ 3.0)*sum/(exp(sum/ 3.0) - 1.0);
 57     }
 58   sum = y[0]+42.0;
 59   if(sum>0){
 60     beta_mNA      = -0.4*exp( -sum/20.0)*sum/
 61                      (exp( -sum/20.0) - 1.0);
 62     }
 63   else{
 64     beta_mNA      = -0.4*sum/
 65                      (1.0 - exp(sum/20.0));
 66     }
 67   m_NA_infinity = alpha_mNA/(alpha_mNA+beta_mNA);
 68   t_m_NA        = 2.0/(alpha_mNA+beta_mNA);
 69   f[1]          = (m_NA_infinity - y[1])/t_m_NA;
 70 
 71   // inactivation parameter for I_NA
 72   double alpha_hNA;
 73   double beta_hNA;
 74   double h_NA_infinity;
 75   double t_h_NA;
 76 
 77   sum = y[0]+55.0;
 78   if(sum<0){
 79     alpha_hNA     = -0.1*sum/
 80                      (1.0 - exp(sum/ 6.0));
 81     }
 82   else{
 83     alpha_hNA     = -0.1*exp( -sum/ 6.0)*sum/
 84                      (exp( -sum/ 6.0) - 1.0);
 85     }
 86   if(y[0]>0){
 87     beta_hNA      =  4.5/(1.0 + exp(-y[0]/10.0));
 88     }
 89   else{
 90     beta_hNA      =  4.5*exp(y[0]/10.0)/(exp(y[0]/10.0) + 1.0);
 91     }
 92   h_NA_infinity = alpha_hNA/(alpha_hNA+beta_hNA);
 93   t_h_NA        = 2.0/(alpha_hNA+beta_hNA);
 94   f[2]          = (h_NA_infinity - y[2])/t_h_NA; 
 95 
 96   I_NA = g_NA_bar*(y[0]-E_NA)*y[1]*y[1]*y[2]; 
 97 
 98   /*===============================================
 99     Transient, Outward Potassium Current
100     ===============================================*/
101   // activation parameter for I_K
102   double t_m_K1;
103   double m_K1_infinity;
104 
105   sum = y[0]+42.0;
106   if(sum<0){
107     m_K1_infinity = exp(sum/13.0)/(exp(sum/13.0)+1.0);
108     }
109   else{
110     m_K1_infinity = 1.0/(1.0+exp(-sum/13.0));
111     }    
112   t_m_K1        = 1.38;
113   f[3]          = (m_K1_infinity - y[3])/t_m_K1;
114 
115   // inactivation parameter for I_K
116   double h_K1_infinity;
117   double t_h_K1;
118 
119   sum = y[0]+110.0;
120   if(sum<0){
121     h_K1_infinity = 1.0/(1.0+exp( sum/18.0));
122     }
123   else{
124     h_K1_infinity = exp( -sum/18.0)/(exp(-sum/18.0)+1.0);
125     }
126   if(y[0]<-80.0)
127     t_h_K1 = 50.0;
128   else        
129     t_h_K1 = 150.0;
130   f[4]          = (h_K1_infinity - y[4])/t_h_K1; 
131 
132   I_K1 = g_K1_bar*(y[0]-E_K)*y[3]*y[4];
133 
134   /*========================================================
135     leakage current
136     ========================================================*/  
137   double I_leak = g_leak_bar*(y[0]-E_L);
138 
139   //Cell Voltage equation
140   static double C_N = .15;
141   double external_current = 0.0;
142 
143   static double delta_t = 1.0;
144   static double t_start = 10.0;
145   static double t_end = t_start + delta_t;
146   static double injection = 20.0;
147   
148   if(t>= t_start && t<= t_end)
149         external_current = injection; 
150   
151 /* 
152   for(int burst=2;burst<51;burst+=50){
153     for(int cnt=0;cnt<1;cnt+=1){
154       sum = burst+.05*cnt;
155       if(t>= sum && t<= sum+delta_t){
156         external_current = injection;
157         }
158       }
159     }      
160 */
161   
162   f[0] = (external_current- I_NA - I_K1 - I_leak)/C_N ;
163 }
syntax highlighted by Code2HTML, v. 0.8.8b

rest.c:

The initial conditions for the activationa and inactivation variables for our gates are a little tricky to caclulate as we discussed eaqrlier. We do this is a startup function rest.c:

rest.c 
  1 #include "simulation.h"
  2 
  3 double rest(double E_M,DOUBLE_VECTOR& p,DOUBLE_VECTOR& q)
  4 {
  5   double sum;
  6   
  7   double g_NA_bar    = p[0];
  8   double g_K1_bar    = p[1];
  9   double g_leak_bar  = p[2];
 10   double E_NA        = p[3];
 11   double E_K         = p[4];
 12   //double E_L         = p[5];
 13   
 14   cout << "In rest calculation: parameters are p = " << p << endl;
 15   
 16   /*----------------------------------------------------------
 17     How do we choose E_L and g_L?
 18     
 19     At equilibrium we must have:
 20     g_total = g_NA_bar+g_K1_bar;
 21     ratio_NA = g_NA_bar/(g_total+g_leak_bar);
 22     ratio_K1 = g_K1_bar/(g_total+g_leak_bar);
 23     ratio_L = g_leak_bar/(g_total+g_leak_bar);
 24     E_M = ratio_NA*E_NA + ratio_K1*E_K + ratio_L*E_L;
 25     
 26     For E_M we must also have
 27 
 28     I_NA = g_NA_bar*(E_M-E_NA)*m_NA_infinity*m_NA_infinity
 29             *h_NA_infinity;    
 30     I_K1 = g_K1_bar*(E_M-E_K)*m_K1_infinity*h_K1_infinity;
 31     sum = I_NA + I_K1; 
 32     g_leak_bar  = -sum/(E_M-E_L);
 33     
 34     So
 35     E_M*(g_total+g_leak_bar) = g_K*E_NA + g_K1*E_K1 + g_leak_bar*E_L;
 36     
 37     Let 
 38     
 39     S1 = g_K*E_NA + g_K1*E_K1 ;
 40     S2 = sum;
 41     
 42     Then
 43     
 44     g_leak_bar = -S2/(E_M-E_L);
 45     
 46     yielding
 47     
 48     E_M*(g_total -S2/(E_M-E_L)) = S1 + (-S2/(E_M-E_L)*E_L;
 49     ----------------------------------------------------------*/
 50 
 51   // activation parameter for I_NA
 52   double alpha_mNA;
 53   double beta_mNA;
 54   double m_NA_infinity;
 55 
 56   sum = E_M+33.0;
 57   if(sum>0){
 58     alpha_mNA     = 0.36*sum/(1.0 - exp(-sum/3.0));
 59     }
 60   else{
 61     alpha_mNA     = 0.36*exp( sum/ 3.0)*sum/(exp(sum/ 3.0) - 1.0);
 62     }
 63   sum = E_M+42.0;
 64   if(sum>0){
 65     beta_mNA      = -0.4*exp( -sum/20.0)*sum/
 66                      (exp( -sum/20.0) - 1.0);
 67     }
 68   else{
 69     beta_mNA      = -0.4*sum/
 70                      (1.0 - exp(sum/20.0));
 71     }
 72   m_NA_infinity = alpha_mNA/(alpha_mNA+beta_mNA);
 73 
 74   // inactivation parameter for I_NA
 75   double alpha_hNA;
 76   double beta_hNA;
 77   double h_NA_infinity;
 78 
 79   sum = E_M+55.0;
 80   if(sum<0){
 81     alpha_hNA     = -0.1*sum/
 82                      (1.0 - exp(sum/ 6.0));
 83     }
 84   else{
 85     alpha_hNA     = -0.1*exp( -sum/ 6.0)*sum/
 86                      (exp( -sum/ 6.0) - 1.0);
 87     }
 88   if(E_M>0){
 89     beta_hNA      =  4.5/(1.0 + exp(-E_M/10.0));
 90     }
 91   else{
 92     beta_hNA      =  4.5*exp(E_M/10.0)/(exp(E_M/10.0) + 1.0);
 93     }
 94   h_NA_infinity = alpha_hNA/(alpha_hNA+beta_hNA);
 95 
 96   double I_NA = g_NA_bar*(E_M-E_NA)
 97                         *m_NA_infinity*m_NA_infinity
 98       *h_NA_infinity; 
 99 
100   /*===============================================
101     Transient, Outward Potassium Current
102     ===============================================*/
103 
104   // activation parameter for I_K
105   double t_m_K1;
106   double m_K1_infinity;
107 
108   sum = E_M+42.0;
109   if(sum<0){
110     m_K1_infinity = exp(sum/13.0)/(exp(sum/13.0)+1.0);
111     }
112   else{
113     m_K1_infinity = 1.0/(1.0+exp(-sum/13.0));
114     }    
115 
116   // inactivation parameter for I_K
117   double h_K1_infinity;
118   double t_h_K1;
119 
120   sum = E_M+110.0;
121   if(sum<0){
122     h_K1_infinity = 1.0/(1.0+exp( sum/18.0));
123     }
124   else{
125     h_K1_infinity = exp( -sum/18.0)/(exp(-sum/18.0)+1.0);
126     }
127 
128   double I_K1 = g_K1_bar*(E_M-E_K)*m_K1_infinity*h_K1_infinity;
129   sum = I_NA + I_K1; 
130   //double g_leak_bar  = -sum/(E_M-E_L);
131   
132   cout << "Given g_leak_bar = 2.0 mu S find E_L" << endl;
133   double E_L = (2.0*E_M + sum)/2.0;
134   cout << "E_L = " << E_L << endl;
135   
136   
137   cout << "Let's check the equilibrium voltage:" << endl;
138   double g_NA_inf = g_NA_bar*m_NA_infinity*m_NA_infinity
139            *h_NA_infinity;
140   double g_K1_inf = g_K1_bar*m_K1_infinity*h_K1_infinity;          
141   double g_total = g_NA_inf+g_K1_inf+g_leak_bar;
142   double ratio_NA = g_NA_inf/g_total;
143   double ratio_K1 = g_K1_inf/g_total;
144   double ratio_L = g_leak_bar/g_total;
145   double E_Mcalc = ratio_NA*E_NA + ratio_K1*E_K + ratio_L*E_L;
146   
147   cout << "Calculated E_M = " << E_Mcalc << endl;
148   
149   /*------------------------------------------------
150     initial activations and inactivations:
151     m_NA_infinity
152     h_NA_infinity
153     m_K1_infinity
154     h_K1_infinity
155     ------------------------------------------------*/
156     
157     cout << "Initial activation and inactivations are" << endl;
158     cout << "m_NA(0) = " << m_NA_infinity << endl;
159     cout << "h_NA(0) = " << h_NA_infinity << endl;
160     cout << "m_K1(0) = " << m_K1_infinity << endl;
161     cout << "h_K1(0) = " << h_K1_infinity << endl;
162     
163     q[0] = m_NA_infinity;
164     q[1] = h_NA_infinity;
165     q[2] = m_K1_infinity;
166     q[3] = h_K1_infinity;
167 
168     E_L = -10.0;        
169     return E_L;
170 }
syntax highlighted by Code2HTML, v. 0.8.8b

2.4.2  C++ Code: The Run-Time Code:

The simulation code is still, of course, not as user friendly as we would like. As discussed in earlier simulation chapters, we control our simulation by editing the files run_plot.c and rkf5_plot.c to choose variables to plot and so forth. Then, we recompile the executible using standard unix make tools. Typically, we manage all of this with a unix Makefile: here is the one we use for this simulation:

MakeApp 
  1 # ------------------------- #
  2 # X/ Motif Paths            #
  3 # ------------------------- #
  4 MOTIFLIB = /usr/X11R6/lib
  5 XLIB = /usr/X11R6/lib
  6 # --------------------------------- #
  7 # Define X libraries location       #
  8 # ----------------------------------#
  9 XLIBS = -lXm  -lXp -lXpm -lXaw -lXt -lXmu -lXext -lX11 
 10 # ---------------------------- #
 11 # Define X include directories #
 12 # ---------------------------- #
 13 XINC  = /usr/X11R6/include
 14 MOTIFINC = /usr/X11R6/include
 15 XINCLUDES = -I$(XINC) -I$(MOTIFINC)
 16 #
 17 MYMAKE = MakeApp
 18 #
 19 # -------------------------- #
 20 # paths for this simulation  #
 21 # -------------------------- #
 22 BASE_PATH = /home/petersj/ProgramsApp
 23 SIM_PATH = $(BASE_PATH)/SimpleHH
 24 #
 25 MAIN_LIBRARY_PATH = /home/petersj/programs/JimLibs
 26 MAIN_INCLUDE_PATH = /home/petersj/programs/JimIncludes
 27 INCLUDE_PATH_APP = $(BASE_PATH)/sim_ode
 28 INCLUDES = $(XINCLUDES) -I$(MAIN_INCLUDE_PATH) \
 29   -I$(INCLUDE_PATH_XGRAPH) -I$(INCLUDE_PATH_APP)
 30 #
 31 LIBRARY_FLAGS_SOURCES = -L$(MAIN_LIBRARY_PATH)
 32 LIBRARY_PATH_X = $(XLIB)
 33 LIBRARY_PATH_MOTIF = $(MOTIFLIB)
 34 LIBRARY_FLAGS_SOURCES = -L$(MAIN_LIBRARY_PATH) \
 35   -L$(LIBRARY_PATH_MOTIF) -L$(LIBRARY_PATH_X)
 36 
 37 SUNLIBS = -ll -lgen
 38 LIBS  = -lm
 39 SHAPE_LIBS = -lbox -lshape
 40 VECTOR_LIBS = -lrkf5parobjects -ldoublevec2 -ldoublevec1 -lintvec2 -lintvec1  
 41 SIMLIBS =  $(VECTOR_LIBS) $(SHAPE_LIBS) 
 42 #
 43 # ----------- #
 44 # Definitions #
 45 # ----------- #
 46 DEFINES =
 47 CC = g++
 48 COMPILE_FLAGS_LINK = -g $(DEFINES) $(INCLUDES)
 49 COMPILE_FLAGS_SOURCES = -c -g $(DEFINES) $(INCLUDES)
 50 #
 51 RANLIB = touch
 52 # ------------- #
 53 # Define target #
 54 # ------------- #
 55 TARGET = ode
 56 OBJECTS = current.o \
 57   sim.o \
 58   popupdrawExposeCB.o \
 59   popupdrawResizeCB.o \
 60   setxcolor.o \
 61   find_named_color.o \
 62   update_color.o \
 63   force_update.o \
 64   popupquitCB.o \
 65   graph_color.o \
 66   axis_color.o \
 67   check_points.o \
 68   rkf5_plot.o \
 69   function.o \
 70   external_force.o \
 71   run_plot.o \
 72   set_plot.o \
 73   handle_button_release.o \
 74   imagepopupplot.o \
 75   set_popupplot.o \
 76   imagepopupcolor.o \
 77   set_axiscolor.o \
 78   set_graphcolor.o \
 79   set_active.o \
 80   expose.o \
 81   input.o \
 82   resize.o \
 83   quit_application.o \
 84   image.o \
 85   rest.o\
 86   gates.o
 87 SOURCES = $(OBJECTS:.o=.c)
 88 
 89 progs:
 90   @make $(TARGET) -f $(MYMAKE)
 91   
 92 $(TARGET): $(OBJECTS) simulation.h
 93   @echo
 94   @echo Creating the Executible $(TARGET)
 95   @echo
 96   rm -f $(TARGET)
 97   $(CC) $(COMPILE_FLAGS_LINK) -o $(TARGET) $(OBJECTS) \
 98   $(LIBRARY_FLAGS_SOURCES) $(SIMLIBS) $(XLIBS) $(LIBS)
 99   
100 # --------------------------- #
101 # 'make' template for OBJECTS #
102 # --------------------------- # 
103 $(OBJECTS): $(@:.o=.c) simulation.h
104   @echo
105   @echo Compiling $@
106   @echo
107   $(CC) $(COMPILE_FLAGS_SOURCES) $(LIBS) $(@:.o=.c)
108 all:
109   @make progs -f $(MYMAKE)
110 #--------------------------------------------------------
111 #      dependencies
112 #--------------------------------------------------------
113 current.o: current.c simulation.h
114 sim.o: sim.c simulation.h 
115 graph_color.o: graph_color.c simulation.h
116 axis_color.o: axis_color.c simulation.h
117 check_points.o: check_points.c simulation.h 
118 rkf5_plot.o: rkf5_plot.c simulation.h
119 function.o: function.c simulation.h 
120 external_force.o: external_force.c simulation.h 
121 run_plot.o: run_plot.c simulation.h
122 handle_button_release.o: handle_button_release.c simulation.h  
123 set_axiscolor.o: set_axiscolor.c simulation.h
124 set_graphcolor.o: set_graphcolor.c simulation.h
125 set_popupplot.o: set_popupplot.c simulation.h
126 imagepopupcolor.o: imagepopupcolor.c simulation.h 
127 imagepopupplot.o: imagepopupplot.c simulation.h
128 set_active.o: set_active.c simulation.h 
129 set_plot.o: set_plot.c simulation.h   
130 expose.o: expose.c simulation.h 
131 input.o: input.c simulation.h 
132 resize.o: resize.c simulation.h 
133 quit_application.o: quit_application.c simulation.h 
134 image.o: image.c simulation.h 
135 scale.o: scale.c simulation.h 
136 rotate.o: rotate.c simulation.h
137 popupdrawExposeCB.o: popupdrawExposeCB.c simulation.h
138 popupdrawResizeCB.o: popupdrawResizeCB.c simulation.h
139 popupquitCB.o: popupquitCB.c simulation.h
140 setxcolor.o: setxcolor.c simulation.h
141 find_named_color.o: find_named_color.c simulation.h
142 force_update.o: force_update.c simulation.h
143 update_color.o: update_color.c simulation.h
144 rest.o: rest.c simulation.h
145 gates.o: gates.c simulation.h
syntax highlighted by Code2HTML, v. 0.8.8b

Hence after editing these two files as you wish, you recompile as follows:

run1.prompts 
 1 [petersj@dagmar SimpleHH]$ make -f MakeApp
 2 make[1]: Entering directory `/home/petersj/ProgramApps/CellModels/SimpleHH'
 3 
 4 Compiling run_plot.o
 5 
 6 g++ -c -g  -I/usr/X11R6/include -I/usr/X11R6/include -I/home/petersj/programs/JimIncludes -I -I/home/petersj/ProgramsApp/sim_ode -lm run_plot.c
 7 g++: -lm: linker input file unused since linking not done
 8 
 9 Creating the Executible ode
10 
11 rm -f ode
12 g++ -g  -I/usr/X11R6/include -I/usr/X11R6/include -I/home/petersj/programs/JimIncludes -I -I/home/petersj/ProgramsApp/sim_ode -o ode current.o sim.o popupdrawExposeCB.o popupdrawResizeCB.o setxcolor.o find_named_color.o update_color.o force_update.o popupquitCB.o graph_color.o axis_color.o check_points.o rkf5_plot.o function.o external_force.o run_plot.o set_plot.o handle_button_release.o imagepopupplot.o set_popupplot.o imagepopupcolor.o set_axiscolor.o set_graphcolor.o set_active.o expose.o input.o resize.o quit_application.o image.o rest.o gates.o \
13 -L/home/petersj/programs/JimLibs -L/usr/X11R6/lib -L/usr/X11R6/lib -lrkf5parobjects -ldoublevec2 -ldoublevec1 -lintvec2 -lintvec1  -lbox -lshape   -lXm  -lXp -lXpm -lXaw -lXt -lXmu -lXext -lX11  -lm
14 make[1]: Leaving directory `/home/petersj/ProgramApps/CellModels/SimpleHH'
15 
16
syntax highlighted by Code2HTML, v. 0.8.8b Then to run the simulation, you type ode at the prompt. You will then see the popup windows we have discussed and see some terminal output in the window you start the simulation in. The output shown here discards some of the X Window diagnostic code that is still present in the simulation for our debugging purposes!

run2.prompts 
 1 [petersj@dagmar SimpleHH]$ ode
 2 
 3 Let's plot the ODE data:
 4 Let's plot action potential:
 5 In rest calculation: parameters are p = 2 2.77075 0.02 57.11 -71.9989 0 
 6 
 7 Given g_leak_bar = 2.0 mu S find E_L
 8 E_L = -59.8051
 9 Let's check the equilibrium voltage:
10 Calculated E_M = -67.352
11 Initial activation and inactivations are
12 m_NA(0) = 9.88698e-05
13 h_NA(0) = 0.987574
14 m_K1(0) = 0.200269
15 h_K1(0) = 0.0585369
16 E_NA        =  57.11
17 E_K         =  -71.9989
18 E_M         =  -60
19 E_L         =  -10
20 g_leak_bar  =  0.02
21 yinit = -60 9.88698e-05 0.987574 0.200269 0.0585369 
22 
23 Entering rkf5_plot code
24 parameters = 2 2.77075 0.02 57.11 -71.9989 -10 
25 
26 parameters = 2 2.77075 0.02 57.11 -71.9989 -10 
27 
28 (t,y[0]) = 0.05,-59.7984
29 (t,y[0]) = 0.1,-59.6003
30 (t,y[0]) = 0.15,-59.4057
31 (t,y[0]) = 0.2,-59.2148
32 (t,y[0]) = 0.25,-59.0273
33 (t,y[0]) = 10.05,-46.8455
34 (t,y[0]) = 20.05,-53.4617
35 (t,y[0]) = 30.05,-53.0657
36 (t,y[0]) = 40.05,-52.9331
37 (t,y[0]) = 50.05,-52.8046
38 (t,y[0]) = 60.05,-52.6794
39 (t,y[0]) = 70.05,-52.5575
40 Finished with PLOT loop 
41 Returned from oderkf5_plot call
syntax highlighted by Code2HTML, v. 0.8.8b

run_plot.c:

The initialization of the simulation itself is handled by the callback run_plot.c which sets the parameters for the ODE solver, the initial conditions for the ion concentrations and so forth.

run_plot.c 
  1 #include "simulation.h"
  2 
  3 void run_plot(Widget w,XtPointer client_data,XtPointer call_data)
  4 {
  5   XmPushButtonCallbackStruct *P = (XmPushButtonCallbackStruct *)call_data;
  6   application_data *T = (application_data *)client_data;
  7 
  8   printf("Let's plot action potential:\n");
  9   int i;
 10   
 11   int size = 5;
 12   double tol,tinit,tend,hinit;
 13 
 14   /*======================================
 15     Allocate memory for state variables
 16     ======================================*/
 17     DOUBLE_VECTOR yinit(size);
 18     DOUBLE_VECTOR y(size);
 19     double C = 14.28;
 20         
 21   /*======================================
 22     Allocate memory for dynamic parameters
 23     p[0] g_NA_bar
 24     p[1] g_K1_bar
 25     p[2] g_leak_bar
 26     p[3] E_NA
 27     p[4] E_K
 28     p[5] E_L
 29     
 30     q[0] = m_NA(0);
 31     q[1] = h_NA(0);
 32     q[2] = m_K1(0);
 33     q[3] = h_K1(0);    
 34     ======================================*/  
 35       
 36     DOUBLE_VECTOR p(6);
 37     DOUBLE_VECTOR q(4);
 38     
 39 /*----------------------------------------
 40     Koch values  
 41   ----------------------------------------*/   
 42     double g_NA_bar     =   2.0;
 43     double g_K1_bar     =  0.120+1.17+0.084+1.2+0.054
 44                           +0.02675+0.116; //2.7708
 45     double g_leak_bar   = .02;
 46     
 47     double NA_O       = 491.0;
 48     double NA_I       =  50.0;
 49     double K_O        =   7.859;
 50     double K_I        = 140.0;    
 51             
 52     double R = 8.31;       //Rydberg's Constant
 53     double DegK = 276.0 + C; //Kelvin Temperature; use C degrees Celsius
 54     double F = 9.649e+4;     //Faraday's constant
 55     double RTF = R*DegK/F*1e+3; 
 56         
 57     double E_NA;
 58     E_NA = RTF*log(NA_O/NA_I);
 59     double E_K;
 60     E_K = RTF*log(K_O/K_I);    
 61     
 62     double E_M = -60.0;
 63           
 64     p[0]  = g_NA_bar;
 65     p[1]  = g_K1_bar;
 66     p[2]  = g_leak_bar;
 67     p[3]  = E_NA;
 68     p[4]  = E_K;
 69     p[5]  = 0.0;
 70     
 71     double E_L = rest(E_M,p,q);
 72     p[5]  = E_L;
 73     
 74     cout << "E_NA        =  " << E_NA << endl;
 75     cout << "E_K         =  " << E_K << endl;     
 76     cout << "E_M         =  " << E_M << endl;
 77     cout << "E_L         =  " << E_L << endl;   
 78     cout << "g_leak_bar  =  " << g_leak_bar << endl;        
 79     
 80   /*======================================
 81     Set up initial conditions
 82     ======================================*/
 83     T->initial = new DOUBLE_VECTOR(size);     
 84     DOUBLE_VECTOR& INIT = *(T->initial);
 85     
 86     INIT[0] =   E_M;
 87     INIT[1] =   q[0];
 88     INIT[2] =   q[1];
 89     INIT[3] =   q[2];
 90     INIT[4] =   q[3];
 91 
 92     yinit = INIT;
 93     y = INIT;
 94     cout << "yinit = " << yinit << endl;
 95 
 96     // y vector assignments
 97   /*===============================================
 98     y[0] = V
 99     y[1] = m_NA
100     y[2] = h_NA
101     y[3] = m_K1
102     y[4] = h_K1
103     
104     plot_type = 0 ==> y[0]  range [-80,100]
105     plot_type = 1 ==> y[1]  range [0,1] 
106     plot_type = 2 ==> y[2]  range [0,1]
107     plot_type = 3 ==> y[3]  range [0,1]
108     plot_type = 4 ==> y[4]  range [0,6.0e-2] 
109     
110     plot_type = 5  ==> alpha_m_NA     range [0,50]
111     plot_type = 6  ==> beta_m_NA      range [0,50]
112     plot_type = 7  ==> alpha_h_NA     range [0,1]
113     plot_type = 8  ==> beta_h_NA      range [0,8] 
114     plot_type = 9  ==> m_NA_infinity  range [0,1]    
115     plot_type = 10 ==> h_NA_infinity  range [0,1]
116     plot_type = 11 ==> m_K1_infinity  range [0,1]
117     plot_type = 12 ==> h_K1_infinity  range [0,0.08]
118     plot_type = 13  ==>t_m_NA         range [0,1]
119     plot_type = 14 ==> t_h_NA         range [0,1]    
120     plot_type = 15 ==> t_m_K1         range [0,1]    
121     plot_type = 16 ==> t_h_K1         range [0,1]                                  
122     ===============================================*/
123     DOUBLE_VECTOR Low_Range(17);
124     DOUBLE_VECTOR High_Range(17);
125     Low_Range[0]    = -80.0;
126     High_Range[0]   = 100.0;    
127     for(i=1;i<5;++i){
128       Low_Range[i]  = 0.0;
129       High_Range[i] = 1.0;
130       }
131     High_Range[4] = 6.0e-2; 
132     for(i=5;i<7;++i){
133       Low_Range[i]  =  0.0;
134       High_Range[i] = 50.0;
135       } 
136     for(i=7;i<17;++i){
137       Low_Range[i]  = 0.0;
138       High_Range[i] = 1.0;
139       } 
140     High_Range[8] = 8.0;
141     High_Range[12] = 0.08;  
142     High_Range[14] = 1.0e+1;          
143     int plot_type;                   
144   /*========================================
145     Set up parameters for ode integration
146     ========================================*/
147     hinit          =   1.0e-6;
148     tinit          =   0.00;
149     tend           =  80.0;
150     tol            =   1.0e-14;
151     plot_type  = 0;
152     double ymin;
153     double ymax;
154     for(i=0;i<17;++i){
155       if(plot_type==i){
156         ymin = Low_Range[i];
157         ymax = High_Range[i];
158         }
159       }
160     T->PLOT->plot_type = plot_type;
161     T->PLOT->x_minimum = tinit;
162     T->PLOT->x_maximum = tend;
163     T->PLOT->y_minimum = ymin;
164     T->PLOT->y_maximum = ymax;    
165     
166     cout << "Entering rkf5_plot code" << endl;
167     cout << "parameters = " << p << endl;
168     oderkf5_plot(yinit,y,p,tol,tinit,tend,hinit,T,w);
169     cout << "Returned from oderkf5_plot call " << endl;
170     
171 }
syntax highlighted by Code2HTML, v. 0.8.8b

rkf5_plot.c:

The actual call to the ODE solver is handled by the utility function rkf5_plot.c which as described in Chapter 6 splits the integration interval into manageable chunks for plotting purposes.

rkf5_plot.c 
  1 #include "simulation.h"
  2 #include <iostream.h>
  3 /*==========================================================*/
  4 /*                                                          */
  5 /* Runga-Fehlberg order 5 method for plotting               */
  6 /*                                                          */
  7 /* yinit          = initial state                           */
  8 /* y              = calculated state                        */
  9 /* u              = control parameter                      */
 10 /* tinit          = initial t                               */
 11 /* tend           = final t                                 */
 12 /* tol            = tolerance for step size changing        */
 13 /* hinit          = initial step size                       */
 14 /*==========================================================*/
 15 void oderkf5_plot(DOUBLE_VECTOR& yinit,
 16                   DOUBLE_VECTOR& y,DOUBLE_VECTOR& p,double tol,
 17             double tinit,double tend,double hinit,
 18             application_data *T,
 19                   Widget w)
 20 {
 21 //open a file for plot data
 22 ofstream fd2;
 23 fd2.open("mK1inf", ios::in);
 24 DOUBLE_VECTOR GATES(14);
 25   /*========================================
 26   Gates[0]  = alpha_m_NA;
 27   Gates[1]  = beta_m_NA;
 28   Gates[2]  = alpha_h_NA;
 29   Gates[3]  = beta_h_NA;
 30   Gates[4]  = m_NA_infinity;
 31   Gates[5]  = h_NA_infinity;
 32   Gates[6]  = m_K1_infinity;
 33   Gates[7]  = h_K1_infinity;
 34   Gates[8]  = t_m_NA;
 35   Gates[9]  = t_h_NA; 
 36   Gates[10] = t_m_K1;   
 37   Gates[11] = t_h_K1; 
 38   Gates[12] = I_NA;
 39   Gates[13] = I_K1;
 40   ============================================*/
 41 int POINTS_PER_UNIT = 20;
 42 int NUMBER_PLOT_POINTS = int(tend)*POINTS_PER_UNIT;
 43 float HPLOT = (tend-tinit)/(float)NUMBER_PLOT_POINTS;
 44 double tfinal;
 45 int i,j,k,plot_diagnostic;
 46 plot_data *PLOT = (plot_data *)(T->POP->Outside_Data);
 47 int plot_type = PLOT->plot_type;
 48 cout << "parameters = " << p << endl;
 49 int FREQ = 200;
 50 
 51   //set diagnostic status level
 52   plot_diagnostic = 0;
 53 
 54   PLOT->draw_axis = 1;
 55   tfinal = tinit + HPLOT;
 56   DOUBLE_VECTOR y0 = yinit;
 57   gates(y[0],GATES,p);
 58   PLOT->first_time = tinit;
 59   if(0<=plot_type && plot_type<=4){
 60     PLOT->first_ordinate = y[plot_type];
 61     }
 62   else if(5<=plot_type && plot_type<17){
 63     PLOT->first_ordinate = GATES[plot_type-5];
 64     }
 65   fd2 << PLOT->first_time << " " << PLOT->first_ordinate << endl;  
 66   for(int count=0; count<NUMBER_PLOT_POINTS; ++ count){
 67     int debug = 0;
 68     int display_counter = 1;
 69     int plot_counter = 1;
 70  
 71     //integrate from tinit to tfinal 
 72     oderkf5_parobjects(y0,y,p,tol,tinit,tfinal,hinit);
 73     gates(y[0],GATES,p);
 74     if( count%FREQ==0 || count < 5){
 75       if(0<=plot_type && plot_type<=4){
 76         cout << "(t,y[" << plot_type << "]) = " << tfinal  
 77              << "," << y[plot_type] << endl;
 78   }
 79       else if(5<=plot_type && plot_type<17){
 80         cout << "(t,GATES[" << plot_type-5 << "]) = " << tfinal  
 81              << "," << GATES[plot_type-5] << endl;      
 82         }
 83       }
 84 
 85     PLOT->second_time = tfinal;
 86     if(0<=plot_type && plot_type<=4){
 87       PLOT->second_ordinate = y[plot_type];
 88       }
 89     else if(5<=plot_type && plot_type<17){
 90       PLOT->second_ordinate = GATES[plot_type-5];
 91       }        
 92     if(plot_diagnostic==1){
 93       printf("first_time = %8.2f first_ordinate = %8.2f\n",
 94       PLOT->first_time,PLOT->first_ordinate);
 95       printf("second_time = %8.2f second_ordinate = %8.2f\n",
 96       PLOT->second_time,PLOT->second_ordinate);
 97       }
 98     fd2 << PLOT->second_time << " " << PLOT->second_ordinate << endl;
 99     T->POP->Image(w,T->POP);
100     PLOT->draw_axis = 0;
101     PLOT->first_time = PLOT->second_time;
102     PLOT->first_ordinate = PLOT->second_ordinate;  
103     //reset yinit
104     y0 = y;
105     tinit = tfinal;
106     tfinal += HPLOT;
107     }//plot points loop
108 fd2.close();    
109 cout << "Finished with PLOT loop " << endl;
110 }
syntax highlighted by Code2HTML, v. 0.8.8b

Run-time Results:

We are now ready to begin our numerical simulations. We will use standard numerical solution techniques for the equation
CM
dV
dt
= - Im + IK1 + INa + IL
tNAm
dmNA
dt
= mNA¥ - mNA
tNAh
dhNA
dt
= hNA¥ - hNA
tK1m
dmK1
dt
= mK1¥ - mK1
tK1h
dhK1
dt
= hK1¥ - hK1

with initial conditions
V(0) = -60.00
mNA(-60) = 9.88698e-05
hNA(-60) = 0.987574
mK1(-60) = 0.200269
hK1(-60) = 0.0585369

For our numerical work, we let our variable vector be y and use the assignments:
y[0] = V
y[1] = mNA
y[2] = hNA
y[3] = mK1
y[4] = hK1

If we run the simulation with the parameters as set in run_plot.c and rest.c we see the voltage traces of Figure 2.8 for three different starting membrane voltages.

Figure 2.8: Simple Hodgkin - Huxley Dynamics:

2.4.3  MatLab Simulation Code:

This MatLab is more complicated and it is worth your while to have studied the previous chapter on how to solve Ordinary Differential Equations in MatLab. Simple code to manage our solution of the Hodgkin-Huxley equations is given below. We simply integrate the equations using a fixed time step repeatedly with a Runga-Kutte method of order 4. This is not terribly accurate but it illustrates the results nicely nevertheless.

SolveSimpleHH.m

The code to manage the solution of the ODE is given below:

SolveSimpleHH.m 
 1 function [tvals,g_NA,g_K,V] = SolveSimpleHH(fname,iename,t0,y0,h,k,n,...
 2                           g_Kmax,g_NAmax)
 3 
 4 %
 5 % We use a simple Runge-Kutta scheme
 6 % using the Matlab function below:
 7 %
 8 %function [tvals,g_Na,g_K,V] = HHFixedRK(fname,iename,t0,y0,h,k,n)
 9 %
10 %        Gives approximate solution to 
11 %          y'(t) = f(t,y(t))
12 %          y(t0) = y0
13 %        using a kth order RK method
14 %
15 % t0     initial time
16 % y0     initial state
17 % h      stepsize
18 % k      RK order  1<= k <= 4
19 % n      Number of steps to take
20 %
21 % tvals  time values of form
22 %        tvals(j) = t0 + (j-1)*h, 1 <= j <= n
23 % 
24 % g_Na = sodium conductance  
25 % g_K  = potassium conductance
26 % V    = membrane voltage 
27 %
28 m_NA = zeros(1,n);
29 h_NA = zeros(1,n);
30 n_K  = zeros(1,n);
31 g_Na = zeros(1,n);
32 g_K  = zeros(1,n);
33 V    = zeros(1,n);
34 
35 [tvals,yvals] = HHFixedRK(fname,iename,t0,y0,h,k,n);
36 
37 %
38 % store values in physical variables
39 %
40 V    = yvals(1,1:n);
41 m_NA = yvals(2,1:n);
42 h_NA = yvals(3,1:n);
43 n_K  = yvals(4,1:n);
44 
45 %
46 %  Compute g_Na and G_K
47 %
48 for i = 1:n
49   u       = m_NA(i)*m_NA(i)*m_NA(i)*h_NA(i);
50   g_NA(i) = g_NAmax*u;
51 end
52 
53 for i = 1:n
54   u       = n_K(i)*n_K(i)*n_K(i)*n_K(i);
55   g_K(i)  = g_Kmax * u;
56 end
syntax highlighted by Code2HTML, v. 0.8.8b

simpleHH.m

The Hodgkin-Huxley dynamics are encoded in the function below: simpleHH.m 
  1 function f = simpleHH(t,y,IE)
  2                     
  3 %    Standard Hodgkin - Huxley Model
  4 %    voltage       mV
  5 %    current       na
  6 %    time          ms
  7 %    concentration mM
  8 %    conductance   micro Siemens
  9 %    capacitance   nF
 10 %    ===============================================
 11 %    y vector assignments
 12 %    ===============================================
 13 %    y(1) = V
 14 %    y(2) = m_NA
 15 %    y(3) = h_NA
 16 %    y(4) = m_K
 17 %
 18 % set size of f
 19 %
 20 
 21 f = zeros(1,4);
 22 
 23 %    ===============================================
 24 %    f vector assignments
 25 %    ===============================================
 26 %    y(1) = V dynamics
 27 %    y(2) = m_NA dynamics
 28 %    y(3) = h_NA dynamics
 29 %    y(4) = m_K dynamics
 30 %
 31 %   ================================================
 32 %    Fast Sodium Current
 33 %    ===============================================
 34 %   ================================================
 35 %    Constants for Equilibrium Voltage Calculations
 36 %    ===============================================
 37 %
 38 
 39 % Rydberg's Constant
 40   R = 8.31;
 41 % Kelvin Temperature; use 9.3 degrees Celsius
 42   T = 9.3;
 43 % Faraday's constant
 44   F = 9.649e+4;
 45 %
 46 % Compute Nernst voltages of E_K and E_Na
 47 % voltage = Nernst(valence,Temperature,InConc,OutConc)
 48 %
 49 
 50 % 
 51 % Sodium
 52 %
 53   NA_O = 491.0;
 54   NA_I = 50.0; 
 55   E_NA = Nernst(1,T,NA_I,NA_O);
 56    
 57 %
 58 % Potassium
 59 %
 60   K_O = 20.11;
 61   K_I = 400.0;
 62   E_K = Nernst(1,T,K_I,K_O);
 63   
 64 
 65 % max conductance for NA
 66   g_NA_bar = 120.0;
 67   
 68 % max conductance for K  
 69   g_K_bar = 3.60;
 70 
 71 %
 72 % activation/inactivation parameters for NA
 73 %
 74 % alpha_mNA, beta_mNA
 75 %
 76   sum = y(1)+35.0;
 77   if sum > 0 
 78     alpha_mNA     = -0.10*sum/(exp(-sum/10.0) - 1.0);
 79   else
 80     alpha_mNA     = -0.10*exp( sum/ 10.0)*sum/(1.0 - exp( sum/ 10.0));
 81   end
 82 %  
 83   sum = y(1)+60.0;
 84   if sum > 0
 85     beta_mNA      =  4.0*exp( -sum/18.0);
 86   else
 87     beta_mNA      =  4.0*exp( -sum/18.0);
 88   end
 89 %  
 90   m_NA_infinity = alpha_mNA/(alpha_mNA+beta_mNA);
 91   t_m_NA        = 1.0/(alpha_mNA+beta_mNA);
 92 %
 93   f(2)          = (m_NA_infinity - y(2))/t_m_NA;
 94 
 95 %
 96 % activation, inactivation parameter for I_NA
 97 %
 98 % alpha_hNA, beta_mNA
 99 %
100   sum = (y(1)+60.0)/20.0;
101   if sum < 0
102     alpha_hNA     =  0.07*exp( -sum);
103   else
104     alpha_hNA     =  0.07*exp( -sum);
105   end
106 %
107   sum = y(1)+30.0;
108   if sum>0 
109     beta_hNA      =  1.0/(1.0 + exp(-sum/10.0));
110   else
111     beta_hNA      =  exp( sum/10.0)/(exp( sum/10.0) + 1.0);
112   end
113 %
114   h_NA_infinity = alpha_hNA/(alpha_hNA+beta_hNA);
115   t_h_NA        = 1.0/(alpha_hNA+beta_hNA);
116   f(3)          = (h_NA_infinity - y(3))/t_h_NA;
117 
118 %
119 % I_NA current
120 %
121   I_NA = g_NA_bar*(y(1)-E_NA)*y(2)*y(2)*y(2)*y(3); 
122 
123 %
124 % activation/inactivation parameters for K
125 %
126 % alpha_mK
127 %
128   sum = y(1)+50.0;
129   if sum > 0
130     alpha_mK     = -0.01*sum/(exp(-sum/10.0) - 1.0);
131   else
132     alpha_mK     = -0.01*exp(sum/10.0)*sum/(1.0 - exp(sum/10.0));
133   end  
134   sum = (y(1)+60.0)*0.0125;
135   if sum > 0
136     beta_mK      = 0.125*exp( -sum);
137   else
138     beta_mK      = 0.125*exp( -sum);
139   end
140   t_m_K          = 1.0/(alpha_mK+beta_mK);
141   m_K_infinity   = alpha_mK/(alpha_mK+beta_mK);
142   f(4)           = (m_K_infinity - y(4))/t_m_K;
143 
144 %
145 % I_K current
146 %
147   I_K = g_K_bar*(y(1)-E_K)*y(4)*y(4)*y(4)*y(4);
148 
149 % 
150 % leakage current
151 %
152   g_leak = 0.30;
153   E_leak = -49.0;
154 
155 %
156 % I_L current
157 %   
158   I_leak = g_leak*(y(1)-E_leak);
159 
160 %
161 % Cell Capacitance
162 %
163   C_M = 1.0;
164   s = feval(IE,t);
165   f(1) = (s - I_NA - I_K - I_leak)/C_M;
syntax highlighted by Code2HTML, v. 0.8.8b

HHFixedRK.m

The code that manages the call to a fixed Runga-Kutte order method from an earlier chapter has been modified to accept an argument which is the name of the function which models the current injection.

HHFixedRK.m 
 1 function [tvals,yvals] = HHFixedRK(fname,iename,t0,y0,h,k,n)
 2 %
 3 %        Gives approximate solution to 
 4 %          y'(t) = f(t,y(t))
 5 %          y(t0) = y0
 6 %        using a kth order RK method
 7 %
 8 % t0     initial time
 9 % y0     initial state
10 % h      stepsize
11 % k      RK order  1<= k <= 4
12 % n      Number of steps to take
13 %
14 % tvals  time values of form
15 %        tvals(j) = t0 + (j-1)*h, 1 <= j <= n
16 % yvals  approximate solution
17 %        yvals(:j) = approximate solution at 
18 %        tvals(j),  1 <= j <= n
19 % fname  name of dynamics function
20 % iename name of current injection function
21 %
22 tc = t0;
23 yc = y0;
24 tvals = tc;
25 yvals = zeros(4,n);
26 yvals(1:4,1) = transpose(yc);
27 fc = feval(fname,tc,yc,iename);
28 for j=1:n-1
29   [tc,yc,fc] = HHRKstep(fname,iename,tc,yc,fc,h,k);
30   yvals(1:4,j+1) = transpose(yc);
31   tvals = [tvals tc];
32 end
syntax highlighted by Code2HTML, v. 0.8.8b

HHRKstep.m

The actual Runga-Kutte code was also modified to allow the name of the current injection function to be entered as an argument.

HHRKstep.m 
 1 function [tnew,ynew,fnew] = HHRKstep(fname,iename,tc,yc,fc,h,k)
 2 %
 3 % fname    the name of the right hand side function f(t,y)
 4 %          t is a scalar usually called time and
 5 %          y is a vector of size d
 6 % yc       approximate solution to y'(t) = f(t,y(t)) at t=tc
 7 % fc       f(tc,yc)
 8 % h        The time step
 9 % k        The order of the Runge-Kutta Method 1<= k <= 4
10 %
11 % tnew     tc+h
12 % ynew     approximate solution at tnew
13 % fnew     f(tnew,ynew)
14 %
15 if k==1
16   k1 = h*fc;
17   ynew = yc+k1;
18 elseif k==2
19   k1 = h*fc;
20   k2 = h*feval(fname,tc+(h/2),yc+(k1/2),iename);
21   ynew = yc + (k1+k2)/2;
22 elseif k==3
23   k1 = h*fc;
24   k2 = h*feval(fname,tc+(h/2),yc+(k1/2),iename);
25   k3 = h*feval(fname,tc+h,yc-k1+2*k2,iename);
26   ynew = yc+(k1+4*k2+k3)/6;
27 elseif k==4
28   k1 = h*fc;
29   k2 = h*feval(fname,tc+(h/2),yc+(k1/2),iename);
30   k3 = h*feval(fname,tc+(h/2),yc+(k2/2),iename);
31   k4 = h*feval(fname,tc+h,yc+k3,iename);
32   ynew = yc+(k1+2*k2+2*k3+k4)/6;
33 else
34   disp(sprintf('The RK method %2d order is not allowed!',k));    
35 end
36 tnew = tc+h;
37 fnew = feval(fname,tnew,ynew,iename);
syntax highlighted by Code2HTML, v. 0.8.8b

RunTime Results:

We ran our simple simulation using this session:

runtime1.prompts 
1 >> [tvals,g_NA,g_K,V] = SolveSimpleHH('simpleHH','IE',0,y,1.0e-3,4,...
2                                         10000,3.6,120);
3 >> plot(tvals,V);
4 >> plot(tvals,g_NA);
5 >> plot(tvals,g_K);
6 >> plot(tvals,g_NA,'g-',tvals,g_K,'g-.');
7 >>
syntax highlighted by Code2HTML, v. 0.8.8b In Figure ( 2.9), we se the membrane voltage versus time curve. This is not as smooth as the one we generated in the X windows simulation because we are using a fixed Runga-Kutte method wihtout adaptive step size control.

Figure 2.9: Membrane Voltage vs Time:

In Figure ( 2.10), we se the ion conductance versus time curves. Note that gK is substantially smaller than gNA.

Figure 2.10: Ion Conductances vs Time:

Previous Contents Next