Coupled biological oscillators and animal gaits

Material Information

Coupled biological oscillators and animal gaits
Saffari-Parizi, Mozhdeh
Publication Date:
Physical Description:
xix, 148 leaves : ; 28 cm

Thesis/Dissertation Information

Master's ( Master of Integrated Sciences)
Degree Grantor:
University of Colorado Denver
Degree Divisions:
Department of Mathematics and Statistical Sciences, CU Denver
Degree Disciplines:
Committee Chair:
Tagg, Randall P.
Committee Co-Chair:
Kimbrough, Doris, R.
Committee Members:
Briggs, William L.


Subjects / Keywords:
Harmonic oscillators ( lcsh )
Nonlinear oscillators ( lcsh )
Van der Pol oscillators (Physics) ( lcsh )
Gait in animals -- Mathematical models ( lcsh )
Bifurcation theory ( lcsh )
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )


Includes bibliographical references (leaves 145-148).
General Note:
Integrated Sciences Program
Statement of Responsibility:
by Mozhdeh Saffari-Parizi.

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
57662666 ( OCLC )
LD1190.L584 2004m S33 ( lcc )

Full Text
Mozhdeh Saffari-Parizi
B.S., University of Colorado at Denver, 1993
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Integrated Science
Physics and Mathematics

This thesis for the Master of Integrated Science
degree by
Mozhdeh Saffari-Parizi
has been approved

Saffari-Parizi, Mozhdeh (M.I.S., Master of Integrated Science)
Coupled Biological Oscillators and Animal Gaits
Thesis directed by Associate Professor Randall P. Tagg
Nature is full of rhythmic phenomena. Mathematical models based on cou-
pled oscillators successfully describe many of these behaviors. In this thesis the
fundamental concepts of linear and nonlinear oscillators are first discussed. Af-
ter a review of the physics of oscillators like the simple harmonic oscillator, the
Duffing oscillator and the Van der Pol oscillator, we consider electrical oscillator
circuits. Then examples of biological oscillators are explored. Finally coupled
identical oscillator models for animal gaits are discussed. The dynamics of indi-
vidual oscillators (internal dynamics) in the network is nonlinear and must be at
least two-dimensional to have a Hopf bifurcation. One example of this type of
dynamics is a network of Van der Pol oscillator which is used to produce primary
gaits such as walk, trot, pace and bound. When coupled together, these oscilla-
tors produce patterns of relative phases that correspond to the different animal
gaits. The patterns emerge through bifurcations as parameters describing the
coupling are varied. We will see that at least eight cells are required to model
the quadruped gaits. We conclude with suggestions for practical demonstration
of such coupled oscillator systems.

This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Randall P. Tagg

To my parents, my husband and all of my family.

I would like to thank my husband, Masoud Asadi-Zeydabadi, who helped me
with mathematics and computer programming in this thesis. This thesis would
not have been possible without his help.

Figures .................................................................... x
Tables.................................................................... xix
1. Introduction............................................................. 1
1.1 Coupled Oscillators and Biological Synchronization..................... 1
2. Fundamental Oscillators.................................................. 9
2.1 Simple Harmonic Oscillator............................................. 9
2.1.1 Simple Harmonic Oscillator Without Damping......................... 9
2.1.2 Simple Harmonic Oscillator with Damping........................... 15
2.1.3 Simple Harmonic Oscillator with Periodic Applied Force............ 17
2.2 Duffing Oscillator.................................................... 23
2.3 Van der Pol Oscillator................................................ 33
2.3.1 Hopf Bifurcation and the Limit Cycle............................... 37 Limit Cycle in Van der Pol Oscillator .......................... 43 Hopf Bifurcation in Van der Pol Oscillator......................... 45
2.3.2 Relaxation Oscillator.............................................. 47
3. Electrical Oscillators.................................................. 51
3.1 Resistor Capacitor (RC) Circuits...................................... 51
3.1.1 Charging the Capacitor............................................. 52
3.1.2 Discharging the Capacitor.......................................... 54

3.2 Relaxation Oscillator
3.2.1 Neon Bulb........................................................ 55
3.2.2 Neon Tube Relaxation Circuit..................................... 56
3.2.3 The 555 Relaxation Oscillator.................................... 59
3.3 Oscillation by Gain and Positive Feedback......................... 60
3.3.1 The Wien Bridge Oscillator....................................... 64
4. Oscillators in Biology............................................ 67
4.1 Excitable Neurons ................................................ 67
4.1.1 Hodgkin-Huxley Model............................................. 69
4.1.2 Fitzhugh-Nagumo Model............................................ 73
4.2 Heart and Pacemaker............................................... 86
4.3 Synchronization of Coupled Biological Oscillator (Fireflies)....... 92
5. Coupled Oscillators and Animal Gaits.............................. 95
5.1 Bipedal Gaits...................................................... 96
5.2 Quadrupedal Gaits.................................................. 96
5.3 Central Pattern Generators ........................................ 98
5.4 Symmetries of Animal Gaits........................................ 100
5.5 Two Coupled Oscillators........................................... 101
5.6 Four Coupled Oscillators......................................... 102
5.7 Mathematical Models of Animal Gaits......................... 103
5.7.1 Four Coupled Van Der Pol Oscillators............................ 103
5.7.2 Eight Coupled Oscillators....................................... 105
5.7.3 Eight Coupled FitzHugh-Nagumo Oscillators for Modelling of Pri-
mary Gaits............................................................ 112

5.7.4 Eight Coupled Van Der Pol Oscillators for Modelling of Primary
Gaits ........................................................... 113
5.7.5 Mathematical Model of Primary and Secondary Gaits by Morris-
Lecar Equations ................................................. 121 Primary Gait................................................... 129 Secondary Gait................................................. 130
6. Conclusion, Suggestions and Applications........................... 139
6.1 Simulation of Pulse Coupling ...................................... 139
6.2 Model of Electronic Pulse Coupling................................. 139
6.3 Application of the Gait Simulation in Animals with Birth Defects . 142
6.4 Conclusion......................................................... 143
References.............................................................. 145

2.1 Mass and spring diagram............................................. 10
2.2 Preservation of phase space area.................................... 13
2.3 Time series, phase trajectory and phase portrait of simple harmonic
oscillator: u = l,o;(0) = 1, (^f)t=o = y(0) = 0; (a) and (b) time
series for x and for y, (c) trajectory in phase space, (d) a phase
portrait consisting of several orbits and the fixed point at the origin. 13
2.4 Simple pendulum diagram.............................................. 14
2.5 RLC circuit diagram.................................................. 15
2.6 Simple oscillator with damping: over damped( 72 > u2), critically
damped (y2 u2) and under damped ('y2 < ui2)...................... 17
2.7 Mass and spring with driving force diagram. Force is f = k(x xo)
where xq = Ad cos (fit)............................................ 18
2.8 Simple oscillator with damping and driving force: (a) under damped
(y2 < lo2), (b) critically damped fy2 = u>2) and (c) over damped
(y2 > lu2) with driving force for F = 1............................ 19
2.9 Time series and phase trajectory for simple oscillator with damping
and driving force (under damped)................................... 20

