1. health and fitness
2. disease

How to Build Spiking CPG Models Using Python
M. Anthony Lewis
Abstract Central pattern generators underly walking, swimming, flying, and crawling in animals. CPGs have been found in both vertebrate and invertebrates. Understanding the function and role of CPGs is critical to understanding movement in animals, and therefore humans. In this brief tutorial, I will show how to build simple
CPGs models using Spiking neurons. I will make the examples concrete by showing how to code these models in Python, a free, widely available and simple to use
computer language. It is hoped that this tutorial is self-contained and is accessible
to researchers and students with a minimum of training in Calculs.
1 Introduction
Cells are the basic building block of multi-cellular organisms. Special cells exist
that are electrically excitable. These cells include muscles cells and in neural cells
(Neurons).
The cell membrane (Fig. 1) can be thought of as a bag containing the cellular
machinery that allows an imbalance in ion concentration to occur inside versus outside the cells. Computationally we do not care about the cells DNA, mitochondria,
cell nuclei etc. We do care about this ion imbalance. The cell membrane itself has a
very high resistance. On either side of the cell membrane is a conductive solution,
thus the cell membrane forms the dielectric of a capacitor. Ion channels penetrate
the cell bilipid (composed of two layers of fat molecules) membrane layer and are
highly selective to specific ion species.
M. Anthony Lewis
Department of Electrical and Computer Engineering, University of Arizona, Tucson, AZ, 85721,
e-mail: malewis@ece.arizona.edu
1
2
M. Anthony Lewis
Fig. 1 The Cell membrane is composed of a bilipid insulating layer and isolate conductive electrolytes inside and outside the cell. Certain ion species, such as potassiam, calcium and sodium are
particularly important to electrical activity in neurons. The cell membrane contains ion channels
that are highly selective to specific ions and can be turned on and off by changes in voltage or
by a special protein (ligand). The cell membrane also contains ion pumps that maintain a specific
ion imbalance between the inside and outside of the cells. Opening ion channels take advantage of
this imbalance to create current flows and hence alter membrane potential. Highly dynamic events
called action potential encode voltage events and propagate this information to thousands of other
cells.
1.1 Information Flows into dendrites and out of axons
Neural cells have two types of processes that are attached to the cell body. The axon
is a long process that can transmit action potentials, discussed below, from the cell
body to other neurons throughout the brain. These axons synapse on the dendrites of
other neurons. These dendrites collect information from thousands of neurons. Each
synapse contributes a bit of current to the cell and results in the alteration of the cell
membrane voltage. Each synapse may have a different coupling strength. A given
cell can either make it more likely or less likely that a cell will generate an action
potential, fire, depending on if it is an excitatory synapse or an inhibitory synapse.
1.2 The Neuron Cell is a Capacitor with a Decision Making
Capability
The membrane/electrolyte system forms a capacitor. The electrical model of an ideal
capacitor is:
dV
C
=i
(1)
dt
How to Build Spiking CPG Models Using Python
3
where C is the capacitance, i is a current flow, and V is the voltage potential measured relative to the inside of the neuron.
Ions flow is dirven by diffusion and voltage driving force. Ion pumps in the cell
membrane maintain a concentration gradient for various ions.
Ion Channels allow charged ions to flow across the cell membrane in a highly
selective way. Given a ion gradient between the inside and outside cell, a diffusion
gradient is set up that causes a net flow of ions from the higher to the lower gradient.
As each ion particle is charged, this causes a current flow. If a voltage potential is
maintained across the membrane, positive ions will be attracted to the negative side
of the cell membrane. This leads to the build up of a concentration gradient.
The Nernst potential specifies the relationship between concentration gradients
and the voltage potential across the membrane:
E=
RT [outside]
ln
zF
[inside]
(2)
where E is the potential, T is the absolute temperature, and R the gas constant,
F is the Faraday constant. For real cells the Ek = [−70.. − 90mv], ENa = [50mv],
ECa2+ = [150mv].
1.3 Neural Models capture the basic dynamics of the cell body and
throw away some details
The difference in neural models comes down the to the equations specifying I in
eqn (1). The trade off is between the complexity of the dynamics, the ability to analyze the systems and the ability to simulate a given model efficiently. The most
well known model is the Hodgkin-Huxley model. Hodkin-Huxley empirically determined their equations [?, ?]:
im = gL (V − EL ) + gk n4 (V − EK ) + gNa m4 h(V − ENa ) + Is
(3)
h,n,m evolve according to complex equations, Is are synaptic currents coming from
other neurons. The equation can be understood in the following way: the values in
parentheses (V − Ex ) are driving forces that will tend to drive the potential of the
neuron to the Nerst reversal potential Ex . In front of each of these driving forces
are weighting factors. These weighting factors determine which driving force will
dominate the equation, and hence the equilibrium value for the cell membrane. If
only gl is active, the cell membrane will relax toward V − > El . When the cell begins
to fire, after the cell membrane has increased in potential, the Na part of the equation
dominates and the membrane potential shoots up toward ENa . Next the potassium
term turns on and resets the neuron.
In general the integration of this equation is difficult, and the analysis is extraordinarily difficult, so researchers have turned to simplified models which capture some
of the features of this system. Here we focus on computationally efficient models.
4
M. Anthony Lewis
1.3.1 Leaky Integrator
One of the simplest models is the so called leaky integrator. This model avoids the
complex dynamics of firing altogether:
cm
dVi
= −gl Vi + Ii
dt
ui = f (vi )
(4)
(5)
N
Ii = ∑ wi j u j
(6)
i=1
where the V is the membrane potential, u is the average firing rate of the neuron
over a small time window. f (v) is a function that transforms membrane voltage into
firing rate.
i
In steady state dV
dt = 0. This implies that gl Vi = Ii . We can often assume that
gl = 1. This equation is an integrator because absent the −Vi term, it is a perfect
integrator. This equation is leaky because absent the input Ii , the voltage Vi “leaks”
back to zero.
This model does not capture the essential spiking characteristic of neurons. All
precise timing information is eliminated. The next step up in complexity is an integrate and fire model.
1.3.2 Integrate and Fire
dVi
= −gL (Vi − EL ) + Ii
dt
i f (V > Vthres )− > V = Vreset
cm
(7)
(8)
N
Ii = ∑ wi j s(u j )
(9)
i=1
In the integrate and fire model, as the membrane approaches a fixed threshold, the
cell ’fires’ and generates a spike. It then resets the membrane voltage. This model is
not sophisticated enough to build interesting models of networks controlling motors
systems. We must introduce spike frequency adaptation:
dVi
= −gL (Vi − EL ) + rm gsra (V − Ek ) + Ii
dt
dgsra
τsra
= −gsra
dt
i f (V > Vthres )− > V = Vreset , gsra − > gsra + ∆ gsra
cm
(10)
(11)
(12)
N
Ii = ∑ wi j s(u j )
i=1
(13)
How to Build Spiking CPG Models Using Python
5
Here we have an adaptation terms. The variable gsra functions much like a spike averager, averaging spikes over an exponential time window. As the number of recent
spike increases, the (V − Ek ) term begins to dominate. This driving force is trying
to shut the neuron off. Hence with more spikes, the neurons spikes less frequently.
1.3.3 Matsuoka Oscillator
The Matsuoka Oscillator [?] another popular oscillator which seems almost like a
continuous time version of the integrate and fire model with adaptation and is given
as:
τi
ui
= −ui − β f (vi ) + ∑ wi j f (u j ) + u0
dt
j6=i
vi
= −vi − f (ui )
dt
f (u) = max(0, u), (i = 1, 2)
τi0
(14)
(15)
(16)
As can be seen, as the neuron fires more, the value vi accumulates (almost like
averaging spikes, only this neuron does not spike). As vi increases, it introduces an
inhibitory term suppressing the “firing rate” of the neuron.
1.3.4 Izhikevich Model
The Hodgkin-Huxley model is capable of a rich range of behavior, but is difficult to
integrate. The leaky integrator, integrate and fire with adaptation, and the Matsuoka
oscillator are easy to integrate by have less richness than the HH.
The best of all-possible-world model was discovered by Izhikevich[?] :
dv
= 0.05v2 + 5v + 140 − u + I
dt
du
= a(bv − u)
dt
i f (v > 30mv)v− > c, u− > u + d
(17)
(18)
(19)
See Fig. 2 for the range of firing behaviors that can be obtained from this model.
These behaviors are controlled by parameters a,b,c, and d.
1.4 Numerical Integration
Next we turn to the practical matter of computing our equations. We can do so by
numerically integrating the equations on a digital computer.
6
M. Anthony Lewis
Fig. 2 Range of behavior available with the Izhikevich neural model. Electronic version of the
figure and reproduction permissions are freely available at www.izhikevich.com
The equations for neurons that we have discussed so far are of the form:
dx
= f (x)
dt
(20)
dx = f (x)dt
(21)
∆ x ≈ f (x)∆t
(22)
we can write:
How to Build Spiking CPG Models Using Python
7
xn+1 = xn + f (x)∆t
(23)
x0 = x(0)
(24)
The use of eqns (23) and (24) is called Eulers Integration. Derivation and use
of Euler Integration is straightforward, and they are completely adequate for most
equations.
However equations such as the Hodgkin-Huxely equation may benefit from a
better integrator such as Runge-Kutta 4th order. The iterative equations are given
as:
h
xn+1 = xn + (a + 2b + 2c + d)
6
a = f (tn , xn )
h
h
b = f (tn + , xn + a)
2
2
h
h
c = f (tn + , xn + b)
2
2
d = f (tn + h, xn + hc)
(25)
(26)
(27)
(28)
(29)
Runge-Kutta is more efficient than Euler Integration. While Runge-Kutta requires
more evaluations per time-step than Euler Integration, we can accurately take much
larger time steps using Runge-Kutta.
However, Euler is often fast enough for small networks and has the advantage
that it implementation is trivial.
1.5 Building Neural Oscillators: Natures coordination and
trajectory generation mechanism
Here we turn to the basic unit of movement generation in vertebrates: the oscillator.
Oscillator circuits can be seen at the core of virtually all models of rhythmic motion generation. They are the basic building blocks of the so-called Central Pattern
Generators, or CPG circuits. These circuits are groups of neurons (in vertebrates
they are in the spinal cord) seen in animals capable of periodic movement including
walking, running, swimming, flying etc.
The basic neural unit is a network oscillator composed of two neuron in mutual inhibition, see Fig. 3. The basic idea of the half-centered oscillator was first
proposed by Brown [?] and was the key element in a theory of central pattern generation, that is, the generation of rhythmic movement by spinal circuits without the
need for a sensory trigger, or a detailed signal from the brain. While this idea was
successfully repressed for period of time, in recent years it has become accepted
since Grillner’s seminal work [?, ?]. In Fig. 5 we see an illustration of the spinal
8
M. Anthony Lewis
Fig. 3 The basic element of the central pattern generator is the Brown half-centered oscillators.
Two neurons are coupled with mutual inhibition. As one neuron fires, it suppresses the other.
Eventually, the firing rate wanes and the other neurons becomes active and the cycle repeats.
Fig. 4 A Brown half centered oscillator constructed using an integrate and fire neurons with adaptation implemented with discrete components (i.e. real resistors, capacitor etc., according to the
equations given above).
circuits for locomotion in lamprey. The lamprey is a animal that first evolved some
550 million years ago and has been relatively undifferentiated since that time. It
thus gives us insight into the earliest solutions of biology. Such circuits, especially
with inhibition across the mid-line to enforce a 180 degree phase shift between the
left and right half body has been highly conserved. Whenever you walk, you take
advantage of the lamprey solution in the left and right alternation of your legs.
1.6 Reflexes and High Level Control
The CPG circuits are modulated by sensory stimuli to adapt the CPG to the animal
and to its immediate environment. Thus, the CPG is not a rigorous general given
orders to a low level control system. Rather it cooperates dynamically with the en-
How to Build Spiking CPG Models Using Python
9
Fig. 5 Basic circuitry controlling locomotion in lamprey, an eel like animal. Adapted from [?].
The C neurons are coupled across the midline with mutual inhibition. They form the basis of the
left-right bending of the lamprey.
vironment via reflexes, for short term adaptation, and descending control from the
brain, for long range anticipatory control and to achieve high level goals.
For details on the biophysics of neurons, the reader is referred to [?]. For more
details on neural models the reader is referred to [?]
2 Notable systems
Here we highlight some notable models in neurobotics/biorobotics by others
1. Case Western group– Beer and colleages have created neural models of locomotion. Starting from simulation, this group has built a robot controlled by CPG
based on the cockroach. This work is highly instructive. See [?, ?, ?, ?, ?] for
details on their work.
2. Neuromechanical model of swimming— Ekeberg has created a series of neuromechanical models of swimming in lamprey. What is interesting here is the
integration of neurons, muscles, body and environment to achieve a complex behavior, swimming [?].
3. TAGA— Through simulation, Taga and colleagues have explored the relationship between dynamics, neural networks, and the environment. [?, ?, ?, ?]
4. Tekken— Tekken is a quadrupedal robot driven by Matsuoka oscillators. It is
capable of rapid walking, and features highly innovative reflex integration with
CPGs, see [?, ?].
3 Conclusion
Neurorobotics is a paradigm which has evolved out of the desire to understand computation in the human brain. Due to advances in the theory of neural computation,
10
M. Anthony Lewis
and the dramatic increase in processing power, Neurorobotics may prove to be capable of creating highly complex behaviors in life-like machines. The market acceptance of such machines is complicated and remains uncertain, but promising.
How to Build Spiking CPG Models Using Python
11
A-1 APPENDIX
Exercises
Study neuron.py (see Appendix A-2) Several simple neuron models are included.
To integrate a particular neuron equation we can use the following commands (this
is inside the file ex1.py):
i m p o r t n e u r o n # i m p o r t t h e n e u r o n . py s c r i p t
Vm = 0 # a v a r i a b l e r e p r e s e n t i n g t h e membrane v o l t a g e a t i m e = 0
Vm t = [ ] # a l i s t c o n t a i n i n g t h e membrane v o l t a t e t h r o u g h t i m e .
params = [ 1 , 1 , 0 ] # These a r e p a r a m e t e r s t h a t can change with each
# integration
t i m e s t e p such as i n p u t c u r r e n t , e t c .
# f o r t h e LPF model o f t h e n e u r o n , t h i s p a r a m e t e r l i s t
# u s e s [ c u r r e n t , membrane c a p a c i t a n c e , N/A]
# n o t i c e here t h a t i t i s the r a t i o of I /C t h a t
# d e t e r m i n e s t h e b e h a v i o r o f t h e n e u r o n . We
# are not concerned here with the c u r r e n t
# d i m e n s i o n ( i . e . 10 c o u l d r e p r e s e n t a p i c o amp ,
# nano amp e t c . )
d e l t a t = .001
# i n t e g r a t i o n t i m e s t e p . 0 . 0 0 1 Sec ( 1 ms )
# i s u s u a l l y a p r e t t y good s t a r t i n g p o i n t
for i in range (1 ,10000):
# we c a n now i t e r a t e f o r 10 s e c t o n
Vm = n e u r o n . e u l e r (Vm, d e l t a t , n e u r o n . l p f , p a r a m s ) # c a l l E u l e r
# i n e g r a t i o n a p p e n d t h e membrane v o l t a g e f o r
Vm t . a p p e n d (Vm)
# ploting later
Note that this listing can be typed directly into ipython or can be run as a script
using:
run ex1 pylab.plot (V mt )
This will result in a plot of the membrane voltage versus time.
12
M. Anthony Lewis
Exercise 1: Neuron as a Low Pass filter The dominatant slow dynamics of a neuron can be thought of as a low pass filter.
A study the script neuron.py. Note the function lpf(v,p). We will use this as a low
pass filter model of a neuron
B Run interactive python with pylab imported : ipython -pylab
C Run ex1.py : [1] run ex1
D now plot the results : pylab.plot(V mt )
E Experiment 1: change the starting membrane voltate Vm in the script to several
different values and re-run the rpogram What do you observe?
F Experiment 2: What happens if you double or half the input current? Does
this change how long the the model takes to change potential from Vm(0) to
(Vm(final)-Vm(0)/2, i.e. 1/2 the way from the initial value to the final value?
Does changing the inital starting voltage affect this time?
Exercise 2: Minimal Model of a spiking neuron Many neurons exhibit spiking
behavior. An integrate and fire model captures both the slow dynamics as well as the
fast dynamics of the spike. To model spiking behavior, we simply reset the neuron
when it reaches given threshold.
A Run exp2.py : [1] run ex2
B Experiment 1: The input current is set to 2. Estimate the number of spikes within
10 seconds (10000 time steps).
C Experiment 2: Double the input current (to 4) . What how many spikes do you
observe in 10 seconds now? Try an input current of 8, 16 and 32. Now make a
plot of input current versus spikes. What conclusion can you make?
D Experiment 3: Is there a minimum threshold for firing?
Exercise 3: Adaptation Certain neurons exhibit adaptation. That is, given a constant current input, the neuron will initially burst rapidly and then fatigue (or adapt).
This adaptation can be exploited to build Brown Half-centered oscillators.
A Run exp3.py : [1] run ex3
B Experiment 1: The input current is set to 2. Estimate the number of spikes within
10 seconds (10000 time steps).
C Experiment 2: Double the input current (to 4) . What how many spikes do you
observe in 10 seconds now? Try an input current of 8, 16 and 32. Now make a
plot of input current versus spikes. What conclusion can you make?
D Experiment 3: Is there a minimum threshold for firing?
Exercise 4: Synapse
In neurons, presynaptic axons connect to postsynaptic dendrites via synapses.
Neurotransmitter is released into the presynaptic cleft, opening up ion channels in
the postsynaptic neuron.
We can create a minimal model of a synapse with the following equation:
I = (Vs −V ) ∗ Ap ∗W
(2)
How to Build Spiking CPG Models Using Python
13
Here Vs is the reveral potential of the ionic channels in the synaptic cleft, Ap takes
on a value of 0 or 1 and W is a coupling strenght.
A Run exp4.py : [1] run ex4
B Experiment 1: The defalult value of the tonic current input is It onic = 10 . We
introduce a new functions; “Neuron.pulses” this function produced a chain of
spikes at regular intervals. This function takes two arguments, a spike count and
a inter-spike interval. In our case, the inter-spike interval is given in ms. Initially,
we can set the intespike interval to a very large number, rendering it ineffective
(e.g. 10000). If you runt the program, the number of spikes produced in the 10
second time interval will be 25. Confirm this by running the program.
C Experiment 2: We have coupled the spikes from pulse to our neuron model. In
ths case, the coupling is inhibitory. Determine the minimum ISI needed to shut
the neuron off. What frequency does this correspond to in Hertz.
D Experiment 3: Try increasing the It onic input to 100. What spike frequency is
needed to turn the neruon off?
E Experiment 4: Lets design 1/2 of the CPG First, lets select a desired firing rate of
the neuron. Lets say that the rate will be 25 spikes per second (250 in a 10 second
period). Set the ISI to a high value to turn it off (e.g. 10000). Now determine the
It onic such that our neuron fires 250 times. At this rate the ISI is 40 ms.
Exercise 5: Brown Half-Centered Oscillator
A Run exp5.py : [1] run ex5
B TBD
14
M. Anthony Lewis
APPENDIX B
Neuron.py Listing
GL = 1 . 0
EL = 0
Ek = −2
Es p = 2
Es m = −2
gtau = 1
syntau = .05
rm = −40
W = 0
d e f l p f ( v , p ) : # e x e r c i s e 1 : Low P a s s f i l t e r
I = p [0]
# input current
Cm = p [ 1 ] # c a p a c
r e t u r n (−GL∗ ( v−EL) + I ) / Cm
def i f n (v , p ) :
I tonic = p [0]
Cm = p [ 1 ]
gsra = p [2]
Ap p = p [ 3 ] # t h i s
# both the
# Ap p i s
Ap n = p [ 4 ] # t h i s
# both the
# Ap n i s
i s a w e i g h t e d a c t i o n p o t e n t i a l and r e p r e s e n t s
sum o f many AP a s w e l l a s t h e i r r e s p e c t i v e w e i g h t s
for p o s i t i v e weights
i s a w e i g h t e d a c t i o n p o t e n t i a l and r e p r e s e n t s
sum o f many AP a s w e l l a s t h e i r r e s p e c t i v e w e i g h t s
for negative weights
W = 1
I c u r r e n t = ( Es m − v ) ∗ Ap n ∗W + ( E s p − v ) ∗ Ap p ∗W
I total
= (−GL∗ ( v−EL) + g s r a ∗rm ∗ ( v−Ek ) + I t o n i c + I c u r r e n t
)
#print I total
return
I t o t a l /Cm
def gsra (g , p ) :
r e t u r n −g / g t a u
def gsyn ( s , p ) :
r e t u r n −s / s y n t a u
How to Build Spiking CPG Models Using Python
def e u l e r ( x , dt , f , p ) :
#
p r i n t x , f (x , p ) , dt
r e t u r n x + f ( x , p )∗ dt
def p u l s e s ( index , f r e q ) :
i f i n d e x%f r e q == 0 :
return 1
else :
return 0
15
16
M. Anthony Lewis
Acknowledgements This tutorial was written at the NSF Telluride Neuromorphic Cognitive Engineering Workshop, 2009.
References
1. Beer, R., Chiel, H.J., Quinn, R., Ritzmann, R.: Biorobotic approaches to the study of motor
systems. Current Opinion in Neurobiology 8, 777–782 (1998)
2. Beer, R., R.D., Q., Chiel, H.J., Ritzmann, R.: Biologically-inspired approaches to robotics.
Communication of the ACM 40(3) (1997)
3. Beer, R.D., Chiel, H.J., Gallagher, J.: Evolution and analysis of model cpgs for walking ii.
general principles and individual variability. Journal of Computational Neuroscience 7, 119–
147 (1999)
4. Brown, T.G.: The intrinsic factors in the act of progression in the mammal. Proc. of the Royal
Soc. Lond., Ser. B. 84, 309–319 (1911)
5. Calvitti, A., Beer, R.D.: Analysis of a distributed model of a leg coordination i. individual
coordination mechanisms. Biological Cybernetics (1999)
6. Dayan, P., Abbott, L.F.: Theoretical Neuroscience: Computational and mathematical Modeling of Neural Systems. The MIT Press, Cambridge, MA (2001)
7. Delcomyn, F.: Neural basis of rhythmic behavior in animals. Science 210, 492–498 (1980)
8. Ekeberg O Grillner, S.: Simulations of neuromuscular control in lamprey swimming. Philos
Trans R Soc Lond B Biol Sci. 354(1385) (1999)
9. Gallagher, J., Beer, R.: Evolution and analysis of dynamic neural networks for agents ingtegrating vision, locomotion, and short-term memory. In: Genetic and Evolutionary Computation Conference (1999)
10. Grillner, S., Zangger, P.: On the central generation of locomotion in the low spinal cat. Exp.
Brain Res. 34, 241–61 (1979)
11. Hille, B.: Ionic Channels of Excitable Membranes. Sinauer, Sunderland (1984)
12. Hodgkin, A.L., Huxley, A.F.: A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117(4), 500–544 (1952)
13. Izhikevich, E.: Simple model of spiking neurons. IEEE Transactions on Neural Networks
14(6), 1569 – 1572 (2003)
14. Kimura, H., Fukuoka, Y., Hada, Y., Takase, K.: Adaptive dynamic walking of a quadruped
robot on irregular terrain using a neural system model. In: ISRR, pp. 88–97. Lorne, Australia
(2001)
15. Kimura, H., Fukuoka, Y., Konaga, K., Hada, Y., Takase, K.: Towards 3d adaptive dynamic
walking of a quadruped robot on irregular terrain by using neural system model. In: IEEE and
RSJ IROS 2001. IEEE, Hawaii (2001)
16. Matsuoka, K.: Mechanisms of frequency and pattern control in the neural rhythm generators.
Biol. Cybern 56, 345–353 (1987)
17. Taga, G.: A model of the neuro-musculo-skeletal system for human locomotion. i. emergence
of basic gait. Biological Cybernetics 73, 97–111 (1995)
18. Taga, G.: A model of the neuro-musculo-skeletal system for human locomotion. ii. real-time
adaptability under various constraints. Biological Cybernetics 73, 113–121 (1995)
19. Taga, G.: A model of the neuro-musculo-skeletal system for anticipatory adjustment of human
locomotion during obstacle avoidance. Biological Cybernetics 78, 9–17 (1998)
20. Taga, G., Yamaguchi, Y., Shimizu, H.: Self-organized control of bipedal locomotion by neural
oscillators in unpredictable environment. Biol. Cybern. 65, 147–159 (1991)