Turk J Elec Eng & Comp Sci, Vol.18, No.4, 2010, c T¨UB˙ITAK
doi:10.3906/elk-0904-20
Washout filter based control for the Hodgkin-Huxley
nerve cell dynamics
Re¸sat ¨Ozg¨ur DORUK
Middle East Technical University-Northern Cyprus Kalkanlı,
G¨uzelyurt, KKTC (TRNC)
e-mail: rdoruk@metu.edu.tr
Abstract
In this work we present a washout filter based control approach for stabilizing the oscillations in HodgkinHuxley
type neurons. The oscillations occur due to the bifurcations that arise from the potassium and sodium
channel conductance deviations. The MATCONT toolbox software environment was used for analysis of the
bifurcation points in conjunction with MATLAB. The simulations of the Hodgkin-Huxley (HH) model at
those points are provided to validate the results obtained from the MATCONT software. Then a washout
filter is proposed to stabilize the oscillations occur at the Hopf bifurcation points where the external current
injection is used to provide the control actuation. A second simulation set is provided to demonstrate the
action of the washout filter.
Key Words: Hopf bifurcation, saddle bifurcation, Hodgkin-Huxley equation, potassium channel, sodium
channel
1. Introduction
The Hodgkin-Huxley nonlinear model [1]of a real neuron (giant squid axon at first) is one of the biggest challenges
in the electrophysiology field in the near history. The cell membrane potential is modeled using the basic circuit
theory and some chemical kinetics. As a result a fourth order highly nonlinear differential equation is obtained.
This model possesses a rich set of bifurcations due to the additional current injection as described by [2] and
external electric fields [3,4]. Bifurcations are problematic in biological neural networks since the diseases such
as Parkinson’s [5], epilepsy [6] and schizophrenia [7, 8] have some connections with the instabilities in the neural
cell dynamics. In addition, various cardiac and muscle diseases possess chaotic behaviors [9] during the course
of a particular disease. Because of that, controlling bifurcations in the neuron dynamics may bring some new
insights in the treatment of some neurological disorders. Local bifurcation analysis is an important branch of
the nonlinear system theory where the properties of the equilibrium points of a nonlinear system changes due
to a single parameter deviation [10]. One of the mostly met bifurcation type is the Hopf type bifurcation where
553
Turk J Elec Eng & Comp Sci, Vol.18, No.4, 2010
the equilibrium point loses its stability as a result of a single pair of purely complex eigenvalues arising from
a parameter drift. There are also other critical points such as the saddle point where there is a real pair of
eigenvalues one of which is on the left hand side and the other is the right hand side of the complex plane (with
equal values). In order to solve the stability issue in such a case, several studies have been proposed. Some
of those methods are the linear delayed feedback [11, 12], nonlinear feedback [13] and washout filters [14, 15,
16]. The advantage of the washout filters is that they do not change the equilibrium points of the uncontrolled
(open loop) nonlinear system. The washout filter is originally a high-pass filter which does not pass steady state
inputs (only allows the transients). An application of the washout filters for the Hodgkin-Huxley type nerve cell
dynamics are presented in [2] for the bifurcations existent in the neuron model where the variable parameter
is the possible external current injection (from the cell environment) into the nerve cell. That attempts to
control the bifurcations by applying an external current stimulus. In this work we present another bifurcation
control mechanism based on the extended washout filter for the bifurcations resulting from the deviations in
the sodium and potassium channel conductances. In the first attempt, the MATCONT toolbox [17] is used to
detect particular Hopf bifurcation points in the Hodgkin-Huxley nerve cell dynamics. Secondly, the points are
simulated in MATLAB in order to verify the oscillations occurred at the detected parameter sets. The approach
of extended washout filter is applied to the nerve cell dynamics by driving the external current input and the
MATLAB simulations are repeated in order to check whether the oscillations decease or not. In addition to the
simulations, a bifurcation analysis of the closed loop against the washout filter coefficients by the MATCONT
package is performed in order to investigate the stability of the HH-Washout Filter combination.
2. Hodgkin-Huxley model
The Hodgkin-Huxley model is a fourth order nonlinear differential equation derived for representing the action
potential dynamics of the giant squid axon. The presented model is explained in detail in the resource [1]. In
this section only the resultant equations are given with the nominal parameter values:
Cm
dV
dt = −gNam3
h (V − VNa) − gKn4
(V − VK) − gL (V − VL) + Iext
dn
dt
= αn (1 − n) − βnn
dm
dt = αm (1 − m) − βmm
dh
dt = αh (1 − h) − βhh
αn = 0.1−0.01V
exp(1−0.1V )−1
, βn = 0.125 exp − V
80
αm = 2.5−0.1V
exp(2.5−0.1V )−1 , βm = 4 exp − V
18
αh = 0.07 exp − V
20 , βm = 1
exp(3−0.1V )+1
(1)
Definitions of the parameters are:
V is the displacement of the cell membrane potential from its resting value in mV;
554
DORUK: Washout filter based control for the Hodgkin-Huxley nerve cell dynamics,
n is a dimensional variable which can vary between 0 and 1 representing the proportion of activating
molecules of the potassium channel;
m is same type of variable as n which represents the proportion of the activating molecules of the sodium
channel;
h same type of variable as n which represents the proportion of the inactivating molecules of the sodium
channel;
Iext denotes the external current injection, in units of μA
cm2 ;
gNa denotes the sodium channel conductance, in units of mS
cm2 ;
gk denotes the potassium channel conductance, in units of mS
cm2 ;
gL denotes the current leakage conductance, in units of mS
cm2 (due to the chloride and other ions);
Cm denotes the membrane capacitance, in units of μF
cm2 ;
VNa denotes the sodium channel resting potential, in units of mV;
VK denotes the potassium channel resting potential, in units of mV;
VL denotes the the level of potential where the leakage current reduces to zero.
The following values were used in the simulations:
gK = 36mS/cm2
gNa = 120mS/cm2
gL = 0.3mS/cm2
VK = −12 mV
VNa = 115 mV
VL = 10.613 mV
Cm = 0.91 μF/cm2
. (2)
In this work the external current injection is set to zero by default. It is to be used as an external stimulus
for stabilizing the neuron with the washout filter. In the next section the Hopf bifurcation analysis results
produced by the MATCONT toolbox are presented for the variations occurred in the sodium and potassium
channel conductance parameters.
3. Bifurcation analysis results
The bifurcation analysis is done by first finding an equilibrium point of the dynamical system given in (1) by
simulating with zero external current(Iext = 0). According to the simulation results the Hodgkin-Huxley model
with the nominal parameters in (2) the following equilibrium point is found:
555
Turk J Elec Eng & Comp Sci, Vol.18, No.4, 2010
V = 0.003620359468
n = 0.31773241094
m = 0.0529551709
h = 0.5959943932.
(3)
Starting the MATCONT continuer using the above as the initial conditions some points are found as the
saddle and Hopf bifurcations. The following results are found by tracing over the potassium or sodium channel
conductance while holding the other one at its nominal value. In the following table, the bifurcation points
found by the MATCONT software is presented.
Table 1. Bifurcation analysis results derived by the MATCONT software.
Case 1 2 3
Parameter & Prop gK = 20.041 gK = 13.8 gK = 7.89
gNa = 120 gNa = 120 gNa = 120
Eigenvalues
λ1 = −4.61532
λ2 = −0.131416
λ3 = j0.360138
λ4 = −j0.360138
λ1 = −4.7913
λ2 = −0.152529
λ3 = 0.152529
λ4 = 0.425965
λ1 = −7.03698
λ2 = −0.331413
λ3 = 0.331413
λ4 = 1.23876
Equilibrium Points
V = 2.6939079
n = 0.3596422
m = 0.072340538
h = 0.49994929
V = 5.8511661
n = 0.4098447
m = 0.10274964
h = 0.389045
V = 25.621793
n = 0.68536241
m = 0.51701021
h = 0.047222083
Type of Condition Hopf Bifurcation Saddle Point Saddle Point
Case 4 5 6
Parameter & Prop gK = 3.8229 gNa = 210.16 gNa = 311.35163
gNa = 120 gK = 36 gK = 36
Eigenvalues
λ1 = −5.57594
λ2 = −0.422727
λ3 = j1.15984
λ4 = −j1.15984
λ1 = −5.04953
λ2 = −0.125772
λ3 = j0.369568
λ4 = −j0.369568
λ1 = −5.71663
λ2 = −0.137406
λ3 = 0.137406
λ4 = 0.64287
Equilibrium Points
V = 35.333876
n = 0.77395955
m = 0.74063787
h = 0.0186268
V = 0.95529707
n = 0.33241241
m = 0.059205632
h = 0.56236738
V = 2.8902874
n = 0.36274397
m = 0.073974298
h = 0.49288875
Type of Condition Hopf Bifurcation Hopf Bifurcation Saddle Point
556
DORUK: Washout filter based control for the Hodgkin-Huxley nerve cell dynamics,
Case 7 8 9
Parameter & Prop gNa = 370 gK = 36 gNa = 539.014 gK = 36 gNa = 1057.516 gK = 36
Eigenvalues
λ1 = −8.15906
λ2 = −0.176176
λ3 = 0
λ4 = 2.70059
λ1 = −18.43536
λ2 = 2.32620
λ3 = 0.33965
λ4 = −0.33965
λ1 = −26.24354
λ2 = j1.61188
λ3 = −j1.61188
λ4 = −0.42966
Equilibrium Points
V = 8.9984957
n = 0.45979811
m = 0.14295607
h = 0.29038163
V = 26.273054
n = 0.692306
m = 0.534022
h = 0.044102
V = 35.693744
n = 0.776670
m = 0.747283
h = 0.018066
Type of Condition Limit Point Saddle Point Hopf Bifurcation
In order to show the bifurcation characteristics of the model in (1) it will be convenient to show the
bifurcation diagrams obtained by the MATCONT software for the varying values of the potassium and sodium
conductance (gK and gNa). These are given in Figures 1 and 2.
The responses of the equation (1) in the conditions presented in Table 1 are given in the Figures 3–11.
In all cases the external current injection is zero.
100
80
60
40
20
0
-20
V
0 500 1000 150
gNa
100
80
60
40
20
0
-20
V
0 160 180 200
gK
20 40 60 80 100 120 140
Figure 1. The bifurcation diagram for sodium conductance
(gNa) variation.
Figure 2. The bifurcation diagram for sodium conductance
(gNa)variation.
It can be said that the bifurcation analysis results produced by MATCONT package for HH model shares
similar outcome with the research in [18]. This situation can be considered as a verification of the MATCONT
algorithms for a high order nonlinear system.
4. Washout filter theory
The washout filters are linear filtering systems that obey a high pass filter property in which only the transient
signals are allowed (since they have an higher frequency) and the steady state inputs are rejected. Its connection
with the nonlinear systems is an interesting point which leads to the preservation of the equilibrium points of
the original system.
557
Turk J Elec Eng & Comp Sci, Vol.18, No.4, 2010
Mathematically, the washout filter is a high pass system which can be presented in transfer function form
as [16]:
G (s) = Y (s)
X(s) = s
s+d
lim
s→0
G (s) → 0
, (4)
where the constant d is the reciprocal of the time constant of the filter which is also an important parameter
for the application. In the state space form,
˙z = x − dz
y = x − dz
. (5)
0 20 40 60 80 100 120 140 160 180 200
-20
0
20
40
60
80
100
120
time (ms)
V(mV)
0 20 40 60 80 100 120 140 160 180 200
-20
0
20
40
60
80
100
120
time (ms)
V(mV)
Figure 3. Results of the open loop simulation for the case
gK = 20.041, gNa = 120.
Figure 4. Results of the open loop simulation for the case
gK = 13.8, gNa = 120.
0 20 40 60 80 100 120 140 160 180 200
-20
0
20
40
60
80
100
120
time (ms)
V(mV)
0 20 40 60 80 100 120 140 160 180 200
-20
0
20
40
60
80
100
120
time (ms)
V(mV)
Figure 5. Results of the open loop simulation for the case
gK = 7.89, gNa = 120.
Figure 6. Results of the open loop simulation for the case
gK = 3.8229, gNa = 120.
558
DORUK: Washout filter based control for the Hodgkin-Huxley nerve cell dynamics,
0 20 40 60 80 100 120 140 160 180 200
-20
0
20
40
60
80
100
120
time (ms)
V(mV)
0 20 40 60 80 100 120 140 160 180 200
-20
0
20
40
60
80
100
120
time (ms)
V(mV)
Figure 7. Results of the open loop simulation for the case
gK = 36, gNa = 210.16.
Figure 8. Results of the open loop simulation for the case
gK = 36, gNa = 311.35.
0 20 40 60 80 100 120 140 160 180 200
-20
0
20
40
60
80
100
120
time (ms)
V(mV)
0 20 40 60 80 100 120 140 160 180 200
-20
0
20
40
60
80
100
120
time (ms)
V(mV)
Figure 9. Results of the open loop simulation for the case
gK = 36, gNa = 370.
Figure 10. Results of the open loop simulation for the
case gK = 36, gNa = 539.014.
where z is the state variable of the system in (4) and yis the output variable. The application of the above
filter to a nonlinear system like (1) requires that a control gain should be introduced to the filter output. In
this case, the above function will be
˙z = Cx − dz
y = K (Cx − dz)
. (6)
Now the output variable ycan be applied directly to the input of the nonlinear system in consideration. In the
above, the state variable x is now considered as a state vector where the measurable variables are fed back to
the filter through a selector matrix, C ∈ R1×n
.
559
Turk J Elec Eng & Comp Sci, Vol.18, No.4, 2010
0 20 40 60 80 100 120 140 160 180 200
-20
0
20
40
60
80
100
120
time (ms)
V(mV)
Figure 11. Results of the open loop simulation for the case gK = 36, gNa = 1056.516.
When a nonlinear system ˙x = f (x, u) is linearized around a certain equilibrium point the linearized
state equation is (including the higher order remaining terms)
˙x = Ax + Bu + h (x, u) . (7)
Interconnection of the above with (6) yields:
˙x = Ax + Bu + h (x, u)
˙z = Cx − dz
u = K (Cx − dz)
. (8)
The resultant equations present the closed loop controller formed by an output feedback. The closed loop state
equation can be written as
⎡
⎣
˙x
˙z
⎤
⎦
˙xc
=
⎡
⎣
A + BKC −BKd
C −d
⎤
⎦
Ac
⎡
⎣
x
z
⎤
⎦
xc
. (9)
In bifurcation theory [15], it is known that only the quadratic and cubic terms affect the stability of them. So
a possible extension of the above design can be helpful. In this case, the control gain is given to the quadratic
and cubic functions of the output variabley. So in this case (8) will be reshaped as
˙x = Ax + Bu + h (x, u)
˙z = Cx − dz
y = Cx − dz
u = K2y2
+ K3y3
. (10)
This should be a better control system to cope with the bifurcations occurring in especially higher order nonlinear
systems like that of the (1). In the next section the application of the above presentation to the Hodgkin-Huxley
equation is described.
560
DORUK: Washout filter based control for the Hodgkin-Huxley nerve cell dynamics,
5. Washout filter control of the Hodgkin-Huxley equation
For the Hodgkin-Huxley equation in (1) only the membrane voltage variable V is measurable because of that
the controller structure of (10) will be in the following form:
Cm
dV
dt
= −gNam3
h (V − VNa) − gKn4
(V − VK) − gL (V − VL) + Iext
dn
dt
= αn (1 − n) − βnn
dm
dt
= αm (1 − m) − βmm
dh
dt
= αh (1 − h) − βhh
˙z = −V − dz
y = −V − dz
Iext = K2y2
+ K3y3
. (11)
For the sake of computational complexity reduction, the quadratic term gain K2 is taken as zero. The
complete system in (11) is analyzed for possible bifurcations using MATCONT toolbox in MATLAB. The
analysis is performed for each bifurcation case presented in Table 1. As expected, the equilibrium potential V
remains constant throughout the analysis (as washout filter should result). Its value is equal to the equilibrium
value of the respective bifurcation condition (The values given in Table 1). During the trace of the parameter
d additional saddle points are detected by the MATCONT software which are shown in Table 2. The new
bifurcations are due to the interaction of the new eigenvalue brought by the washout filter with the original
eigenvalues in the respective bifurcation condition (i.e. the ones given in Table 1). The trace of the parameter
K3 doesn’t result any critical condition (neither a Hopf and saddle point).
It should be noted here that the value of K3 is kept at 10. For holding d = 0.001and tracing K3 between
-1000 and 1000, MATCONT does not provide any single parameter bifurcation results. The 7th
case of Table 1
(limit point or LP) cannot be solved by MATCONT as it has an eigenvalue exactly at the origin. The algorithm
gets stuck and cannot converge to a solution.
6. Simulation of the Washout Filter – HH model combination
First simulations performed used d = 0.001 and K3 = 10 as the washout filter parameters. The first results
show that this combination presents stable responses for all of the cases presented in Table 1 except the 7th
.
Graphical representations of the associated results are presented in Figures 12–20. In the graphs, both the
membrane potential response and the required electrical stimulation current levels are shown. In addition,
the initial transient behavior of membrane potential and stimulation current are shown in order to prove that
there is no initial jump in the responses and current profile as the original graphs (in the left column) give this
impression to the readers at a first look. To indicate this, the abbreviation IR (initial region) is used next to
the title of the subplots.
561
Turk J Elec Eng & Comp Sci, Vol.18, No.4, 2010
Table 2. The results of the MATCONT analysis for the closed loop.
Case 2 2 3
Value of d d = 0.1518 d = 0.4270 d = 0.3314
Eigenvalues
λ1 = −0.1518
λ2 = −4.79144
λ3 = 0.42696
λ4 = 0.15180
λ5 = −0.15254
λ1 = −0.4270
λ2 = −4.79144
λ3 = 0.4270
λ4 = 0.15180
λ5 = −0.15254
λ1 = −0.3314
λ2 = −7.03698
λ3 = −0.331413
λ4 = 0.331413
λ5 = 1.23876
Type of Condition Saddle Point Saddle Point Saddle Point
Case 6 8
Value of d d = 0.1374 d = 0.3397
Eigenvalues
λ1 = −0.1374
λ2 = −5.71663
λ3 = −0.137406
λ4 = 0.137406
λ5 = 0.64287
λ1 = −0.3397
λ2 = −18.43536
λ3 = 2.32620
λ4 = 0.3397
λ5 = −0.3397
Type of Condition Saddle Point Saddle Point
0 0.5 1 1.5 2
x 10 4
0
1
2
3
V(mV)
Membrane Potential
0 0.5 1 1.5 2
x 10 4
-6
-4
-2
0
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current
0 1 2 3 4 5
-4
-2
0
2
V(mV)
Membrane Potential - IR
0 1 2 3 4 5
-10
-5
0
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current - IR
Figure 12. The response of the closed loop WF-HH system for gK = 20.041, gNa = 120.
562
DORUK: Washout filter based control for the Hodgkin-Huxley nerve cell dynamics,
0 0.5 1 1.5 2
x 10 4
0
2
4
6
V(mV)
Membrane Potential
0 0.5 1 1.5 2
x 10 4
-6
-4
-2
0
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current
0 1 2 3 4 5
-5
0
5
V(mV)
Membrane Potential - IR
0 1 2 3 4 5
-10
-5
0
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current - IR
Figure 13. The response of the closed loop WF-HH system for gK = 13.8, gNa = 120.
0 5 10
x 10 4
0
10
20
30
V(mV)
0 5 10
x 10 4
-15
-10
-5
0
5
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current
0 0.5 1 1.5 2
-5
0
5
V(mV)
Membrane Potential - IRMembrane Potential
0 0.5 1 1.5 2
-15
-10
-5
0
5
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current - IR
Figure 14. The response of the closed loop WF-HH system for gK = 7.89, gNa = 120.
563
Turk J Elec Eng & Comp Sci, Vol.18, No.4, 2010
0 5 10
x 104
0
10
20
30
40
V(mV)
Membrane Potential
0 5 10
x 104
-40
-20
0
20
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current
0 0.5 1 1.5 2
-5
0
5
V(mV)
Membrane Potential - IR
0 0.5 1 1.5 2
-30
-20
-10
0
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current - IR
Figure 15. The response of the closed loop WF-HH system for gK = 3.8229, gNa = 120.
0 2000 4000 6000
0
0.5
1
V(mV)
Membrane Potential
0 2000 4000 6000
-10
-5
0
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current
0 0.5 1
-1
-0.5
0
0.5
V(mV)
Membrane Potential - IR
0 0.5 1
-10
-5
0
5
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current - IR
Figure 16. The response of the closed loop WF-HH system for gK = 36, gNa = 210.16.
564
DORUK: Washout filter based control for the Hodgkin-Huxley nerve cell dynamics,
0 5000 10000 15000
0
1
2
3
V(mV)
Membrane Potential
0 5000 10000 15000
-10
-5
0
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current
0 0.5 1
-1
0
1
2
V(mV)
Membrane Potential - IR
0 0.5 1
-10
-5
0
5
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current - IR
Figure 17. The response of the closed loop WF-HH system for gK = 36, gNa = 31.35163.
0 1 2 3
x 10
4
0
5
10
V(mV)
Membrane Potential
0 1 2 3
x 10
4
-15
-10
-5
0
5
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current
0 0.5 1
0
2
4
6
8
V(mV)
Membrane Potential - IR
0 0.5 1
-10
0
10
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current - IR
Figure 18. The response of the closed loop WF-HH system for gK = 36, gNa = 370.
As it can be understood from Figure 18, Case 7 of Table 1 cannot be stabilized with the current
configuration. This is an expected result as washout filters are not recommended for the models having an
equilibrium point with a complete zero eigenvalue. However, it is interesting to note that a combination of
565
Turk J Elec Eng & Comp Sci, Vol.18, No.4, 2010
d = 0.00001, K3 = 1000 leads to a stable result as shown in Figure 21. Nevertheless, the settling time increases
considerably when compared to the other cases of Table 1 (other than 7). Another observation up to this point
is the increase in the absolute level of the stimulation current with the increasing sodium channel conductance.
0 1 2 3
x 10 4
0
10
20
30
V(mV) Membrane Potential
0 1 2 3
x 10
4
-40
-20
0
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current
0 0.5 1
0
5
10
V(mV)
Membrane Potential - IR
0 0.5 1
-40
-20
0
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current - IR
Figure 19. The response of the closed loop WF-HH system for gK = 36, gNa = 539.014.
0 1 2 3
x 10 4
0
10
20
30
40
V(mV)
Membrane Potential
0 1 2 3
x 10 4
-300
-200
-100
0
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current
0 0.5 1
0
5
10
15
V(mV)
Membrane Potential - IR
0 0.5 1
-300
-200
-100
0
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current - IR
Figure 20. The response of the closed loop WF-HH system for d = 0.00001, K3 = 1000.
566
DORUK: Washout filter based control for the Hodgkin-Huxley nerve cell dynamics,
0 2 4 6 8
x 107
0
5
10
V(mV)
Membrane Potential
0 2 4 6 8
x 10 7
-10
0
10
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current
0 0.5 1
0
2
4
V(mV)
Membrane Potential - IR
0 0.5 1
-10
-5
0
5
time (ms)
Iext(uA/cm2)
Electrical Stimulation Current - IR
Figure 21. Stabilized case 7 of Table 1, with the selection gK = 0.00001, gNa = 1000.
In Table 2, it is stated that the value of d leads to a new saddle point as the additional pole created by
the washout filter interacts with the original eigenvalues. In Figure 22, the membrane potential response of the
closed loop at gNa = 311.351 with d = 0.1374, K3 = 10 is given for convenience. This example corresponds to
the information given in the 4th
column of Table 2.
0 50 100 150 200 250 300 350 400
0
0.5
1
1.5
2
2.5
3
3.5
time (ms)
V(mV)
Figure 22. The closed loop response corresponding to Case 6 of Table 1 with d = 01374, K3 = 10.
As expected the new saddle point is unstable with an oscillation at 40 Hz whereas the original oscillation
corresponding to the information in the 6th
column of Table 1 has a frequency of 59 Hz. The frequencies
are not that much different however the amplitudes are considerably different as easily inferred from Figures
8 and 22. When the value of parameter d is increased the amplitude of oscillations increases. In Figure 23,
567
Turk J Elec Eng & Comp Sci, Vol.18, No.4, 2010
an example response is given for d = 0.5. Note that there is a considerable level of increase in the oscillation
amplitude. With the increasing K3 the oscillation amplitude seems to decrease however the closed loop may
not be stabilized. This effect is demonstrated in Figure 24 as the amplitude drops under 1 mV when K3 is
increased to the value K3 = 1000. Lastly, with a decreasing d, the amplitude of oscillation is decreasing and
after some value the oscillation ceases. As a result, it can be said that washout filter can either stabilize an
equilibrium point (oscillations terminated) or change its amplitude (or frequency). To have a stable result out
of a washout filter application on a HH type nerve fiber one should choose the value of filter time constant d
as small as possible, such as d = 0.001. The value of K3 is then adjusted according to the transient behavior
and level of stimulation current. However, too high values should be avoided so as not to experience solver
difficulties.
0 50 100 150 200 250 300 350 400
0
5
10
15
20
25
time (ms)
V(mV)
0 50 100 150 200 250 300 350 400
0
0.5
1
1.5
2
2.5
3
3.5
time (ms)
V(mV)
Figure 23. The closed loop response corresponding to
Case 6 of Table 1 with gK = 0.05, gNa = 10.
Figure 24. The closed loop response corresponding to
Case 6 of Table 1 with d = 0.5, K3 = 1000.
7. Discussion and conclusion
In this research, a single parameter bifurcation analysis of the oscillating Hodgkin-Huxley type nerve fiber
and a cure for the instabilities are presented. The purpose is to stop the oscillations which correspond to the
repetitive firings in the view of biophysics. First of all, possible single parameter bifurcations are detected using
the MATLAB based MATCONT numerical continuation package. The potassium and sodium channel have
several Hopf and Saddle points at various levels of conductances. The first results show that, the increasing
level of sodium and decreasing level of potassium conductances lead to at least two Hopf bifurcations which are
the main cause of oscillations in the fiber. Not only the Hopf points but also the saddle points can be the causes
of oscillations in the HH type nerve fibers. To control the oscillations a first order washout filter is integrated to
the HH model which will adjust the level of the stimulation current automatically by measuring the membrane
potential. The washout filter is of extended type which interacts with the cubic nonlinearities in HH model
which may be affiliated with the bifurcations. The simulation results show that the time constant of the filter
should be chosen as small as possible to have a stable outcome. Increased level of the filter time constant can
lead to additional saddle points due to the interaction of the filter pole with the eigenvalues corresponding to
the model. In addition, the amplitude of the oscillations is quite sensitive to the value of the filter time constant.
568
DORUK: Washout filter based control for the Hodgkin-Huxley nerve cell dynamics,
With its increasing value, the amplitude of the oscillation increases. The output of the filter is amplified by a
gain which affects especially the duration of the transient behavior. Higher the gain, smaller is the amplitudes
of the oscillations. Gain values too high may lead to solver difficulties during application. Thus, as with the
filter time constant, the designer should select the gain’s value as small as possible. The stimulation current
variation increases with the increasing sodium channel conductance as understood from the graphical results.
In the literature, it is indicated that the washout filters are not suitable for the cases with a complete zero
eigenvalue. In this research, however the HH model is stabilized by the washout filter by selecting an extremely
small filter time constant and a higher gain. This is one of the most interesting results of this work.
A study focused on controlling a nerve fiber can be useful in brain neurophysiologic studies as repetitive
firings correspond to the bifurcations in the basic and advance HH type models. However, the behavior of the
nerve fibers may be different in a realistic situation as the cells operate in a network where each cell affects the
other’s activities. Therefore, bifurcation conditions can be very different. However, this study may constitute a
prototype for similar research on the networked cells. This is one of the topics for future research. In addition
the nonlinear analysis of the washout filter-HH model combination can be another direction for a related future
study.
References
[1] HODGKIN, A.L. and HUXLEY, A.F., 1952. A Quantitative Description Of Membrane Current And Its Application
To Conduction And Excitation In Nerve, Journal Of Physiology-London, 117 (4):, pp. 500-544
[2] WANG, J., CHEN, L. and FEI, X., 2007. Bifurcation control of the Hodgkin-Huxley equations. Chaos, Solitons and
Fractals, 33(1), pp. 217-224
[3] LI, H.-., CHE, Y.-., GAO, H.-., DONG, F. and WANG, J., 2008. Bifurcation analysis of the hodgkin-huxley model
exposed to external DC electric field, IEEE 22nd IEEE International Symposium on Intelligent Control, ISIC 2007.
Part of IEEE Multi-conference on Systems and Control, art. no. 4450897, pp.271-276
[4] JIANG, W., TSANG, K.M. and HUA, Z., 2004. Hopf bifurcation in the Hodgkin-Huxley model exposed to ELF
electrical field. Chaos, Solitons and Fractals, 20(4), pp. 759-764.
[5] MUSTAFA, I.H., IBRAHIM, G., ELKAMEL, A., ELNASHAIE, S.S.E.H. and CHEN, P., 2009. Non-linear feedback
modeling and bifurcation of the acetylcholine neurocycle and its relation to Alzheimer’s and Parkinson’s diseases.
Chemical Engineering Science, 64(1), pp. 69-90
[6] PEREZ VELAZQUEZ, J.L., CORTEZ, M.A., CARTER SNEAD III, O. and WENNBERG, R., 2003. Dynamical
regimes underlying epileptiform events: Role of instabilities and bifurcations in brain activity. Physica D: Nonlinear
Phenomena, 186(3-4), pp. 205-220
[7] HEIDEN, U.A.D., 2006. Schizophrenia as a dynamical disease. Pharmacopsychiatry, 39, pp. S36-S42
[8] SCHMID, G.B., 1991. Chaos theory and schizophrenia: Elementary aspects. Psychopathology, 24(4), pp. 185-198.
[9] GOLDBERGER, A.L., 1988. Nonlinear dynamics, fractals, and sudden cardiac death: New approaches to cardiac
monitoring. Journal of Electrocardiology, 21.
[10] CRAWFORD, J.D., 1991. Introduction to bifurcation theory. Reviews of Modern Physics, 63(4), pp. 991-1037
569
Turk J Elec Eng & Comp Sci, Vol.18, No.4, 2010
[11] BRANDT, M.E. and CHEN, G., 1997. Bifurcation control of two nonlinear models of cardiac activity. IEEE
Transactions on Circuits and Systems I: Fundamental Theory and Applications, 44(10), pp. 1031-1034.
[12] CHEN G, YU X., 1999, On time-delayed feedback control of chaotic dynamical systems. IEEE Trans Circ Syst I
46, pp.767–72
[13] ABED EH, WANG HO, CHEN RC., 1994. Stabilization of period doubling bifurcations and implications for control
of chaos. Physica D, 70, pp. 154–64.
[14] WANG, H.O. and ABED, E.H., 1995. Bifurcation control of a chaotic system. Automatica, 31(9), pp. 1213-1226.
[15] WANG, H.O., CHEN, D.S. and BUSHNELL, L.G., 2000. Dynamic feedback control of bifurcations, 2000, pp.
1619-1624.
[16] HASSOUNEH, M.A., LEE, H.-. and ABED, E.H., 2004. Washout filters in feedback control: Benefits, limitations
and extensions, 2004, pp3950-3955.
[17] DHOOGE, A., GOVAERTS, W. and KUZNETSOV, Y.A., 2003. MATCONT: A MATLAB package for numerical
bifurcation analysis of ODEs. ACM Transactions on Mathematical Software, 29(2), pp. 141-164.
[18] HORNG, T-L., HUANG, M-W., 2006, Spontaneous Oscillations in Hodgkin-Huxley model, Journal of Medical and
Biological Engineering, 26(4), pp. 161 - 168
570