ar
X
iv
:q
-b
io
/0
41
10
29
v1
[
q-
bio
.N
C]
1
2 N
ov
20
04
Neurokinematic Modeling of Complex Swimming Patterns of the
Larval Zebrafish
Scott A. Hill a,b, Xiao-Ping Liu c, Melissa A. Borla c, Jorge V. Jose´ a,b,∗, and
Donald M. O’Malley a,c
aCenter for Interdisciplinary Research in Complex Systems, Northeastern University, Boston MA USA
bDepartment of Physics, Northeastern University, Boston MA USA
cDepartment of Biology, Northeastern University, Boston MA USA
Abstract
Larval zebrafish exhibit a variety of complex undulatory swimming patterns. This repertoire is controlled by
the 300 neurons projecting from brain into spinal cord. Understanding how descending control signals shape
the output of spinal circuits, however, is nontrivial. We have therefore developed a segmental oscillator model
(using NEURON) to investigate this system. We found that adjusting the strength of NMDA and glycinergic
synapses enabled the generation of oscillation (tail-beat) frequencies over the range exhibited in different larval
swim patterns. In addition, we developed a kinematic model to visualize the more complex axial bending patterns
used during prey capture.
Key words: zebrafish, swimming, locomotion, CPG, spinal cord
1. Introduction
The spinal cords of vertebrate animals contain segmental oscillators, or central pattern generators
(CPGs), that can produce rhythmic movements. The operations of these spinal CPGs are best under-
stood in lower vertebrates, such as lamprey and Xenopus, where they are used for undulatory swimming
(Roberts, 1998; Buchanan, 1999; Grillner, 2003). Control signals descending from brainstem to spinal cord
have also been studied extensively in these and other lower vertebrates, such as goldfish and zebrafish.
Their lack of a corticospinal tract avoids a degree of complexity that is present in mammals. The func-
tioning of descending control systems in higher vertebrates has been difficult to understand. Studies of
lower vertebrates should reveal conserved principles by which these systems operate. The larval zebrafish
takes vertebrate simplicity to an extreme: the decreased numbers of neurons allows exact identification of
∗ Corresponding author: jjv@neu.edu.
Preprint submitted to Neurocomputing 7 February 2008