2.10 The ratio of long-time response amplitude to driving amplitude,
A/(^)2, versus driving frequency, Q (top), phase shift, (ip), versus
driving frequency, 0 (bottom), for simple harmonic oscillator with
damping. For to = 1,F = 1,7 = (0.1,0.2)...................... 22
2.11 Time series and phase trajectory for Duffing oscillator with damping:
to = 1,7 = 0.1, e = 0.01 and x(0) = l,x'(0) = v(0) = 0; (a) time
series and (b) phase trajectory.......................................... 25
2.12 Phase trajectory of Duffing oscillator with driving force: u = 1,7 =
0.1, e = 0.01, 0 = a; + e, F = 1 and x(O) = 1,2/(0) = f(0) = 0. ... 26
2.13 Numerical solution of equation (2.46) and the approximate solutions
from equations (2-49) and (2.51), for Duffing oscillator with driving
force: (u = 1,7 = 0.1, e = 0.01,0 = + e, F = 1 and x(0) =
1,^(0) =u(0) = 0).................................................... 28
2.14 Duffing oscillator (weak forcing): a = versus (A = £Lzmj for
Ad = 0.5,u = 1,7 = 0.1, e = 0.01..................................... 30
2.15 Duffing oscillator (Moderate forcing): a = ^ versus (A = Q=^)
for Ad = 1, cu = l,7 = 0.1,e = 0.01.................................. 31
2.16 Duffing oscillator (strong forcing): a = ^ versus (A = ^=sl) for
Ad = 5, to = 1,7 = 0.1, e = 0.01,0 = u + e........................... 32
2.17 Time series (top) and phase trajectory (bottom), critical control pa-
rameter, ec = 0: (left) e = 0.2 < 0, stable spiral and (right)
e = 0.2 > 0, unstable spiral......................................... 41

2.18 The approximate solution by averaging method (equation (2.103))
and numerical solution by Runge-Kutta method for Van der Pol os-
cillator: (e = 0.2 and x(0) = 1, ^'(O) = i'(O) = 0,)...................... 42
2.19 Limit cycle of Van der Pol oscillator: (e = 0.2 and top: x(0) =
l,a;'(0) = i/(0) = 0, bottom: x(0) = ^'(O) = 'u(O) = 2)................. 44
2.20 Hopf bifurcation in Van der Pol equation (time series (top) and phase
trajectory (bottom), critical control parameter, ec = 0: (Left) e =
0.2 < 0, stable spiral and (right) c = 0.2 > 0, stable limit cycle. . 46
2.21 Time series for relaxation oscillator: (a) p = l,a;(0) = 1, ^'(O) =
v(0) = 0, (b) p = 5, x(0) = l,a/(0) = v(0) = 0 and (c) p =
10, x(0) = 1, x'(0) = v(0) = 0............................................ 48
2.22 Limit cycle of relaxation oscillator (initial condition inside the limit
cycle), (a) p = 1, x(0) = lj-X^O) = f/(0) = 0, (b) p == 5, x(0) =
ljX^O) = 1/(0) = 0 and (c) p = 10, a;(0) = 1,£'(0) = i>(0) = 0. ... 49
2.23 Limit cycle of relaxation oscillator (initial condition outside the limit
cycle), (a) p = l,x(0) = l,a/(0) = u(0) = 2, (b) p = 5, x(0)
1. a/'(0) = v(0) = 3 and (c) p 10, x(0) = 1, x'(0) = t/(0) 5. . 50
3.1 Resistor capacitor (RC) circuit............................................ 51
3.2 Charge and current versus time for charging resistor capacitor (RC)
circuit...................................................................... 53
3.3 Charge and current versus time for discharging resistor capacitor
(RC) circuit................................................................. 54
3.4 Neon tube circuit............................................................. 56
3.5 Voltage V versus current I for neon bulb...................................... 57

3.6 Voltage V versus time t for neon bulb.............................. 58
3.7 The 555 oscillator circuit......................................... 60
3.8 Voltage V versus time, t for 555 oscillator........................ 61
3.9 General configuration of a system with gain and feedback........... 62
3.10 A system capable of self-sustained oscillation at frequency to0 where
G{uiq)H{u)q) = 1.................................................. 63
3.11 Wien Bridge oscillator. For comparison to figure (3.10) we note that
y' = y+-y_........................................................ 64
4.1 Schematic diagram of Hodgkin-Huxley model.......................... 70
4.2 x and y Nullclines for Fitzhugh-Nagumo Model: (a = 0.7, b = 0.8, c =
3 and z = 0)...................................................... 75
4.3 Time series and phase trajectory for Fitzhugh-Nagumo model: (a =
0.7, b = 0.8, c = 3 and z = 0).................................... 76
4.4 Time series and phase trajectory for Fitzhugh-Nagumo model: (a =
0.7, b = 0.8, c = 3 and z 0A).................................. 77
4.5 Time series and phase trajectory for Fitzhugh-Nagumo model: (a =
0.7, b = 0.8,c = 3, z = 0.336,) and initial condition (0.25, 0.1). . 78
4.6 Time series and phase trajectory for Fitzhugh-Nagumo model: (a =
0.7, b 0.8,c = 3,z = 0.34) and initial condition (0.6,1.4). .. . 80
4.7 Time series and phase trajectory for Fitzhugh-Nagumo model: (a =
0.7, b = 0.8, c = 3, z = 0.34j and initial condition (0.6, 0.2). ... 81
4.8 Time series and phase trajectory for Fitzhugh-Nagumo model: (a =
0.7, b = 0.8, c = 3, z = 0.34,) and initial condition (0.4,0.7). 82

4.9 Time series and phase trajectory for Fitzhugh-Nagumo Mmodel:
(a = 0.7, b = 0.8, c 3,z = 0.34,) and initial condition (0.6, 0.1). 83
4.10 Time series and phase trajectory for Fitzhugh-Nagumo model: (a =
0.7, b = 0.8, c = 3, z = 0.34,) and initial condition (0.6, 0.15). . 84
4.11 Phase portrait for Fitzhugh-Nagumo model: (a = 0.7,6 = 0.8, c =
3,z 0.34:) and initial conditions (0.6,1.4), (0.6, 0.2), (0.4,0.7)
and (0.6, 0.1)................................................. 85
4.12 Time series (top), phase trajectory (middle) and an enlargement of
the phase trajectory (Bottom) for Fitzhugh-Nagumo model: (a =
0.7,6 = 0.8, c = 3, ^ = -0.25;.................................. 87
4.13 Fitzhugh-Nagumo model with excitation term: (a = 0.7,6 = 0.8, c =
3, z = 0.2b) and E = 6 at time t = 12,16,25,33,38, initial condi-
tion (1,0) and time step dt = 0.02........................... 88
4.14 A typical electrocardiogram (ECG) configuration................. 90
4.15 Fitzhugh-Nagumo model with periodical sinoatrial node signal term:
(a = 0.7,6 = 0.8, c = 3, z = 0.33) (a) without sinoatrial signal,
(b)with sinoatrial signal and (c) Sn(t) = 0.35sin(0.4t)....... 91
5.1 Relative phase for quadrupedal gaits............................. 97
5.2 Two identical coupled oscillator (in phase and out of phase patterns),
(a) in pase pattern, (b) out of phase patterns.................. 101
5.3 Five possible networks for four coupled oscillators. Front is shown
by arrow........................................................ 102
5.4 Four coupled cell network: (a) bidirectional and unidirectional, (b)
bidirectional. Figure (a) and (b) are equivalent................ 104

5.5 Cartoon horse-camel and camel-horse by Hossein Moradi-Khalaj and
Yahya Movahedi-Parizi (top right) horse-camel trot, (top left) horse-
camel pace, (bottom right) Camel-horse trot, (bottom left) camel-
horse pace....................................................
5.6 Zig-zag network for eight coupled cell oscillators.............
5.7 Criss-cross network for eight coupled cells oscillators........
5.8 Model ofpronk by FitzHugh-Nagumo equation: a = 0.02, b = 0.2, c =
0.5, a = 0.01, £ = 0.014,7 = 0.025,5 = 0.02, (27(0), ^(0)) =
(0.06,0.04) and (27(0), y;(0)) = (0,0) for 2 < i < 8..........
5.9 Model of pace by FitzHugh-Nagumo equation: a = 0.02,5 = 0.2, c =
0.44, a = 0.025, (3 = 0.02,7 = -0.01,5 = 0.012, (0), yi(0)) =
(0.06,0.04) and (27(0), ^(0)) = (0,0) for 2 < i < 8...........
5.10 Model of bound by FitzHugh-Nagumo: a = 0.02,5 = 0.2, c
0.44, a = -0.01, P = -0.0102,7 = 0.025,5 = 0.02, (27(0), jft(O)) =
(0.06,0.04) and (27(0), j/i(0)) = (0,0) for 2 < i < 8..............
5.11 Model of trot by FitzHugh-Nagumo equation: a = 0.02,5 = 0.2, c =
0.44, a = -0.02, (3 = -0.002,7 = -0.025,5 = -0.015, (27(0), 2/1(0)) =
(0.06,0.04) and (27(0), y,(0)) = (0,0) for 2 < i < 8............ 117
5.12 Model of jump by FitzHugh-Nagumo equation: a = 0.02,5 = 0.2, c =
0.44, a = 0.02, /? = 0.01,7 = 0.025,5 = 0.015, (27(0), yx(0)) =
(0.06,0.04) and (27(0), 2/i(0)) = (0,0) for 2 < i < 8........... 118
5.13 Model of walk by FitzHugh-Nagumo equation: a = 0.02,5 = 0.2, c =
0.44, a -0.01,/? = 0.0102,7 = -0.025,5 = 0.02, (27(0),yi(0)) =
(0.06,0.04) and (27(0), ^(0)) = (0,0) for 2 < i < 8............. 119

5.14 Model of pronk by Van der Pol equation: p = 0.5, a = 0.02, (3 =
0.02,7 = 0.02,5 = 0.02,£ = r) = 0, (x1(0),y1(0)) = (0.06,-0.04) and
(Xi(0)^(0)) = (0,0) for 2 < i < 8........................................ 122
5.15 Model of pace by Van der Pol equation: p = 0.5, a = 0.02, /3 =
0.02,7 = 0.02,8 = -0.02, £ = p = 0, (rn(0),yi(0)) = (0.06,-0.04)
and (xi(0),yi(0)) = (0,0) for 2 < i < 8.................................. 123
5.16 Model of bound by Van der Pol equation: p = 0.5, ct = 0.02,/? =
-0.02,7 = 0.02, S = 0.02,£ = 77 = 0, (xi(0),j/i(0)) = (0.06,-0.04)
and (Xi(0),yi(0)) = (0,0) for 2 < i < 8.................................. 124
5.17 Model of trot by Van der Pol equation: p = 0.5, a = 0.02,/? =
-0.02,7 = -0.02, <*= -0.02, f = ?7 = 0, (zi(0), 2/i(0)) = (0.06,-0.04)
and (xi(0),yi(0)) = (0,0) for 2 < i < 8. .......................... 125
5.18 Model of jump by Van der Pol equation: p = 0.5, a = 0.02,/? =
-0.02,7 = 0.02,5 = 0.02, £ = r) = 0, (zi(0),2/i(0)) = (0.06,-0.04)
and (rci(0),2/i(0)) = (0,0) for 2 < i < 8................................ 126
5.19 Model of walk by Van der Pol equation: p = 0.5, a = 0.02,/? =
-0.02,7 = -0.02, 5 = -0.02, £ = 7 = 0, (m(0), yi(0)) = (0.06, -0.04)
and (xi(0),yi(0)) = (0,0) for 2 < i < 8.............................. 127
5.20 Model of pronk by Van der Pol equation (phase trajectory): p =
0.5, a = 0.02,/? = 0.02,7 = 0.02,5 = 0.02, £ = r] = 0, (zi(0),2/i(0)) =
(0.06,-0.04) and (xi(0),yi(0)) = (0,0) for 2 < i < 8................. 128

5.21 Model of jump by Van der Pol equation for reference foot of left hind
limb and left fore limb; left panel: (p = 0.5, a = 0.02, (3 = 0.02,7 =
0.02,8 = 0.02, £ = rj = 0) and right panel: (p = 0.5, a = 0.02, f3
0.02,7 = 0.02,5 = 0.02, £ = 7/ = 0). For both left and right panels:
(07(0), 2/1(0)) = (0.06,-0.04) and (2^(0), 2/* (0)) = (0,0) for2 5.22 Model of canter by Morris-Lecar equation: a = 0.17,/? = 0.2,7 =
-0.9,5 = 1, £ = 0,7/ = 0(27(0), 2/i(0)) = (0.4,0.3) and (27(0), 2/*(0))
= (0,0) for 2 < i < 8............................................... 133
5.23 Phase trajectory for model of canter by Morris-Lecar equation: a =
0.17, (3 = -0.2,7 = -0.9,5 = -1,£ = 0,7/ = 0(27 (0), 2/1 (0)) =
(0.4,0.3) and (27(0), 2/i(0)) = (0,0) for 2 < i < 8. ............... 134
5.24 Model of transverse gallop by Morris-Lecar equation: a = 0.09,/? =
-0.8,7 = -0.3,5 = -0.08, £ = 0,t/ = 0(27(0), 2/1(0)) = (0.061,0.047)
and (21^(0), 2/i(0)) = (0,0) for 2 < i < 8.......................... 135
5.25 Phase trajectory for model of transverse gallop by Morris-Lecar equa-
tion: a = 0,09,13 = 0.8,7 = 0.3,5 = 0.08, £ = 0,7/ = 0
(27(0), 2/i(0)) = (0.061, 0.047) and (xi(0),yi(0)) = (0, 0)for2 8................................................................... 136
5.26 Model of rotary gallop by Morris-Lecar equation: a = 0.02,/? =
0-23,7 = 0.4,5 = 0.7,£ = 0,7/ = 0, (27(0),2/i(0)) = (-0.79, O.48),
(22(0), 2/2(0)) = (0, 0.56) and (2^(0), 2/i(0)) = (0,0) for 3 < i < 8. . 137

5.27 Phase trajectory for model of rotary gallop by Morris-Lecar equation:
a = 0.02,0 = 0.23,7 = -0.4,5 = -0.7, f = 0, rj = 0, (i(0), yi(0)) =
(-0.79,0.48), (£2(0), 2/2(0)) = (0,0.56) and (xi(0), yi{0)) = (0,0) for
3 < i < 8.................................................................... 138
6.1 555 circuit as unit cell in the model of pulse coupled oscillator. . 140
6.2 The model of pulse coupled oscillators. Each square represents a
circuit like the one shown in figure (6.1)............................. 141
6.3 Two views of an animal with three legs: Museum of Natural History
and Technology, Shiraz University, Iran, (by Samira Moradi-Khalaj). 142

5.1 Possible assignments of cell to legs in zig-zag network....... 109
5.2 Possible assignments of cell to legs in criss-cross network... 109
5.3 The phase shift for primary gaits in term of xi...................... 110
5.4 Symmetry of primary gaits............................................ 110
5.5 Parameter values of primary gaits in Fitzhugh-Nagumo model. ... 113
5.6 Parameter values of primary gaits in Van der Pol model............... 121
5.7 Parameter in the case of linear synaptic coupling leading to stable
Hopf bifurcations................................................... 130
5.8 Parameters of primary gait for synaptic coupling Morris-Lecar model:
£ = ?7 = 0................................................... 131
5.9 Secondary gaits with linear diffusive coupling in Morris-Lecar. . . . 131
5.10 Secondary gaits with linear synaptic coupling in Morris-Lecar. . . . 132
5.11 Nonzero initial conditions in Morris-Lecar........................... 132

1. Introduction
1.1 Coupled Oscillators and Biological Synchronization
The purpose of this thesis is to guide the selection of experimental oscil-
lators that will be used to create arrays of coupled dynamical systems whose
behavior is governed by symmetries in the couplings. In particular, the thesis is
motivated by recent work on animal gaits. Transitions to different patterns of
footfalls can be understood on broad grounds by the effect of symmetry on an
oscillator network modelling the neurophysiology of the animals. An experimen-
tal robotic horse is under construction whose leg motions will be orchestrated by
coupled electronic or fluidic oscillators. The choice of oscillators for this robot
will be governed by several considerations: experimental realizability, precision
of control of the couplings, the availability of matching mathematical models,
and potential for comparison to biological oscillators. These critical factors will
be examined in detail for various oscillator models in this thesis.
Following this introduction, chapter 2 concentrates on oscillator fundamen-
tals and mathematical prototypes in order to develop the language for describing
and evaluating the dynamics of individual oscillators. Chapter 3 then investi-
gates specific electronic circuits that might conveniently be used in coupled-
oscillator experiments. In chapter 4, there is a shift of emphasis in order to
describe some of the more basic theoretical models for biological oscillators.
This gives a perspective for experiments, even if a specific circuit implemen-
tation does not attempt to mimic most of the biological details. Indeed, one

expectation to come from the study of nonlinear dynamics is that the overall
behavior is independent of many of the details and depends instead on such
broad features as symmetry and the number of elements interacting with one
another. Chapter 5 continues with the biological orientation by describing the
animal gait problem and some corresponding oscillator models. Concludes by
outlining principal features one would like to mimic in an electronic coupled-
oscillator model of central pattern generators.
First, however, we give an overview of the coupled oscillator problem as it
arises in a broad variety of physical and biological systems. In the seventeenth
century, Huygens invented the pendulum clock [31]. He watched two clocks that
were hanging side by side. He noticed that the two pendulums were swinging
in perfect synchrony. They swung for hours without any broken steps. Then he
tried disturbing them and within half an hour they regained synchrony. Huygens
believed that the clocks must somehow influence each other by air movement or
vibrations in their common support. In other observations he moved the clocks
to opposite sides of the room. The clocks gradually fell out of step, one losing five
seconds a day relative to the other. His observations initiated the mathematics
of the theory of coupled oscillators [38]. Recently this investigation has been
reproduced using modern laboratory equipment [4].
Coupled oscillators can be found in nature. Examples are the pacemaker
cells in the heart, insulin-secreting cells in the pancreas, and the neural networks
in the brain and spinal cord that control the rhythmic behaviors of breathing,
running, and chewing. Another example is the synchronization of congregations
of fireflies [38]. In each case, an element isolated from the system can continue

to oscillate in its own rhythm. This rhythm is entirely dependent on proper-
ties of the elemental subsystems itself. Typically, the elemental oscillators are
autonomous nonlinear dynamical systems, meaning that their behavior is not
determined by an external time-dependent forcing. In fact, the systems are ca-
pable of self-sustained oscillation. For example one can isolate a firefly from
other insects, put it into a place with a constant temperature and illuminance
and observe that the lone firefly produces rhythmic flashes. The systems that
generate these rhythms are dissipative and thus need some energy source to sus-
tain their oscillations. By contrast, ah example of non self-sustained oscillation
is the tidal flow caused by the daily variation of gravitational forces due to the
Moon. These oscillations would not happened if there were no Moon to provide
external forcing [31].
One class of natural oscillator is called a relaxation oscillator. In rhythmic
fashion, some quantity slowly grows to a threshold level and then suddenly resets
to a lower value. Then the process repeats. This behavior occurs, for example,
in spontaneously firing neurons. If a current is injected into the cell, the electric
potential slowly changes until the cell potential reaches a threshold of about
50mV. At this point the electric potential suddenly spikes for a duration of
about two milliseconds. After discharging, the cell resets to about 70rnV.
Such integrate and fire oscillators are very important in various problems of
neuroscience [31].
Another example is an experimental situation in which an isolated heart
preserves its ability to rhythmically contract in vitro, i.e. outside the body and
decoupled from the bodys nervous system. If the heart is arterially perfused

with a physiologic oxygenated solution, the activity of the primary pacemaker,
the sinoatrial node, continues to trigger the beats. The pacemaker is composed of
a population of rhythm-generating cells that synchronize their action potentials.
This system is different from the heart in vivo because the frequency of the heart
in vitro is higher and the variability of the interbeat intervals is lower than those
of the heart in vivo. Thus the heart is a highly complex system that can behave
as a self-sustained oscillator. In the body, the normal heartbeat is influenced
by the the whole cardiovascular system, including all the control loops of the
autonomic neural system and the brain centers involved in the regulation of the
heart rate. This system is also affected by other physiological rhythms, such
as respiration. Nonetheless, we see that the heart is capable of producing its
rhythm by itself in spite of these perturbations [31].
Peskin worked on a schematic model of the hearts pacemaker (10,000 cells
of the sinoatrial node) in order to understand how they synchronize their indi-
vidual electrical rhythms to generate a heartbeat [30]. He suggested that the
pacemaker is similar to a large number of coupled identical oscillators In the
model each oscillator is an electrical circuit with a capacitor in parallel with a
resistor. A constant current causes the voltage across the capacitor to increase
steadily. The voltage reaches a threshold where the capacitor discharges, and
the voltage drops to zero (firing and discharging of the pacemaker). In the
model, each oscillator affects the others only when it fires by causing a fixed
jump in the voltage of every oscillator coupled to the one that fired. Peskin
stated that the system would always eventually become synchronized, and it
would synchronize even if the oscillators were not quite identical. Peskin mod-

elled only two coupled oscillators. He used a variation of Poincares concept of
a mathematical equivalent to stroboscopic photography. Using the pulse of one
oscillator (A) as a strobe, he took the other oscillators (Bs) voltage every time
A fired. Peskin found a formula for the change in Bs voltage between snapshots.
Based on the formula, if the voltage is less than a certain critical value, it will
decrease until it reaches zero, whereas if it is larger, it will increase. Eventually
B will synchronize with A. If Bs voltage is almost equal to the critical voltage,
then it stays poised at criticality and the oscillators fire about half a cycle out
of phase from each other. This equilibrium is unstable however. Note that a
general observation seen in both Peskins model and in Huygens experiments is
the following: When two identical oscillators are coupled, there are two possi-
bilities: synchrony with a phase difference Of zero, and anti-synchrony with a
phase difference of one half period [38].
Natural oscillators may be coupled only to a few immediate neighbors or they
may be coupled to all the oscillators in a larger community. In such systems of
coupled oscillators, synchrony is the most familiar mode of organization. For ex-
ample, thousands of male fireflies gather in trees at night and flash on and off in
unison to attract females. Fireflies are a paradigm of a pulse coupled oscillator
system. They interact only when one sees the sudden flash of another and shifts
its rhythm accordingly. In models of such behavior Strogatz proved that oscil-
lators started at different times will always become synchronized [25],[38],[36].
Inspired by Peskins work, Strogatz wrote a computer program to do nu-
merical experiments with a large number of pulse-coupled oscillators. Mirollo
and Strogatz proved mathematically that the system always becomes synchro-

nized for any number of oscillators and for almost all initial conditions. If one
oscillator kicks another over threshold, they will remain synchronized forever
In contrast to the firefly behavior bioluminescent algae show desynchroniza-
tion when removed from daily light variations. J. Woodland Hastings and his
colleagues at Harvard University have found that when a tank full of Gonyaulax
is kept in constant dim light, a circadian glow rhythm is observed with a period
close to 23 hours. As time goes by, however, the rhythm gradually damps out.
The individual cells continue to oscillate, but they drift out of phase because of
differences in their natural frequencies [38].
Do we account for synchronization of some systems and desynchronization
in others? Consider each element in a system to be an oscillator with a preferred
amplitude and frequency of oscillation. This preferred behavior, to which a lone
oscillator will return even if it is momentarily perturbed, is called a limit cycle.
Winfree pointed out that if oscillators are weakly coupled, they remain close to
their limit cycles at all times. Thus amplitude variations could be neglected and
only the spread or decrease of relative phases need be considered. He found that
the systems behavior depends on the width of the frequency distribution. If the
spread of frequencies is large compared with coupling, the system always lapses
into incoherence, as if it were not coupled at all. But if the frequency spread
decreases below a critical value, part of the system spontaneously freezes into
synchrony [44].
Finally, a brief mention is made of the animal gait problem here, with more
detail to come in chapter 5. Quadruped gaits follow patterns of four oscillator

systems. They illustrate how the symmetry of a system oscillating in phase can
be broken. For example when a rabbit bounds, it moves its front legs together
and then its back legs. There is a phase difference of zero between the two front
legs and of one half of a cycle between the front and back legs. In a giraffes
pace the front and rear legs on each side move together. In a horses trot the
diagonal legs move in synchrony. When an elephant ambles it lifts each food in
turn, with phase differences of one quarter at each stage. When gazelles pronk,
all legs move in synchrony. Biologists believe that the nervous system must
produce the different patterns of locomotion. They hypothesize the existence
of networks of coupled neurons called central pattern generators, which act as
oscillators that control the gait more or less independently from signals from the
brain [34], [19].
Most animals have several gaits. For example horses walk, trot, canter and
gallop. Mathematical models show that the same central pattern generator
circuit can produce all of an animals gaits. Only the strength of the couplings
among neural oscillators must vary.
An interesting question concerns 3-legged gaits. When three identical os-
cillators are coupled in a ring they can be phase-locked in four basic patterns.
(1) three in synchrony, (2) each one third out of phase with the others, (3) two
in synchrony and one with random phase shift and (4) two half a cycle out of
phase and one twice as fast. An example of this last case is a person using a
walking stick [38].
The variety of coupled oscillators systems discussed above suggests that a
rich phenomenology should be achieved with laboratory coupled oscillator cir-

cuits. By incorporating significant features of the models of various biological
systems, one hopes to mimic the spatial and temporal patterns discussed above.
This may lead to new insights into these models and to the corresponding biolog-
ical systems. Laboratory systems may also lead to novel methods of preparing
and controlling coupled networks and to the discovery of dynamics whose bio-
logical counterparts are yet to be identified.

2. Fundamental Oscillators
In this chapter we will study the mathematics of several fundamental types
of oscillators. We will start with the simple harmonic oscillator as the basis
of several key ideas. Most notable about the simple harmonic oscillator is its
linearity: if x (t) and x2(t) represent solutions that describe the oscillators po-
sition as a function of time, then an arbitrary combination AiX]_{t) + A2x2(t) is
also a solution. While suitable for modelling many types of physical systems,
especially if their oscillations are small, the linear simple harmonic oscillator
model must be extended to nonlinear models in order to describe a much wider
range of phenomena. Fundamental mathematical models for nonlinear oscilla-
tors include the Duffing and Van der Pol oscillators, whose behavior will also be
described in this chapter.
2.1 Simple Harmonic Oscillator
The motion of a mass on spring, figure (2.1), is one example of the simple
harmonic oscillator. The motion of the simple pendulum in the approximation
of small angle is another example.
2.1.1 Simple Harmonic Oscillator Without Damping
The force on a mass attached to a spring is proportional to the displacement
x from the equilibrium point but acts in the opposite direction. We have
fa = ~kx, (2.1)
where k is the spring constant. From Newtons second law we have
m = kx

3 ^
x = 0
Figure 2.1: Mass and spring diagram.
where m is the mass of the object. Define the angular frequency as
uj = \l
Then we have the equation of a simple harmonic oscillator
+ U! X = 0.
The solution of this equation is sinusoidal with a period of
T =
in the general form of
x = A sin tot + B cos u>t
= C sin (not + ip).
The constants (A and B) or (C and ip) are determined by initial conditions.

Further insight comes from re-examining the dynamics from a geometric
viewpoint. The second order equation (2.4) is equivalent to the system of two
first order equations
y (2.7)
UJ2X (2.8)
These describe motion of a state (x,y) in an abstract plane called the phase
plane. Plotting y versus x for any initial condition in the phase plane makes a
curve that is called a trajectory. An accumulation of trajectories for a variety of
initial conditions provides a comprehensive view of the dynamics that is called
a phase portrait.
Note that (^, describes a vector that is tangent to the trajectory at
point (x,y). An accumulation of such vectors is similar to the velocity field of
a moving fluid: here the fluid consists of points in the phase plane. Linking the
vectors gives streamlines that correspond to trajectories. The mapping of points
of the phase plane onto other points in the same plane is called the flow of this
dynamical system.
The fixed (equilibrium) point of this system can be determined by putting
the right side of the equation (2.7) and (2.8) equal zero. The only fixed point
for this system is the origin (0,0).
All other points belong to trajectories that form closed curves This can be
seen as follows. From equation (2.6) and (2.7) we have

Eliminating time between (2.6) and (2.9) by summing suitable multiples of
the squares of each side of the equations gives
c2 (Cu>y
= sin 2(ut + ip) + cos2(a)t + (p)
= 1.
Thus the phase space (trajectory) for a simple harmonic oscillator is a closed
orbit which is an ellipse. Its equation given above can be rewritten as
^4-^ = 1
where a = C and b = Cu. The center of the ellipse is at the fixed point which
in this case is the origin and the length of its axes along the x and y axes are 2a
and 2b respectively. Each orbit represents a constant total energy. This means
the system of a mass and spring without damping is a conservative system. One
of the important properties of a phase portrait in a conservative system is area
preservation. This means that all of the points in an area, at time t\ occupy
an equal area, A2, at time t2 (figure (2.2)). Another property of phase portraits
of any kind of system (not just conservative systems) is that any trajectory never
crosses itself. This is the result of uniqueness of the trajectory in a deterministic
dynamical system.
Suppose the initial condition is (:r(0), y(0)) = (1,0). Then the constants in
equation (2.6) can be determined to be (A = 0,B = 1) or (C l,(p = 7r/2).
Let the natural frequency lj = 1. Then the trajectory is a circle with radius of
1 and its center is located at the origin. Figure (2.3) shows the time series and
phase space trajectory of this example.
The dynamics of motion of a pendulum is another example of the funda-

Figure 2.2: Preservation of phase space area.
Figure 2.3: Time series, phase trajectory and phase portrait of simple har-
monic oscillator: u> = l,o;(0) = 1, (fjr)t=o = y(0) = 0; (a) and (b) time series
for x and for y, (c) trajectory in phase space, (d) a phase portrait consisting of
several orbits ana the fixed point at the origin.

Figure 2.4: Simple pendulum diagram.
mental oscillators. The equation of a simple pendulum is given by
d29 o.
- + ur. sm 6 = 0
where 0 is the angular position and the angular frequency is
Here g is the acceleration of gravity and l is the length of the pendulum (figure
(2.4)). For small angles we have
Thus in the small angle approximation the equation of a pendulum becomes a
simple harmonic equation
j2 a
^+u29 = 0. (2.15)
A key point for both types of linear oscillators is that the frequency does not
depend on the amplitude of the oscillation.

2.1.2 Simple Harmonic Oscillator with Damping
Consider a third prototype for the harmonic oscillator: an RLC circuit
composed of a resistor R, capacitor C, and inductor L, arranged in a loop figure
(2.5). The charge on the capacitor is q and the current flowing in the circuit is
Figure 2.5: RLC circuit diagram.
i. Applying Kirchoffs Law of voltages to the loop, we have
Q di
iR L = 0.
C dt
Now the current i around the loop is related to the charge q on the capacitor by
Thus, differentiating equation (2.16) once gives
i % _ _ L*i = o.
C dt dt dt2
Substituting for ^ from equation (2.17) gives
-1. % R----------L = 0.
C dt dt2

Rearranging terms gives
d2i Rdi 1 .
dt? + TJt + lc1 =
Let the variable i be replaced with x, define 2y = ^, and identify u2 = to
d2x dx 9
dF + 2T*+1JI = 0' (2'21)
This is the equation for a simple harmonic oscillator with damping. The angular
frequency u> is the natural frequency of the simple harmonic oscillator. The
second term, 2y^|, is a damping term. In the solution of equation (2.21) we
have three cases [32].
The first case is known as over-damped motion (y2 > u2) and its solution is
x = e-7* (Aeat + Be~at) where a = a/t2 ^2> (2.22)
and where A and B are constants that can be determined by initial conditions.
The second case is critically damped motion (y2 = uj2). The solution of
equation (2.21) in this case is
x = e~'rt(A + Bt), (2.23)
where A and B can found from initial conditions.
The last case is called under-damped or damped oscillatory motion (y2 <
to2). In this case equation (2.21) has the solution of
x = e7t (A sin pt + B cos fit) with fi = \/uj2 y2. (2.24)
Again, the constants A and B are determined by initial conditions. Equation
(2.24) can be written as
x = Ce~jt cos (fit ip) with fi = \/u2 ~ y2, (2.25)

where C = y/A2 + B2 is the amplitude and

Figure 2.6: Simple oscillator with damping: over damped(72 > oj2), critically
damped ('y2 = to2) and under damped (y2 < to2).
and 2 damping is large enough so that there is no oscillation. In the last case,
damping is not as strong as in the two other cases so there are oscillations about
the equilibrium. However the amplitude decreases with time. See figure (2.6)
for plots of displacement, x, versus time, t, for all three cases. (In the circuit
model, x actually represents the current i).
2.1.3 Simple Harmonic Oscillator with Periodic Applied Force
Now suppose a periodic external force is applied to the simple harmonic
oscillator. One way of doing this is shown in figure (2.7). The equation of
motion becomes
+ 2777 +w2x = Fcos(fi£), (2.26)
dV dt

Figure 2.7: Mass and spring with driving force diagram. Force is f = k(x xo)
where xo = Aa cos(£lt).
where is the angular frequency of the applied force and F is a force per unit
mass (thus having units of acceleration). In the set up shown in figure (2.7)
The general solution to (2.26) consists of two parts: a transient that is a
solution of the unforced equation (2.21) plus a particular solution to the forced
problem. The transient to (2.26) depends on the magnitude of j2 compared to
u2 according to the three cases discussed above. Note that in all three cases the
transient solution goes to zero after a sufficient period of time. The particular
solution of (2.26) is given by [32]
x = A(Q) cos(Qt where the amplitude
y/(Q2 u2)2 + 472^2
tan 0 = rpry 0 < (f) < 7T. (2.29)
The particular solution is known as the steady-state solution because after a
short time the motion of forced oscillation is described by this form. Figure (2.8)
shows that as time goes on a solution starting from arbitrary initial conditions

converges to the steady-state (particular) solution. A time series and phase
trajectory for a simple oscillator with damping and driving force for the under-
damped case is shown in figure (2.9).
Notice in equation (2.29) that the phase 0 for small Q, passes through
7t/2, where = tv, and approaches 7r as Q. becomes much larger than u>. The
Figure 2.8: Simple oscillator with damping and driving force: (a) under
damped (j2 < u2), (b) critically damped (q2 = u2) and (c) over damped
(72 > oj2) with driving force for F = 1.
maximum amplitude occurs at the resonant drive frequency (see (2.10))
aR = Ji-2(Tf

Figure 2.9: Time series and phase trajectory for simple oscillator with damping
and driving force (under damped).

A common practice is to define a parameter Q by
Q s <2-31>
An oscillator whose damping coefficient 7 is small relative to the natural fre-
quency u> is said to be a high-Q oscillator. In terms of this parameter,
nR = aV1~ 2Q2' (2'32)
For high-Q oscillators, the resonant frequency Qr is at the natural frequency
ClR ~ w (Q > 1). (2.33)
On the other hand, a strongly damped system has no maximum response am-
plitude when
Q < A- (2.34)
The amplitude at resonance is given by
and thus, for large Q,
A~Q(5) (g>i}- (2-36)
An examination of figure (2.10) shows that there is a peak in the plot of ampli-
tude A versus drive frequency fh Usually, the energy of an oscillator is propor-
tional to the square of the amplitude. We measure the width of the peak by the
spread in frequency (Q+ fi!_) where
42(Si+) = = \a\. (2.37)

For large Q oscillators, this full width at half maximum energy is given by
(0>>1) (2'38)
Finally, the Q of the oscillator also measures the rate of energy loss due to the
Figure 2.10: The ratio of long-time response amplitude to driving ampli-
tude, Aj(£)2, versus driving frequency, 0 (top), phase shift, ( ing frequency, £1 (bottom), for simple harmonic oscillator with damping. For
u = 1,F= 1,7= (0.1,0.2).
dissipative term per radian of oscillation (i.e. for a time interval At such that
Q.A t = 1):
Average stored energy
Energy lost per radian of oscillation
= Q.

2.2 Duffing Oscillator
Now we consider nonlinear oscillator models. An equation of a system of a
mass and nonlinear spring without damping is the Duffing equation
d2x 3
mz- kx kx ,
where m is the mass and k and k are the linear and nonlinear spring constants
respectively. We may view this as a spring that, for k > 0, becomes stiffer as
the spring is stretched. Correspondingly, for k < 0, the equation represents a
spring that becomes weaker as it is stretched. Dividing this equation by m we
+ urx + ex3 = 0
where w = y | is the natural frequency of the linear system of mass and spring
when, k = 0, and e = ^ presents the nonlinearity of the system. As a physical
example of a Duffing oscillator, consider the pendulum equation (2.12), +
uj2 sin 9 = 0. The Taylor expansion of sin 6 is
e3 e5 e7
SinW = e-3! + 5r7! +
Using the Taylor expansion of sin 0 up to the order 93 in the pendulum equation,
we get
= (2.43)
This is a special case of the Duffing equation with a soft spring. In the case of
damping the general Duffing equation is
d2x ^ dx 9 o
r + 27 + UJ2X + ex3 = 0.
dV dt

A time series and phase trajectory of the Duffing equation with damping is given
in figure (2.11).
The Duffing equation (2.41) can be shown to have a periodic solution x =
Acos(un(A)t). Here the nonlinear frequency uin depends on the amplitude A.
For small e, Poincare-Lindstedt method [37], [23], give the following relationship
between frequency and amplitude
OJn(A) ~ OJ +
3 eA2
As e 0, the nonlinear term ex3 disappears and the frequency of the periodic
solution returns to the amplitude independent value oj.
Adding a periodic external force to this system gives:
+ 2q^ + u?x + ex3 = F cos (fit)
at2 at
A phase trajectory of the Duffing oscillator with damping and periodic driving
force is shown in figure (2.12). Several different approximation methods have
been developed to describe solutions to the Duffing equation that have the same
period 2tt/£1 as the forcing term. See, for example, Stoker, Hagedorn, and
Nayfeh, [35], [18], [28].
Atlee Jackson, [23], describes one approach where the nonlinear ex3 term is
dropped but the natural frequency to is replaced by the amplitude dependent
frequency un(A) given by equation (2.45). The resulting linear forced oscillator
equation is
^ = F*m- (2-47)
We assume a solution
x(t) = Acos(flt 0), (2.48)

Figure 2.11: Time series and phase trajectory for Duffing oscillator with damp-
ing: u> = 1,7 = 0.1, e = 0.01 and :r(0) = 1, x'(0) = ^(O) = 0; (a) time series and
(b) phase trajectory.

Figure 2.12: Phase trajectory of Duffing oscillator with driving force: ui
1,7 = 0.1, e = 0.01, D = l> + e, F = 1 and x(0) 1, ^'(O) = w(0) = 0.

where A is the same amplitude in the expression for uin(A). Then from (2.28)
the steady-state solution of equation (2.47) is given by
x =
cos (Qt (ft).
y/H,{A) O2)2 + 472D2
This solution is an approximate solution of the Duffing equation. Self-consistency
requires that
A = =. (2.50)
((" + T?)2 2) +VA1
Squaring both sides and rearranging gives
w+^V-S22) +472Q2] -F2 = 0.
Roots of this equation give the approximate amplitude. Solving equation (2.51)
numerically for A and equation (2.29) for (f> then provides the approximate solu-
tion for the Duffing equation with driving force. The solution can be determined
numerically by the Runge-Kutta method for integrating differential equations
like equation (2.46). Figure (2.13) gives the comparison of the numerical solu-
tion of equation (2.46) and the approximate solutions from equations (2.49) and
(2.51). Let us assume that the frequency of driving force Q is near the natural
n = u + eA. (2.52)
Then Q2 up to first order in e is
O2 = (cu + eA)2
u>2 + 2ueA. (2.53)

Figure 2.13: Numerical solution of equation (2-46) and the approximate solu-
tions from equations (2-49) and (2.51), for Duffing oscillator with driving force:
(uj = 1,7 = 0.1, e = 0.01, Q, = cu + e,F=l and a;(0) = 1, £'(0) = i>(0) = 0).

Using equation (2.45) and again keeping terms to first order in e, we have
o/ .. 2 3 A?
w + "g^"2we-
U(A) n2) k (|i) a)(2^)=
Using these approximations in equation (2.51), we obtain a cubic equation for
-----(2u;e)2 + 4q2(a; + eA)2
= F2.
By letting a = 3A2/(8u>), and P = 3F2/(32e2w3) then we have
(a A)2 -j-
2 / 2q(a; + eA)N 2
0 = 0.
In this equation, a measures the square of the amplitude of oscillation, A is
proportional to the difference of the driving frequency and natural frequency
(f2 u>) and f3 measures the square of the applied force. Figures (2.14), (2.15)
and (2.16) show a versus A for small ft (F = 0.5), moderate (5 (F 1) and large
(3 (F = 5). For small /?, no nonlinearity shows and the maximum amplitude
happens at A = 0 where the driving frequency is equal to natural frequency
(U = uj). For moderate P, nonlinearity causes the maximum to occur at larger
value of A. For large P, we observe a bent turning curve as a result of three
real roots of equation (2.57) for a for a certain range of A. In the case of
multiple roots, the solution for the largest and the smallest (amplitude) a are
stable and the middle value for the amplitude is not stable [23]. A key result
of this bistability of amplitudes (i.e. region where there are two different stable
solution) is the phenomena of hysteresis and amplitude jumps. If the forcing

Figure 2.14: Duffing oscillator (weak forcing): a = ^ versus (A = -) for
Ad = 0.5,a; = 1,7 = 0.1, e = 0.01.

Figure 2.15: Duffing oscillator (Moderate forcing): a = versus (A = )
for Ad = 1,uj = 1,7 = 0.1, e = 0.01.

Figure 2.16: Duffing oscillator (strong forcing): a = versus (A = Q-j£)
for Ad = 5, u = 1,7 = 0.1, e = 0.01, = u> + e.

frequency £1 starts out well below to (A is negative) and is steadily increased,
the amplitude climbs to a value at point a in (2.16). Further increase in Cl
causes a sudden jump to a smaller amplitude at point b and then the amplitude
continuously decreases further as frequency Cl grows. Next, if Cl is decreased
from a value much larger than to the amplitude grows until a point c is reached
which is to the left of B. Here a further decrease in to causes the amplitude
to jump to point d. The key feature is that the downward and upward jumps
occur at different frequencies and there is a difference in solution paths for Cl
increasing and D decreasing. This is what is meant by hysteresis.
In summary, the Duffing oscillator can be considered a natural extension
of the simple harmonic oscillator to systems where the restoring force is no
longer linear in displacement. Key new features are an amplitude dependent
frequency of oscillation for the unforced system and a bent tuning curve for the
forced system. In fact, other phenomena, including solutions with three times
the drive period (called subharmonics) and chaos can occur for this system. See
Strogatz [37] for further discussion.
2.3 Van der Pol Oscillator
We saw with the Duffing oscillator that new features were introduced by
nonlinearity of the restoring force. Another modification would be to make the
damping term nonlinear. The way this might be done is not as obvious as
the nonlinearity introduced by stiffening or weakening springs in the Duffing
Suppose we think of the oscillator as a particle of mass m moving in a

potential V(x). The restoring force Fr in the absence of damping is given by
dV, .
Fr = ~^x)-
For the simple harmonic oscillator,
V(x) =
The mechanical energy of the system is
E = -rax1 + V(x),
where x = ^. The time rate of change of this mechanical energy is
dE d /1 ,2 ... A
dT = 7t{ 2mx +v(x))
1 dV .
= -m(2xx) + -zx
2 dx
Now we can think of this rate of energy change as a function of the state
(x, x) of the system:
= G(x, x)
For symmetric oscillators it would be reasonable to expect this function to be
even in both x and x. The undamped simple harmonic oscillator corresponds
to the trivial case
G(x,x) = 0, (2.63)
while the simply damped system equation (2.21) corresponds to
G(x,x) = 2m'yx1.

To see this, recall that for (2.21), where Fr = kx = muFx,
mx Fr = 2m^x
Then equation (2.61) gives
dE , ...
= {2m'yx)x
G(x,x) = 2m'yx2.
Here we assume 7 > 0. This energy dissipation is independent of position,
preserving the linearity of the system. Another kind of energy loss due to dry
friction is described by
G(x, x) = fimg sign(i)i
1, x < 0
< 0, x = 0
1, x > 0
This dissipation is independent of position, but introduces nonlinearity into the
system via the function sign(x). The equation of motion is
mx + = fJ,mg sign(x). (2.70)
Further details on this dynamics may be found in [24].
Now let us consider extending the general dissipation term to include posi-
tion dependence
G(x, x) = (C + C'x2)2m/yx2, (2-71)

where C and C' are two constants. Equation (2.64) corresponds to C = 1
and C' = 0. Suppose instead that we have some means of putting energy into
the system (negative dissipation) when x is small enough, but that this energy
exchange reverses once x exceeds some threshold x = a. We describe this by
G(x, x) = ^1 ^ ^ 2m/yx2
Now for x < a and 7 > 0, G is positive and thus ^ is positive.
We relate the rate of energy dissipation G(x, x) to a dissipative force f(x, x)
f(x, x)x = G(x, x). (2.73)
\ G(x,x) J{x,x) = X (2.74)
The dynamics including this dissipative force is
mx + = fix, x).
Substituting G(x,x)/x for f(x,x) gives
dV G(x,x)
mx + = v 1
dx x
- dV n / fx\2\
Choosing the simple harmonic oscillator potential
V(x) = ^kx2,
we have after some rearrangement:
mx 2m/y ^1 ^ ^ x + kx = 0

Dividing through by m and using u2 = k/m,
x 27(1 ())2x + lu2 = 0.
If we use a length scale a and time scale ^ then in dimensionless form,
the equation becomes
where e = ^ = x and = x. This is the Van der Pol equation.
While a mechanical realization of this type of energy exchange may not be
easy to implement, it turns out that this equation does describe certain types of
electronic oscillators. Here we examine the Van der Pol equation as an abstract
dynamical system. Moreover, we will find it useful to examine behavior when
physical ingredient of negative dissipation, that is, of being able to feed energy
into the system. We will find that a new type of dynamics becomes possible,
namely self-sustained oscillation.
2.3.1 Hopf Bifurcation and the Limit Cycle
Let us consider the dynamics near the fixed point at (x, v) = (0,0). What
happens to tiny disturbances to this fixed point? For short times, such distur-
bances may grow or decay according to dynamics that neglects all terms higher
than linear order in the disturbances. We can re-write the Van der Pol equation
(2.81) as two first order difrential equation:
e < 0 as well as for e > 0. A key question concerns what happens due to the new

Thus we examine a linearized version of the above equation
= v
= ev x
The term ex2 has been excluded in this approximation. The linearized system
may be written in matrix form
dx dt 0 1 X
dv \ dt J l -l£ J i J
Suppose we found eigenvectors of the above matrix such that:
f n A / \ / \
z = A x
v-le/ \v) {vj
where A represent the two eigenvalues of the matrix. The significance of the
() notation will be seen below. Look at the special case where the initial
conditions conform to one of these eigenvectors.
a;(0)^ f \ x
v^()y lV)
In fact, a solution to system satisfying this initial condition is
We see this by first noting that
= ext
= AeAt
/ \
( \

But also
f01\/ A
\v(t) J
= e
/ n \
0 1 x
= eAtA^
/ \
Consequently we see that
0 1
1 £
\ ( _.A
( ,nA ( \ f \
X(0) = eA x X
1^(0)) {vj {v)
Thus we have a solution matching the prescribed initial conditions.
Now an arbitrary initial condition can be expressed in terms of the two
eigenvectors (except in special cases where they are not linearly independent):
( ,-A ( \ (
*(0) 1 x+ ,
= a+ + a_
^(0) J v+J V

The corresponding solution is
( x(t)^ = a+ex+t ( \ x+ + one*-4 ( \ X-
yV(t)j {v+j v->
So it remains only to find the eigenvalues and eigenvectors of the Van der Pol
equation. The eigenvalues satisfy
^0 A 1 ^
T e Aj
which yields the characteristic equation
= 0,
A?. eA + 1 0

The eigenvalues are roots of this quadratic equation:
A+ =
e \/e2 ~ 4
corresponding eigenvectors are in fact
values of e satisfying
/ \ /
-2 < e < 2
Now for the range of
we will have complex eigenvalues
A+ = i
. \/4 e2
The time dependence is accordingly
M M (i)
X = e£t/2 (a+)e^)t + (a-)e~iw{£)t
(a+J w
where there is an oscillatory component of frequency
w(£) = V 1(l)'
1 for small e.
The solution decays exponentially for e < 0 and grows exponentially for e > 0
(figure (2.17)). We conclude that e = 0 marks a special transition point, or
bifurcation, where the disturbances to the fixed point can change from a damped
to an exponentially growing solution. As e increases through 0, the fixed point
is unstable and is replaced at least for short times by growing oscillatory
terms. As we will discuss below, the amplitude of the oscillatory mode eventually
saturates as it becomes too large to ignore the nonlinearity of the original Van

der Pol equation. In fact, the fixed point solution is replaced by a periodic
motion of specific amplitude which we call a limit cycle. The system is said to
have undergone a Hopf bifurcation when this transition from a time-independent
solution to limit cycle motion occurs.
Figure 2.17: Time series (top) and phase trajectory (bottom), critical control
parameter, tc = 0: (left) e = 0.2 < 0, stable spiral and (right) e 0.2 > 0,
unstable spiral.
It can be shown that an approximate solution to the Van der Pol equation
for e >C 1 with averaging method is given by [37],
- -= cos t.
Vl + 3e-£t
The approximate solution (equation (2.103)) and numerical solution by Runge-
Kutta method for Van der Pol oscillator are compared in figure (2.18).

Figure 2.18: The approximate solution by averaging method (equation (2.103))
and numerical solution by Runge-Kutta method for Van der Pol oscillator: (e =
0.2 and x(0) = 1, ^'(O) = w(0) = 0)-
42 Limit Cycle in Van der Pol Oscillator
Before proceeding further with the Van der Pol equation we make some
general remarks about limit cycles and Hopf bifurcations. Cyclic or periodic
response is possible in a system of two or higher dimensions. The trajectory of
the system in this case in the phase space is a closed loop. If the trajectories in
the neighborhood of the closed loop attract toward or repel from the loop then
the loop is called a limit cycle. A limit cycle is an isolated closed trajectory.
Isolated means that the trajectories in the neighborhood of limit cycle are not
closed: they either spiral toward or away from the limit cycle [37]. The phase
trajectory of a Van der Pol oscillator in figure (2.19) shows a trajectory from
outside and one from inside of the limit cycle spiraling toward the limit cycle.
If neighboring trajectories approach the limit cycle it is called an attracting
or stable limit cycle. If neighboring trajectories move away from the limit cycle,
it is known as repelling or unstable limit cycle. The stable limit cycle shows a
self oscillation, meaning that system can oscillate without an external periodic
driving force.
Let us restrict the problem to two dimensions. In two dimensions the system
linearized about a fixed point (the origin) has two eigenvalues. The eigenvalues
for a limit cycle are complex. If the real parts of eigenvalues are positive the fixed
point becomes unstable and is replaced with a limit cycle. But if the eigenvalues
have negative real parts then the fixed point is stable. If the linearized system
has pure imaginary (no real part) eigenvalues then the fixed point is called a

Figure 2.19: Limit cycle of Van der Pol oscillator: (e = 0.2 and top: x(0)
ljO^O) = w(0) = 0, bottom: x(O) = 1, ^'(O) = i;(0) = 2).

center and the trajectory in a Poincare section is a closed orbit for any initial
condition and this case is not called a limit cycle. One example of the center
fixed point is the simple harmonic oscillator which has a sinusoidal time series
and its trajectory in phase space is an elliptical orbit.
The following second order equation is a generalization of the Van der Pol
equation and is known as the Lienard equation:
§+/()!+*() = 0, (2.104)
where f(x) is even function and g{x) is odd function. It can describe a system
of mass and spring with nonlinear damping coefficient f(x) and nonlinear spring
force g(x). This equation has a unique, stable limit cycle under appropriate
hypotheses on f(x) and g(x) [37]. We focus on the Van der Pol oscillator as one
kind of Lienard system since it has a wide range of application. By comparing
the Van der Pol equation with the Lienard equation we obtain f(x) and g(x)
for Van der Pol equation
f(x) = e(x2 1) (2.105)
g(x) = x. (2.106) Hopf Bifurcation in Van der Pol Oscillator
As we have seen the limit cycle occurs in a system that is at least two
dimensional. The periodic solutions of equations linearized around a fixed point
have eigenvalues that are complex. A bifurcation occurs when the sign of the
real part of the eigenvalues changes. The appearance and disappearance of a
limit cycle is called a Hopf bifurcation [20], [37].

Suppose the eigenvalues of a system depend on a control parameter ji. As-
sume a fixed point of the system for control parameter lower than a critical value
(j,c is a spiral attractor (stable spiral) which has negative real eigenvalues. If for
control parameter larger than pc the stable fixed point becomes unstable and a
limit cycle appears then a supercritical Hopf bifurcation occurs.
Figure 2.20: Hopf bifurcation in Van der Pol equation (time series (top) and
phase trajectory (bottom), critical control parameter, ec = 0: (Left) e = 0.2 <
0, stable spiral and (right) e = 0.2 > 0, stable limit cycle.
A Hopf bifurcation occurs at e = 0 because the real part of the eigenvalues
are negative for e < 0 and positive for e > 0. A Hopf bifurcation for the Van der
Pol oscillator is shown in figure (2.20). For e < 0 the fixed point (x = 0, v = 0) is
stable and any trajectories converge to the fixed point and they are stable spirals.
For e > 0 the fixed point is unstable and the system has a limit cycle which is
stable. Any trajectories from inside or outside or the limit cycle converge to the

limit cycle. A Hopf bifurcation occurs at ec = 0 which is known as its critical
2.3.2 Relaxation Oscillator
The relaxation oscillation is a kind of limit cycle with fast and slow portion
orbit on phase space. Let us write the Van der Pol equation again
+ (2.107)
The above Van der Pol equation for fj, ~^> 1 is a strongly nonlinear system which
show a relaxation dynamics (figures (2.21)-(2.23)) [37].
The period of relaxation oscillation can be determined approximately by
analytical methods. The period of Van der Pol for yx 1 is given approximately
by [37],
T ji(3 2 In 2). (2.108)
The time series of relaxation oscillators for different values of p are given in
figure (2.21). We observe that for large /x, equation (2.108) gives a better ap-
proximation. The phase trajectories of relaxation for different p and the initial
conditions inside the limit cycle and outside limit cycle are given in figures (2.22)
and (2.23).

Figure 2.21: Time series for relaxation oscillator: (a) /x = l,:r(0) = ^^'(O)
u(O) = 0; (b) fj, = 5,z(0) = 1,^(O) = v(0) = 0 and (c) /x = 10,:r(0) = 1, ^(O)
v(O) = 0.

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
Figure 2.22: Limit cycle of relaxation oscillator (initial condition inside the
limit cycle), (a) /a = l,:r(0) = 1, ^(O) = d(0) = 0, (b) (a = 5,a;(0) = l^^O) =
u(0) = 0 and (c) /i = 10, x(0) = 1, ^'(O) = u(0) = 0.

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
Figure 2.23: Limit cycle of relaxation oscillator (initial condition outside the
limit cycle), (a) jj, = l,a:(0) = ljZ^O) = f(0) = 2, (b) fj, 5,o;(0) = 1, rc'(O) =
t>(0) = 3 and (c) = 10, z(0) = 1, ^'(O) = ^(0) = 5.

3. Electrical Oscillators
In this chapter we will consider some electrical oscillators. First we review
some basic concepts on charging and discharging capacitive loads. This will be
important for understanding certain types of relaxation oscillator circuits which
later will serve a basis for integrate and fire models of neuron behavior. Then
we will see the neon and 555 circuits as relaxation oscillators and at the end the
Van der Pol Oscillator is discussed.
3.1 Resistor Capacitor (RC) Circuits
Suppose a circuit consists of a capacitor with capacitance C, a resistor of
resistance R, a battery with emf E and a switch as in figure (3.1). In this circuit
the resistor and capacitor are in series.
Figure 3.1: Resistor capacitor (RC) circuit.

3.1.1 Charging the Capacitor
As the switch is closed, at position a in figure (3.1), the capacitor starts to
be charged. This means that current begins to flow across the capacitor as a
result of an electric field across the resistor. The current increases the charge q
on the plates of capacitor and the potential difference across the capacitor is
Ve = q/C. (3.1)
As the voltage reaches the potential difference of the battery (E) there is no
voltage drop across the resistor and the current becomes zero. The capacitor is
fully charged and it is in equilibrium. By applying the loop rule in the circuit
we have
E-iR-^ = 0.
Current is the rate of change of charge:
Substituting (3.3) into equation (3.2), we obtain an equation for charging the
Rd RTt + c~K
Using equation (3.1) this can be presented in terms of the potential difference
of the capacitor Vc
CR^f + VC = E.

The amount of charge at the beginning is zero. Using this fact as the initial
condition, the solution of equation (3.4) is
q = CE( 1 e~t/RC) . (3.6)
The current during the charging is
4 Hi) ^ (3'7)
Figure (3.2) shows the charge and current versus time for charging the RC
Figure 3.2: Charge and current versus time for charging resistor capacitor
(RC) circuit.
The potential (voltage) across the capacitor is
VC=^ = E{ 1 e~t/RC) . (3.8)

The capacitive time constant is
tc = RC. (3.9)
Charge on the capacitor reaches 63% of the final value (CE) after a time t = rc.
3.1.2 Discharging the Capacitor
Assume the capacitor is fully charged to a potential of Vo- Now if the switch
is opened, at position b in figure (3.1), then the battery is not in the circuit.
Therefore the equation for discharging is
Figure 3.3: Charge and current versus time for discharging resistor capacitor
(RC) circuit.
Rft + h = 0- (3-10)
Equation (3.10) in terms of the capacitors voltage is
RC- + FC = 0. (3.11)

The solution of equation (3.10) is
Q = ?oe
where qo = CVq.
The potential across the capacitor is
K = Ve-'sc.
The current is
The charge and current are plotted versus time for discharging RC circuit in
figure (3.3).
3.2 Relaxation Oscillator
This section considers the neon and 555 oscillators as two examples of re-
laxation oscillators.
3.2.1 Neon Bulb
A neon tube includes two electrodes separated by a gap of a few millimeters
in an atmosphere of neon. The neon tube is a nonlinear device whiclndoes not
obey Ohms law. A circuit with a neon bulb is shown in figure (3.4). When the
voltage across a neon tube exceeds a certain value the neon gas breaks down
into an ionized plasma that conducts electric current and gives off light. An
electron free in the gas acquires enough energy between collisions with atoms
to ionize an atom on the next collision. This gives rise to a cascade of released
electrons that in turn collide to produce more ions and electrons. This is known

as a charge avalanche. Collisions from the cascade of electrons excite atoms to
produce light as the atoms return to their ground state.
The R and C combine to build the voltage across the tube with time in
a manner described in the previous section. When the breakdown threshold is
reached the tube begins conducting and discharges the capacitor. Eventually
at a lower threshold the cascade of electrons shuts off. The cycle repeats itself,
producing a free running oscillator.
Figure 3.4: Neon tube circuit.
3.2.2 Neon Tube Relaxation Circuit
The neon bulb circuit, figure (3.4), has a resistor with the large resistance
R, a switch, a capacitor with the capacitance of C and a neon bulb. Fig-
ure (3.5) shows the nonlinear voltage-current characteristic of the bulb and the
corresponding voltage versus time is shown in figure (3.6). When the switch
is closed, almost no current flows through the bulb and the capacitor slowly

charges (from point 0 to A) until the potential across the capacitor reaches its
maximum at the firing voltage (V/) at point A and the tube turns on. Auto-
matically switching from charging through resistor R (at point A) to discharging
through neon bulb (at point B) happens instantaneously. When the bulb turns
on, its resistance drops and the capacitor rapidly discharges through the neon
tube because its resistance is low. As the capacitor discharges (from B to C),
the potential across the capacitor decreases until the extinction voltage (Ve)
is reached and the tube turns off. As is shown in figures (3.5) and (3.6) from
point C (end of discharging) to poirit D (beginning of recharging) occurs instan-
taneously. When the bulb turns off, its resistance increases and the capacitor
begins to recharge and the limit cycle is in operation. Note that if the value of
Figure 3.5: Voltage V versus current I for neon bulb.
the charging resistor R is less than some critical value then it is possible that
the bulb stays on [12]. Figure (3.6) is the time series and since I = % =

Figure 3.6: Voltage V versus time t for neon bulb.
figure (3.5) is a phase trajectory which is a limit cycle for a relaxation oscillator.
From figure (3.6) we can find the period of relaxation oscillator. The voltage
across the charging capacitor is
VC = E( 1 e-'/).
Let t = te the time correspond to extinction voltage. Then
E Ve .
= e
For the firing voltage we have
E____Yl .-tt nno)
Therefore we have
te RC In

tf = RC In
The period of oscillation is
t = tf te = RC In
Let Ve = (1/3)E and Vf = (2/3)E then
t = tf te = RC ln(2)
For R = 100.KT2 and C = 3.5/j,F we get a period of r = 0.24s and a frequency
of 4.1 HZ.
3.2.3 The 555 Relaxation Oscillator
The 555 chips produce a stable oscillation [22].. Figure (3.7) shows a 555
circuit. The capacitor starts to be discharged when power is applied. The 555
has an internal switch that closes when the voltage exceeds \VCC and opens when
the voltage drops to |Vcc, where Vcc is the power supply voltage [15]. The 555
is triggered and therefore the output goes high and the discharging is stopped.
The charging through Ra and Rb occurs again. When the voltage across C
reaches §Vcc, the threshold input is triggered again and as a result the output
goes low and then the discharging through Rb occurs until reaches This
makes an oscillator between |Vcc and |Vcc- The time of charging is
Tc = ln{2)(RA + RB)C (3.22)
= 0.693 (Ra + RB)C. (3.23)

Figure 3.7: The 555 oscillator circuit.
The discharging time is
Td = H2)RBC (3.24)
= 0.693RbC. (3.25)
Then the period of the oscillator is
T = TC + Td (3.26)
= ln(2){RA + 2RB)C (3.27)
= 0.693(i?A + 2 Rb)C. (3.28)
Figure (3.8) shows the responses for trigger (V2) and out (V3).
3.3 Oscillation by Gain and Positive Feedback
The relaxation oscillators discussed in the previous section depend on the
coupling of charge integrators to a threshold mechanism for switching the circuit

Figure 3.8: Voltage V versus time, t for 555 oscillator.

into and out of a discharge state. Oscillation occurs by repeated integration and
discharge organized by the nonlinear switching. Another class of oscillators uses
amplification (gain) and positive feedback to sustain oscillation. This is shown
schematically in figure (3.9).
Figure 3.9: General configuration of a system with gain and feedback.
We first assume linear behavior of all elements of this system. Suppose a
sinusoidal signal of frequency u is applied at the terminal marked V' given by
V' = Via + H(u)Vout.
Here H{uS) is a frequency dependent weighting factor. Now the role of the
amplifier is to multiply V' by a large gain factor G{oj) that may also be frequency
yout = G(co)V'. (3.30)
We can eliminate V' from equations (3.29) and (3.30) to get
Vout = G{oj) (t>in + H(u)V0ut) . (3.31)
Solving for Voul, we get
1 G(u>)H(uj)

The expression multiplying Vin is called the loop gain of the system:
Gloop(w) = l-G(u)H(ujy
What happens if we can construct a system so that at a particular frequency
u>o we have G(ujo)H(u>q) = 1? As we see
lim Gioapito) > oo. (3.34)
Heuristically we argue that Vout can grow to a finite level even in the absence
of Vin. That is, the system shown in figure (3.10) is capable of self-sustained
oscillation at frequency lo0. In practice, the feedback element has a nonlinear
Figure 3.10: A system capable of self-sustained oscillation at frequency too
where G(u>0)H(u>0) = 1.
response that depends on the output amplitude
H = H(oj,V0ut). (3.35)
If we arrange for GH to fall below unity when Vout exceeds some saturation
level, Vsat, the system will set up a limit cycle oscillation at an amplitude of
= Vsat. (3.36)

This heuristic argument, common to many electronics textbooks, bypasses a
detailed analysis of the dynamical behavior of the system by assuming a priori
- that the signals are sinusoidal. In actual behavior, practical oscillator circuits
can show the full richness of nonlinear systems including non-sinusoidal limit
cycles and when coupled together to form a system of a sufficient number of
degrees of freedom chaos.
3.3.1 The Wien Bridge Oscillator
To proceed further, we consider a specific implementation of figure (3.10)
known as the Wien bridge oscillator. The circuit is shown in figure (3.11).
Figure 3.11: Wien Bridge oscillator. For comparison to figure (3.10) we note
that V' = V+-V-.
A detailed analysis of this circuit may be found in [12]. The salient features
are the following:

1. The network of resistors and capacitors forms a frequency selective feed-
back. The condition for net positive feedback and sustained oscillation
occurs when
i?2 = 2i?io (3.38)
2. The resistor Ri is made of a material whose resistance increases with the
voltage drop across the resistor. Here Rio is the resistance of R\ when V'
goes to zero. Typically this is the tungsten filament of an incandescent
light bulb in series with an ordinary resistor. Self-heating of the filament
causes its resistance to increase.
Ri = Rw +:c*V/2. (3.39)
This provides the nonlinear mechanism for saturating the output voltage
A dynamical analysis of this circuit, [12], shows that the output voltage is in
fact governed by the Van der Pol equation when resistor Ri has the dependence
given in equation
d2V0 \Rz-2R10-2aV?
dt2 [ R10 + R2
where an assumption has been made that
Rio + R2 olV2 (3.41)
This means that the change in resistance Ri due to self-heating is small compared
to the sum of i?io and R^-
+ uJqVo 0

A key feature of this equation is that the onset of oscillations is governed
by the bifurcation parameter.
i?2 2Rio
e = 3-----------
Rio + i?2
By literally turning the knob that adjusts R2 one obtains oscillation when
i.e. when i?2 > ^R\o- The amplitude of oscillations is very sensitive to the
AR2 by which i?2 exceeds 2A10.
£ > 0

4. Oscillators in Biology
4.1 Excitable Neurons
The nervous system is a network that performs communication between
cells. Nervous tissue has two kinds of cells, which are known as neurons and
glial cells. Neurons are for transmitting information. Glial cells surround the
neurons and protect them. Neurons have four parts: dendrites, a cell body
(soma), an axon, and synaptic terminals. Dendrites have a large number of
branches (on the order of 100,000). These branches receive input signals from
other cells and transmit them to the soma. There is a mechanism for summing
up input signals, which can individually either act to excite or inhibit response
by the neuron. The cell body (soma) contains the nucleus, which holds the
genetic information. The axon transmits messages away from the cell body.
Synaptic terminals are at the end of axon. They release chemical messengers
upon stimulus by electrical signals from the axon. Myelin is formed from plasma
membranes of Schwann cells to create an insulating cover for the axons. The
spaces between Schwann cells are called nodes of Ranvier and are important for
ion exchange necessary for propagating signals down the axon.
There are three types of neurons: afferent neurons (sensory), efferent neu-
rons (motor), and interneurons. Afferent neurons carry messages from the or-
gans or sensory cells into the central nervous system. Efferent neurons transmit
messages from the central nervous system to muscles or glands. Interneurons
connect neurons to each other and are located in the central nervous system.

In a neurons resting condition, its intracellular medium is negatively-
charged to a potential of -70 mV relative to the extracellular medium. At the
resting potential no travelling electrical signal is produced. The resting poten-
tial depends on the concentration of ions on both sides of the membrane and
permeability of the membrane for different ions. In the rest state the concen-
tration of Na+ in the extracellular medium is higher, while the concentration of
K+ inside the cell is higher. These concentrations are regulated by membrane
channels (gates) that pass specific types of ions and by mechanisms for pump-
ing ions against established concentration gradients. Most electrical activity is
generated by the Na+ and K+ gates. The permeability of the membrane for K+
is normally larger than for Na+ [14], [41].
An action potential is a reversal of the potential across the membrane that
occurs rapidly and temporarily in response to a stimulus. During the action
potential first sodium and then potassium channels open. Sodium ions move
into the cell via sodium gates and generate the initial change in polarization
that produces an action potential. As a result the charge inside becomes more
positive than outside. Then potassium ions move to the outside to reset the
potential. At the last stage sodium ions are returned back to the outside and
potassium ions are pumped into the cell to reset the ionic concentration on both
sides of the cell membrane and return to the original (resting) condition. Other
ions, such as a Ca2+, can play a role in modulating this process.
The action potential occurs as a threshold phenomenon. A stimulus causes
some change in membrane potential. The stimulus must be strong enough for
the membrane potential to reach a threshold and then the action potential will

be produced. A stronger stimulus does not increase the amplitude of the ac-
tion potential. During the action potential another action potential cannot be
generated: this inhibition during and immediately after triggering of the action
potential is called the refractory period. The action potential in one part of a
membrane acts as a stimulus to cause another action potential in an adjacent
part of the membrane by a local current. As a result of the sequence of trig-
gering of action potentials in adjacent parts of the cell membrane the action
potential propagates down the axon. Finally the signal is sent to other neurons
via synapses [41]. >
4.1.1 Hodgkin-Huxley Model
A. L. Hodgkin and A. F. Huxley built a mathematical model for neuron
excitation [21], [3], [27] based on an ingenious series of experiments using the
giant axon of a squid. They won the Nobel Prize in physiology/medicine for
their work.
Hodgkin-Huxley modelled the electrical current across a patch of nerve mem-
brane by an electrical circuit shown in figure (4.1). The input current can charge
the membrane capacitance or can move into the cell through sodium, potassium
and leakage gates. The leakage current stands for other channels that are not
described and consists mostly of Cl- ions.
Looking at figure (4.1) Kirchhoffs law implies
i(t) = ic(t)+Y,h(ty (4-i)

Inside Cell
I(t) \
T e, ~Te T e
I L | K | Na
Outside Cell
Figure 4.1: Schematic diagram of Hodgkin-Huxley model.
Here I(t) is input current, Icif) is capacitors charging current and h(t) is
the sum of the ionic currents as follows:
= I Naif) + Inf) + (4.2)
where lNa(t)-, h<[t) and /l(£) are currents through the sodium, potassium and
leakage channels respectively. Let C be the membrane capacitance and V the
voltage across the membrane. Then we have
Ic(t) = C-^. (4.3)
c^ = ~Y,m + m- (4-4)
In the Hodgkin-Huxley model the current through each channel is given by
iNa = gNa^hiy Eno),
h = gc,(V-BL),

where gxa, Qk and are the maximum conductivities of sodium, potassium and
leakage channels respectively. It is assumed that the conductivity of leakage
gate is constant at value gl- However, the conductivity of the sodium and
potassium channels depend on both voltage and time through the parameters
m, h and n. The variables m, h and n are called gating variables and describe
the percentage of a given type of channel that is open. They are functions of V
and t (see below). The parameters E^a, Ex and El are reversal potentials that
represent the thermodynamic equilibrium ionic concentrations through Nernst-
type relations. For example,

where [Na+]i is the Na+ concentration inside the cell, [Na+}0 is the Na+ con-
centration outside the cell, kg is Boltzmans constant, T is the absolute temper-
ature, Z is the number of ionizations (+1 for sodium) and e is the magnitude
of the electron charge. Finally the ionic current can be expressed as:
A(t) = ll!:Ujnh(V EN) + gKn4(V EK) + gL(V EL).
In this model the empirical parameters for conductance and reversal poten-
tial are given by [3]
ENa = 55mV, Ek = 72 mV, EL = 49.387 mV, C = 1 [iF/cm2, (4.10)
9Na - 120 mS/cm2, gx = 36 mS/cm2, gL = 0.3 mS/cm2, (4.11)
where mS/cm2 is millisiemens per square centimeter (conductivity units). Sub-
stituting equation (4.9) into (4.4) we get the final form of the Kirchhoff law for

this model. The final Kirchhoff equation and the gate variable equations are a
system of four coupled nonlinear differential equation as follows:
C = -gNam3h(V ENa) gKnA
= am(V)(l m) (3m(V)m
= an(V)(l n) (3n(V)n
^ = ak(V)(l -h.)-fa(V)h.
(V -EK) + gL(V -Et)
The equations for m, n and h are first-order kinetic equations for describing
opening/closing gate of the form [3]
Close ^ Open,
where a and /5 are corresponding rate constants. The empirical functions a and
/3 are given by
<*n{V) = -
3.5 + 0.11/
1 e(3.5+0.1V)
0.5 + 0.011/
e-(5+0.1V) 1
ah(V) = 0.07e(6O+y)/2
Pm(V) = 4e-^0+v^18
/3n(V) = 0.125e-(6O+v)/8 (4.21)
Ph{V) e_(3+o.iv) _|_ ^' (4.22)

Note that these equations and their parameters are different from the original
Hodgkin-Huxley model since they are based on redefinition of the membrane
potential according to current usage [3].
The Hodgkin-Huxley model accurately captures many features of the be-
havior of a neuron, especially the squid axon. For detailed comparisons, see
Weiss, chapter (4), [42]. However, it is a complex model and more intuition
about its behavior can be obtained using a simplified model that still captures
the qualitative features. This is discussed in the next section.
4.1.2 Fitzhugh-Nagumo Model'
The Fitzhugh-Nagumo model describes the action potential of a neuron
[13], [23], [3]. This model is a simplification of the Hodgkin-Huxley model. The
system of first order nonlinear differential equations that describes this model is
given by
This equation reduces to a kind of Van der Pol equation for a = b = z = 0. Here
2 is the input current and it corresponds to /, x represents both the membrane
potential V and the gate variable m, and y plays the role of the two gate variables
n and h in Hodgkin-Huxley equation. These combinations of V and m into a
single variable x and of n and h into a single variable y are justified by observing
that in each pair, the original variables tend to track one another over similar
time scales and might be considered to be transformations of a single variable.
{x a + by)

The system parameters, a, b and c are constants and satisfy
1 < a < 1, 0 < 6 < 1, b < c2. (4.25)
We can used the following data which satisfy these conditions.
a = 0.7, b 0.8, c 3.
The fixed point of the system corresponds to the rest state and is given by
"3 N /dx x
= 0
y + x- + zj=0
x-a + by = 0 ( = 0 ] ,
This is a cubic equation for x and it has only one solution for the above param-
eters. Thus there is only one fixed point for the system.
Another point of view is obtained by using the above equations to define
nullclines of the system
x3 ( dx
y ~z x z ( nullcline where
3 y dt
The first one is a cubic equation and the second one is a linear equation. The
fixed point is the intersection of these two equations which can be found easily
by graphing them on the (x, y) phase plane. For the given values of parameters
and z = 0 the fixed point is (Xf = 1.20, yj = 0.625), (see figure (4.2)). It
should be noted that this geometrical configuration appears to be a common
feature for excitable media, of which neurons are one example.
Figure (4.3) shows the response of the system for z = 0. The dynamics has
one stable fixed point. The fixed point for z = 0.4 is unstable and the system
1 a
(nullcline where ^ = 0
V dt

Figure 4.2: x and y Nullclines for Fitzhuqh-Naqumo Model: fa = 0.7,6
0.8, c = 3 and z = 0).

Figure 4.3: Time series and phase trajectory for Fitzhugh-Nagumo model:
(a = 0.7, b 0.8, c = 3 and z = 0).

Figure 4.4: Time series and phase trajectory for Fitzhugh-Nagumo model:
(a = 0.7, b = 0.8, c = 3 and z = 0A).

Figure 4.5: Time series and phase trajectory for Fitzhugh-Nagumo model:
(a = 0.7,b = 0.8,c = 3,z = 0.336) and initial condition (0.25, 0.1).

has one stable limit cycle (figure (4.4)). The fixed point for the case of z = 0.4
is (0.9066,-0.2582). The fixed point is stable for 0.3456 < z. Stable fixed
points correspond to a non-firing state. For 0.336 < 2 all trajectories with
different initial states converge to the fixed point. Figure (4.5) shows that for
^ = 0.336 the trajectory is attracted to the fixed point at (0.96348, 0.32935).
The interesting case is the region of 0.3456 < z < 0.336. In this interval
there are two limit cycles. The interior limit cycle is unstable so any trajectory
that starts from inside this cycle goes away from this limit cycle and converges
to the fixed point. The external limit cycle is stable and thus any trajectory
between the two limit cycles moves away from the internal one and is attracted
to the external limit cycle. A trajectory that starts from outside the exterior
limit cycle approaches this limit cycle. For example the system for 0.3456 <
z = 0.34 < 0.336 has a fixed point at (xf 0.9601, yf = 0.3251). Figures
(4.6)-(4.10) show the phase trajectory for different initial conditions and the
phase portrait is shown in figure (4.11). These figures show the behavior of
system for different trajectories relative to the two limit cycles and the fixed
In summary for z < 0.3456 the system has one stable limit cycle and the
fixed point is unstable. Between 0.3456 < z < 0.336 the system has two
limit cycles, where the interior one is unstable and the exterior one is stable.
For 0.336 < z dynamics has no limit cycle and the fixed point is stable. At
z 0.3456 a Hopf bifurcation occurs due to the disappearance of a limit cycle
and at z 0.336 another Hopf bifurcation occurs through the disappearance

Figure 4.6: Time series and phase trajectory for Fitzhugh-N agumo model:
(a = 0.7, b = 0.8, c = 3, z = 0.34,) and initial condition (0.6,1.4).

Figure 4.7: Time series and phase trajectory for Fitzhugh-Nagumo model:
(a = 0.7, b = 0.8, c = 3,z = 0.34) and initial condition (0.6, 0.2).