many cell types in both brainstem and spinal cord (Kimmel et al., 1982, 1985; Liu and Westerfield, 1988;
Bernhardt et al., 1990); this in turn provides major experimental and modeling advantages.
One important aspect of lower vertebrate locomotion is the regulation of swimming speed, which is often
correlated with the frequency of alternating left and right contractions of the axial muscles of the trunk or
tail, termed tail-beat frequency (TBF). In a steady-swimming fish TBF is generally equal to the oscillation
frequency of the spinal CPGs, so by understanding the modulation of CPG frequency in spinal cord we
can understand a major element of the control of swim speed. In goldfish, lamprey, and other fishes,
swim frequency can be modulated by stimulation of the midbrain locomotor region or bath application
of NMDA; It is also known that serotonin, acetylcholine, dopamine, and other neurotransmitters also
influence the CPG’s oscillation frequency. However, the identities and locations (in brainstem or spinal
cord) of the cells involved in swim frequency control are unknown. Larval zebrafish may provide insights
into this problem because they exhibit distinct swim patterns that span a broad range of tail beat
frequencies, ranging from 25 to 75 Hz (Budick and O’Malley, 2000). A critical unknown is whether or
not the distinct larval swim patterns (slow, burst, and capture) are controlled by distinct control systems
and/or specific neurotransmitters(Borla et al., 2002; O’Malley et al., 2003). Modeling the different types
of swimming behavior using a combined neural and kinematic model can shed light onto the presence
and relevance of different motor control elements.
The functioning of spinal networks that underlie locomotion in fishes and tadpoles has been extensively
modeled (see e.g. Roberts and Tunstall, 1990; Dale, 1995, 2003; Grillner, 2003). We created a zebrafish
neural model, based on previous Xenopus spinal network models, of the CPGs in the larval zebrafish spinal
cord, incorporating known properties of the oscillators underlying swimming (Tunstall et al., 2002). The
spinal interneuron types in zebrafish (Hale et al., 2001) are likely homologous to those in both Xenopus
and lamprey (Fetcho, 1992). We explored the control of tail-beat frequency, and found that by altering
the strengths of AMPA, NMDA and glycinergic-like synapses (all known to be present in lower vertebrate
spinal cords), we were able to generate TBFs that spanned the range of speeds observed during burst
and slow swimming behaviors. We also created a simple mechanical or “kinematic” model (that can be
driven by the neural model) to visualize how the spinal neural activity might be transformed into larval
behaviors. Our ultimate goal in creating this neurokinematic model is to introduce a tool for testing
theories of descending motor control in the larval zebrafish.
2. Methods
Figure 1 shows the structure of our neural model, which is a simplification of the larval spinal cir-
cuitry (Hale et al., 2001) incorporating the minimal elements required to generate rhythmic, propagating
alternating activity. Each segmental oscillator consists of two neurons, with NMDA and AMPA-like au-
tosynapses, which are connected by glycinergic synapses. An oscillator can be triggered with a single
excitatory pulse to one cell, followed by a single pulse to the other cell several milliseconds later (7 ms
in our simulations). Individual oscillators are connected into a 25-segment chain (corresponding to the
approximately 25 segments in the zebrafish spinal cord), with nearest-neighbor descending excitatory
synaptic connections; the entire chain can be started merely by triggering the head segment.
We use cells with the standard Hodgkin-Huxley channels provided by the NEURON software package:
persistent potassium (g¯K = 0.036 S/cm2, EK = −77mV), transient sodium (g¯Na = 0.12 S/cm2, ENa =
50mV), and leak (gL = 0.0003 S/cm2, EL = −54.3mV) channels. Our synapses are modeled using the
difference of two exponentials, with rise time τ1(= 1ms) and fall time τ2. See Dayan and Abbott (2001,
p. 182) for details. The three types of synapses we used were AMPA (τ2 = 6ms, E = 0mV), NMDA
(τ2 = 80ms, E = 60mV), and glycine (τ2 = 2ms, E = −80mV).
2

Fig. 1. Schematic diagram of the neural oscillator model. The squares represent neurons, and each line indicates one of three
types of synaptic connection: glycinergic (fast inhibitory), AMPA-like (fast excitatory), or NMDA-like (slow excitatory).
Only the first two segments of a 25-segment chain are shown.
Our kinematic model is meant to translate the neural signals from the spinal circuitry into observable
kinematic behaviors. To do so, we use a relatively simple transformation of putative neural output into
effects on the radius of curvature of a line segment representing the trunk of the larva. This model does
not take into account the full physics (e.g., elasticity and hydrodynamics) of the situation; nevertheless,
it does approximate the observed axial kinematics recorded experimentally with a high-speed camera.
Suppose that each segment x is receiving a neural signal Fs(x, t) at time t from the spinal circuitry:
Fs > 0 for a signal to the right side of the segment, Fs < 0 for a signal to the left. We suppose that this
signal is integrated through exponential synapses, so that the signal passing to the muscles is
Fm(x, t) =
t
∫
−∞
Fs(x, t′)[e(t
′−t)/τ2 − e(t
′−t)/τ1 ]dt′, (1)
where τ1 and τ2 are the growth and decay time constants of the synapse (we use τ1 = 6ms and τ2 = 8ms
in our calculations.) Assuming that the muscle contracts linearly as a function of Fm, it follows that the
radius of curvature of segment x is
R(x, t) = W (x)/Fm(x, t) (2)
W (x) is a function describing the stiffness, or resistance to bending, of segment x. It takes into account
the width of the body (the tail is more flexible, and bends more than the rostral trunk).
We can feed the output of our neural model directly into this system, but it is also useful to introduce
a (more) artificial signal. We build this signal out of three components: an oscillatory signal Fosc(x, t),
which is a series of delta functions propagating caudally; a bending signal Fbend, which is a tonic signal
applied to one side of all segments at once; and a rostral stiffening signal, which reduces the signal to
the first xinh rostral segments by a factor finh. The first two pieces are responsible for swimming and
turning, respectively. The stiffening signal is seen in larvae during prey capture, where the fish keeps its
head relatively still while adjusting its orientation with its tail. We can write our artificial neural signal
as
Fs(x, t) = (Fosc(x, t) + fbend)×
{
finh, x ≤ xinh
1, x > xinh
, (3)
where fbend is a constant and
3

Fig. 2. Instantaneous TBF in a single segmental oscillator, for several values of gNMDA. Each point represents the iTBF for
the left neuron; the corresponding points from the right side gives near-identical results. Oscillators approach a steady-state
frequency only after an initial “wind-up” period.
Fosc(x, t) = fosc



1, x ≡ 2piνt (mod λ)
−1, x ≡ 2piνt + λ/2 (mod λ)
0, otherwise
(4)
The parameter ν is the tail-beat frequency, and λ is the length (in number of segments) of the resulting
wave which propagates down the fish.
3. Results and Discussion
We began by studying a single segment of the neural model shown in Fig. 1. A transient trigger
pulse to this model initiates sustained alternating activity, mimicking Xenopus spinal cord, where a
transient stimulus produces a sustained bout of swimming (Roberts and Tunstall, 1990). The first larval-
zebrafish specific task was to generate the range of tail-beat frequencies (TBFs) used in different swimming
behaviors. The frequency of oscillation was calculated on a half-cycle basis (termed “instantaneous”
TBF in Borla et al., 2002). Figure 2 shows how instantaneous TBF varies with time, with each trace
representing the frequency profile for a specific value of the NMDA synaptic conductance. In all cases,
there is an initial transient period, 15 to 20 cycles long, where the TBF increases by about a third before
settling into an indefinite steady state with constant frequency. We subsequently re-inspected behavioral
sequences, and found varying degrees of wind-up in some, but not all, experimentally recorded swim bouts
(Borla et al., 2002). It is unknown whether the biological and simulational wind-up behaviors are related,
but the phenomenon lasts much longer in the model than in the lab.
To evaluate the range of oscillator or tail-beat frequencies that might be generated with this simple
model, we tested different combinations of synaptic strengths (conductances). Varying the strengths of
the NMDA and glycinergic synapses was found to alter the steady-state oscillator frequency (Fig. 3):
increasing the NMDA conductance increased TBF, whereas increasing the glycinergic conductance de-
creased it. The AMPA conductance can also affect TBF, but to a lesser extent. Different combinations
of conductances could give rise to the entire range of TBFs observed in different larval swim patterns,
from the slow swim (25 to 40 Hz) to burst swim (45 to 75 Hz; Budick and O’Malley, 2000). This is just
4

Fig. 3. A density plot showing the steady-state frequency of a single oscillator as we modify the synaptic conductances
gNMDA and gINH, fixing gAMPA = 10−4. All conductances are in µS, and all frequencies are in Hz. The missing regions in
the upper-right and lower-left corners represent non-oscillatory states.
Fig. 4. Rastergram for a chain of 25 segmental oscillators. Each horizontal line corresponds to a single oscillator; ticks drawn
below each line are for action potentials on the left, and above for those on the right. Parameters: gAMPA = 10−4 µS,
gNMDA = 6× 10−4 µS, gINH = 10−2 µS, gdescending = 10−2 µS.
one potential means of varying oscillator frequency (or TBF), implemented in a reduced system (i.e. a
two-cell model), but it illustrates how a minimal model of larval zebrafish spinal cord, with just a few
essential conductances, can give rise to a range of outputs relevant to the larval behaviors. In this model,
each crossed inhibitory signal results in one post-inhibitory rebound firing of a single action potential, but
to more completely capture slow and burst swims we will need to incorporate mechanisms that regulate
the strength of output of the motoneuron pools.
We next extended the model to a 25-segment chain of oscillators, creating a kind of artificial spinal
cord that should facilitate quantitative analyses of the influences of descending signals on spinal network
activity. The firing pattern of the model spinal cord is illustrated in Fig. 4. The firing of segment num-
ber 1 shows an alternating left-right pattern over the duration of the simulation. In successively more
caudal segments, the firing is delayed, matching the general pattern of undulatory swimming. As in the
single segment model, each segment along the chain showed a frequency wind-up period. In this model a
fixed intersegmental time-delay (nominally a synaptic delay) was used to establish the phase relationship
between segments, but a more realistic model might include ascending and descending connections of
varying strengths and lengths, as in lamprey (Kotaleski et al. 1999, and see Dale 2003).
The kinematic model of the larval trunk is a complementary tool for exploring theories of descending
motor control. Preliminary work (Hill et al., 2003) assumed that the fish was equally flexible from trunk
5

Fig. 5. Comparison of a 7-day old larval zebrafish with our kinematic model during a slow swim (left) and a J-turn (right).
to tail (W (x) = 1 in Eq. 2), with partial success. By accounting for the increased flexibility of the tail
(e.g. by making W (x) linear), however, we can better simulate the trunk kinematics observed during
slow swims and J-turns (Fig. 5). J-turns are unique locomotive maneuvers that contribute to the larval
prey-capture behavior (Borla and O’Malley, 2002). They require repetitive, asymmetric and far-caudal
contractions of axial musculature. This might in principle be achieved by sending an asymmetric excitatory
brainstem signal exclusively to far-caudal spinal cord; however, a survey of the spinal outputs of zebrafish
reticulospinal neurons revealed no neurons with the requisite arborization pattern (Gahtan and O’Malley,
2003). Alternatively, neurons that selectively arborize in rostral spinal cord (which were observed) might
be activated bilaterally to stiffen the rostral musculature. In conjunction with such signals, other neurons
that arborize along the entire length of spinal cord could then generate far-caudal contraction. In our
model, this is implemented as a bilateral “inhibitory” rostral signal, but this is kinematically equivalent
to stiffening the rostral end of the larva by bilateral excitation of rostral motoneurons.
Future Directions
As we extend the neural model’s capabilities to produce dynamically varying bend amplitudes, we
hope to produce an increasingly realistic “artificial spinal cord” that can be used to test ideas of how the
larval locomotive repertoire is generated. The kinematic model is useful here because it provides a first
approximation of how complex neural model outputs might affect axial kinematics. Further extensions
of this model should incorporate known constraints of the larval CNS, such as the diversity of brainstem
and spinal systems involved in swimming and turning behaviors (O’Malley et al., 1996; Hale et al., 2001;
Ritter et al., 2001), and the wide distribution of activity during escape behaviors (Bosch et al., 2001;
Gahtan et al., 2002). By incorporating such constraints, the combined neurokinematic model should
become increasingly useful in generating experimentally-testable hypotheses of descending and spinal
control of vertebrate locomotion.
4. Acknowledgements
Support for this work was provided by NIH-NS37789 (DMO) and by the Center for Interdisciplinary
Research in Complex Systems (CIRCS) at Northeastern University (SAH, JVJ).
5. References
R.R. Bernhardt, A.B. Chitnis, L. Lindamer, J.Y. Kuwada, Identification of spinal neurons in the embry-
onic and larval zebrafish, J. Comp. Neurol. 302 (1990) 603-616.
M.A. Borla and D.M. O’Malley, High-speed imaging of tracking swims used in the larval zebrafish prey
capture behavior, Soc. Neurosci. Abs. 28 (2002) 361.18.
6

M.A. Borla, B. Palecek, S.A. Budick, and D.M. O’Malley, Prey capture by larval zebrafish: evidence for
fine axial motor control, Brain Behav. Evol.. 60 (2002) 207-229.
T.J. Bosch, S. Maslam, B.L. Roberts, Fos-like immunohistochemical identification of neurons active during
the startle response of the rainbow trout, J. Comp. Neurol. 439 (2001) 306-314.
J.T. Buchanan, Contributions of identifiable neurons and neuron classes to lamprey vertebrate neurobi-
ology, Prog. Neurobiol. 63 (1999) 441-466.
S.A. Budick and D.M. O’Malley, Locomotive repertoire of the larval zebrafish: swimming, turning, and
prey capture, J. Exp. Biol. 203 (2000) 2565-2579.
P. Dayan and L.F. Abbott, Theoretical Neuroscience, MIT Press, Cambridge MA, 2001.
N. Dale, Experimentally derived model for the locomotor pattern generator in the Xenopus embryo,
J. Physiol. 489.2 (1995) 489-510.
N. Dale, Coordinated motor activity in simulated spinal networks emerges from simple biologically plau-
sible rules of connectivity, J. Comput. Neurosci. 14 (2003) 55-70.
J.R. Fetcho, The spinal motor system in early vertebrates and some of its evolutionary changes. Brain Be-
hav. Evol. 40 (1992) 82-97.
E. Gahtan and D.M. O’Malley, Visually-guided injection of identified reticulospinal neurons in zebrafish,
J. Comp. Neurol. 459 (2003) 186-200.
E. Gahtan, N. Sankrithi, J.B. Campos, and D.M. O’Malley, Evidence for a widespread brainstem escape
network in larval zebrafish, J. Neurophysiol. 87 (2002) 608-614.
S. Grillner, The motor infrastructure: from ion channels to neuronal networks, Nature Rev. Neurosci. 4
(2003) 573-586.
M.E. Hale, D.A. Ritter, and J.R. Fetcho, A confocal study of spinal interneurons in living larval zebrafish,
J. Comp. Neurol. 437 (2001) 1-16.
S.A. Hill, M.A. Borla, J.V. Jose´, and D.M. O’Malley, Modeling the neural control of zebrafish locomotive
behaviors, Soc. Neurosci. Abs. 29 (2003) 278.10.
C.B. Kimmel, S.L. Powell, and W.K. Metcalfe, Brain neurons which project to the spinal cord in young
larvae of the zebrafish, J. Comp. Neurol. 205 (1982) 112-127,
C.B. Kimmel, W.K. Metcalfe, and E. Schabtach, T-reticular interneurons: a class of serially repeating
cells in the zebrafish hindbrain, J. Comp. Neurol. 233 (1985) 365-376.
J.H. Kotaleski, S. Grillner, and A. Lansner, Neural mechanisms potentially contributing to the interseg-
mental phase lag in lamprey: I. Segmental oscillations dependent on reciprocal inhibition, Biol. Cybern.
81 (1999) 317-330.
D.W. Liu and M. Westerfield, Function of identified motoneurons and coordination of primary and sec-
ondary motor systems during zebrafish swimming, J. Physiol. (Lond.) 403 (1988) 73-89.
D.M. O’Malley, Y.-H. Kao, and J.R. Fetcho, Imaging the functional organization of zebrafish hindbrain
segments, Neuron 17 (1996) 1145-1155.
D.M. O’Malley, Q. Zhou, and E. Gahtan. Probing Neural Circuits in the Zebrafish: A Suite of Optical
Techniques. Methods 30 (2003) 49-63.
D.A. Ritter, D.H. Bhatt, and J.R. Fetcho, In vivo imaging of zebrafish reveals difference in the spinal
networks for escape and swimming movements, J. Neurosci. 21 (2001) 8956-8965.
A. Roberts, S.R. Soffe, E.S. Wolf, M. Yoshida, and F.-Y. Zhao, Central circuits controlling locomotion in
young frog tadpoles, Annals NY Acad. Sci. 860 (1998) 19-34.
A. Roberts and M.J. Tunstall, Mutual re-excitation with post-inhibitory rebound: a simulation study on
the mechanisms for locomotor rhythm generation in the spinal cord of Xenopus embryos, Euro. J. Neu-
rosci. 2 (1990) 11-23.
M.J. Tunstall, A. Roberts, S.R. Soffe, Modelling inter-segmental coordination of neuronal oscillators:
synaptic mechanisms for uni-directional coupling during swimming in Xenopus tadpoles, J. Com-
put. Neurosci. 13 (2002) 143-158.
7