The spinal cord is a column of nerve tissue that connects the brain to the rest of the body in vertebrate animals. Nerve cells in the spinal cord, called neurons, help to control and coordinate the body’s movements. As the spinal cord develops, new neurons are born and new connections are made between neurons and muscles, resulting in more coordinated and skillful movements as time goes on.

Zebrafish, for example, display body-bending maneuvers called coils within 24 hours of the egg being fertilized. Next, bursts of swimming movements emerge, which are driven by sporadic tail beats. These tail maneuvers become more consistent as the fish develops, and eventually result in smooth movements called beat-and-glide swimming. The groups of spinal cord neurons that appear at each stage of zebrafish development have been characterized, but it remains unclear how newly formed circuits (groups of neurons recently connected to each other) work together to produce swimming maneuvers.

To answer this question, Roussel et al. simulated changes in the spinal cord that help zebrafish acquire new swimming movements as they grow. The computer models encoded neural circuits based on cell populations identified in experimental studies, and replicated swimming behaviors that emerge during the first few days of zebrafish development. Simulations tested how specific neural circuits generate the characteristic swimming movements that represent key developmental milestones in zebrafish.

The results showed that adding new neurons and more cell-to-cell connections led to increasingly sophisticated swimming maneuvers. As the zebrafish spinal cord matured, the fish were better able to control the pace and duration of their swimming movements. Roussel et al. also identified specific patterns of neural activity linked to particular maneuvers. For example, tail beats switch direction when neurons on one side of the spinal cord excite neurons on the opposite side. This activity, which becomes more rhythmic, also needs to be exquisitely timed to produce and coordinate the right motion.

Roussel et al.’s modelling of developmental milestones in growing zebrafish provides insights into how neural networks control movement. The computer models are among the first to accurately reproduce swimming behaviors in developing zebrafish. More experimental data could be added to the models to capture the full range of early zebrafish movements, and to further investigate how maturing spinal cord circuits control swimming. Since zebrafish and mammals have many spinal neurons in common, further research may aid our understanding of movement disorders in humans.

Movements made in the early stages of development can be critical for the survival of many species. The escape response seen in various fish and amphibians is one such example of a vital movement present at early developmental stages (Domenici and Hale, 2019). However, the nervous system’s control of movement does not come fully formed but matures as the nervous system develops (Favero et al., 2014). This maturation enables a broader repertoire of movements to arise. During this process, new neurons are born and subsequently integrated into neural circuits that are newly formed or refined, presumably leading to the emergence of progressively more coordinated and skillful maneuvers. Determining how the assembly of new circuits leads to the emergence of new movements can provide valuable insights into the role of distinct neurons or circuits in motor control.

The maturation of swimming in developing zebrafish has been well described at both the ethological and the cellular levels (Drapeau et al., 2002; McLean and Fetcho, 2009). Single strong body bends on one side of the body, also known as coils, emerge during the first day of development at around 17 hr post-fertilization (hpf) as the earliest locomotor behavior (Saint-Amant and Drapeau, 1998). Single coils are quickly followed by double coils (i.e., two successive coils, one for each side of the body) at around 24 hpf (Knogler et al., 2014). Touch-evoked swimming appears around 27 hpf as coiling begins to subside. Spontaneous swimming movements emerge around 2–3 days post-fertilization (dpf) (Saint-Amant, 2010). The first swimming movement zebrafish exhibit is burst swimming characterized by long (1 s long) but infrequent episodes of tail beats. Burst swimming is then replaced by beat-and-glide swimming characterized by shorter (several hundreds of milliseconds long) but more frequent episodes. In both cases, swim episodes consist of repetitive left-right alternating, low-amplitude tail beats that propagate from the rostral toward the caudal end of the fish body and are generated at 20 to 80 Hz (Budick and O'Malley, 2000; Buss and Drapeau, 2001).

During this rapid series of transitions between locomotor maneuvers, populations of spinal neurons are progressively generated, starting with primary motoneurons (MNs) at about 9 hpf. Subsequently, spinal MNs and interneurons (INs) are generated in stereotyped spatiotemporal birth orders (Kimmel et al., 1994; Myers et al., 1986; Satou et al., 2012). Two successive waves of axogenesis occur in the embryonic spinal cord (Bernhardt et al., 1990). The first wave occurs around 16–17 hpf. It includes axon growth in primary MNs that innervate red and white muscle fibers at early developmental stages (Buss and Drapeau, 2000). Primary MNs enable coiling and escape movements (Kimmel et al., 1995; Saint-Amant and Drapeau, 2000). Several spinal INs that are also important for early movements extend their axons along with primary MNs. These include Ipsilateral Caudal (IC) INs that are thought to play an essential role in driving the rhythm of early locomotor behavior due to their endogenous bursting activity (Tong and McDearmid, 2012). The second wave of axon growth occurs at around 23–25 hpf. It involves axon growth in secondary MNs involved with slower movements (Liu and Westerfield, 1988) and spinal IN populations that include excitatory and inhibitory, ipsilaterally and contralaterally, and ascending and descending projecting subtypes (Bernhardt et al., 1990; Higashijima et al., 2004). The progressive generation of new populations of spinal neurons and continued axonal growth coincides with the expansion of the zebrafish locomotor repertoire. This timing suggests that incorporating spinal circuits into existing locomotor circuits underlies the acquisition of novel locomotor maneuvers.

We have recently provided evidence that the maturation from coiling to later stages of swimming is accompanied by an operational switch in how spinal locomotor circuits generate the rhythm underlying tail beats. Specifically, we demonstrated that spinal circuits transitioned from relying upon pacemakers with endogenous bursting properties during coiling toward depending upon network oscillators whose rhythm is driven by excitatory and inhibitory synapses (Roussel et al., 2020). In light of these and earlier findings describing the composition and maturation of spinal locomotor circuits, we sought to generate computational models that replicate developmental locomotor movements of the zebrafish. We iteratively constructed models for several locomotor movements by incorporating specific spinal populations, shifts in relative connection strength, and changes in the firing behavior of neurons. While computational modeling has generated invaluable insights into the function and mechanisms of spinal locomotor circuits of several species (Ausborn et al., 2019; Bicanski et al., 2013; Danner et al., 2019; Ferrario et al., 2018; Hull et al., 2016; Kozlov et al., 2014; Sautois et al., 2007), there is to our knowledge no such model for the developing zebrafish spinal cord. Here, we build some of the first computational models of the zebrafish spinal locomotor circuit that can accurately reproduce predominant locomotor behaviors during early zebrafish development. In the process, we test theories about the possible contributions of specific neural circuits and spinal populations to locomotor movements in zebrafish and identify untested hypotheses on the operation of spinal locomotor networks in developing zebrafish.

We aimed to model how new locomotor movements may emerge from the integration of spinal INs and the modification of synaptic weights and firing behavior over the first few days of development in the zebrafish. Our approach was to build an initial model based upon previously reported experimental observations of spinal circuits when the first locomotor movements emerge in zebrafish around 1 dpf. We then successively built upon this initial model to replicate several locomotor maneuvers of the developing zebrafish.

The models were composed of single-compartment neurons whose firing dynamics were determined by a small set of differential equations (Izhikevich, 2007). The firing of MNs was converted to muscle output. This output was used to estimate body angle and locomotor activity during simulations (Figure 1). The composition of each model depended on the developmental stage and the locomotor movement to be generated.

Figure 1

Download asset Open asset

Single coiling (>17 hpf) results from unilateral gap junction coupling

Coiling, which is already observed at 1 dpf, is characterized by a single strong, slow (hundreds of milliseconds in duration) tail beat on one side of the body followed by a return to resting position (Saint-Amant and Drapeau, 1998). Coiling events are relatively infrequent, reaching a maximum frequency of 1 Hz around 20 hpf (Saint-Amant and Drapeau, 1998). Previous studies have established that this behavior is generated by a spinal circuit relying primarily on gap junctions (i.e., electrical synapses) (Saint-Amant and Drapeau, 2001). It has been proposed that rostrally located IC pacemaker spinal neurons (Tong and McDearmid, 2012) drive periodic depolarizations (PDs) of ipsilateral MNs via electrical synapses (Drapeau et al., 2002; Saint-Amant and Drapeau, 2001). Glycinergic synaptic bursts (SBs) are observed in MNs during contralateral coiling events (Saint-Amant and Drapeau, 2001). These SBs have been proposed to arise from contralaterally projecting glycinergic neurons (Saint-Amant and Drapeau, 2000) but are not responsible for any action potential firings or coiling movements (Saint-Amant and Drapeau, 2001). Applying a gap junction blocker, heptanol, but not glutamatergic and glycinergic antagonists, suppressed spinal activity responsible for coiling (Saint-Amant and Drapeau, 2000).

Network description

(Figure 2A) Based on the experimental observations reported above, the model for single coiling consisted of rostrocaudal chains of electrically coupled spinal neurons driven by a kernel of five recurrently connected pacemakers (IC neurons). One chain consisted of 10 MNs. The other chain consisted of 10 contralaterally projecting commissural inhibitory neurons. Neurons from the V0d population are active during large amplitude movements such as escapes (Satou et al., 2020), and so we assumed V0ds were the commissural inhibitory neurons active during coiling, which is another large amplitude movement. We selected an IC kernel size of five as a trade-off between computational simplicity and robustness of the kernel to the failure of firing of a small number of cells. Similarly, the size of the coiling model was set to 10 somites. Thus, each model somite represents approximately three biological somites. This choice was made as a trade-off between computational simplicity and recreating the kinematics of coiling fish.

Figure 2 with 4 supplements see all

Download asset Open asset

Figure 2—source data 1

P-values for sensitivity testing in single coiling model.

https://cdn.elifesciences.org/articles/67453/elife-67453-fig2-data1-v2.xlsx

Download elife-67453-fig2-data1-v2.xlsx

IC neurons have been reported to project caudally through multiple somites (Bernhardt et al., 1990). Therefore, in addition to their recurrent connections, each IC formed electrical synapses with several rostral MNs and V0ds (the first four of each ipsilateral chain in our model). Electrical coupling between many populations of early born spinal neurons has been previously demonstrated, including between IC and MNs (Saint-Amant and Drapeau, 1998). Coupling between IC neurons and commissural inhibitory neurons has not been demonstrated yet. We based this electrical coupling between ICs and V0ds on the fact that glutamatergic blockers do not block glycinergic SBs present at this stage (Saint-Amant and Drapeau, 2001), suggesting that gap junctions mediate the activation of V0ds underlying these glycinergic bursts. Gap junction weights are found in Table 1 .

Table 1

Coiling Beat-and-glide	MN	IC	V0d dI6	V0v	V2a
MN Single coiling Double coiling Beat-and-glide (all models)	0.1 0.07 0.005
IC Single coiling Double coiling	0.04 0.03	0.001 0.0001
V0d or dI6 Single coiling Double coiling Beat-and-glide (all models)	0.01 0.0001 0.0001	0.05 0.05	0.04 0.04 0.04
V0v Double coiling Beat-and-glide (all models)	0.0001 0.005	0.0005		0.05 0.05
V2a Double coiling Beat-and-glide (all models)	0.005 0.005	0.15			0.005 0.005

The connectivity within the chains was identical for both MN and V0d chains. Each neuron in a chain formed electrical synapses with its three nearest rostral and caudal neighbors within the same chain. There was also electrical coupling across the two ipsilateral chains as MNs formed gap junctions with the three nearest rostral and three nearest caudal V0ds and vice-versa. Paired recordings of MNs and V0ds at this stage have yet to be published. Our assumption that MNs and V0ds are electrically coupled at this stage was based upon the widespread electrical coupling between ipsilateral spinal neurons (Saint-Amant and Drapeau, 2001). To reproduce the glycinergic bursts observed in MNs at this stage (Saint-Amant and Drapeau, 2001), V0ds projected to contralateral MNs. Thus, V0ds formed glycinergic synapses with contralateral MNs and ICs. The reversal potential of glycinergic synapses (E_gly) is depolarized during development (Ben-Ari, 2002) and was set to −45 mV in the single coiling model (see Table 2). All V0ds sent ascending projections to contralateral ICs. V0ds projected to contralateral MNs within five to six segments so that the ith V0d projected to all contralateral MNs between the i−5 and i+5 segments. Chemical synaptic weights are found in Table 3.

Table 2

Chemical synapse	Erev	τr	τf	Vthr
Glutamatergic	0	0.5	1.0	−15
Glycinergic	−45, –58, −70*	0.5	1.0	−15

1. ^*For single and double coiling and beat-and-glide swimming models, respectively.

Table 3

Pre-synaptic	Post-synaptic
	MN	IC	V0d dI6	V0v	V2a	V1	Muscle
MN Single coiling Double coiling Beat-and-glide (all models)							0.015 0.02 0.1
V0d Single coiling Double coiling	0.3 2.0	0.3 2.0			2.0
dI6 Beat-and-glide (base) Beat-and-glide (bursting V2a) Beat-and-glide (all tonic neurons)	1.5 1.5 1.5		0.25* 0.25* 0.25*		1.5 2.0 1.5
V0v Double coiling Beat-and-glide (base) Beat-and-glide (bursting V2a) Beat-and-glide (all tonic neurons)		0.175			0.4 0.75 0.4
V2a Double coiling Beat-and-glide (base) Beat-and-glide (bursting V2a) Beat-and-glide (all tonic neurons)	0.5 0.5 0.5		0.5 0.75 0.5	0.04 0.3 0.275 0.25	0.3 0.3 0.3	0.5 0.5 0.5
V1 Beat-and-glide (base) Beat-and-glide (bursting V2a) Beat-and-glide (all tonic neurons)	1.0 1.0 1.0		0.2 0.2 0.2	0.1 0.1 0.1	0.5 0.5 0.6

1. ^*Scaled by a random number selected from a gaussian distribution with mean of 1 and variance of 0.1.

Each neuron was modeled as a single compartment neuron with subthreshold and suprathreshold membrane potential dynamics described by a small set of differential equations (Izhikevich, 2007). These equations have nine parameters: a, b, c, d, and V_max (which respectively represent the time scale of the recovery variable u, the sensitivity of u to the subthreshold variation of V, the reset value of V after a spike, the reset value of u, and the action potential peak), and k, C, V_r, and V_t (coefficient for the approximation of the subthreshold part of the fast component of the current-voltage relationship of the neuron, cell capacitance, resting membrane potential, and threshold of action potential firing). Parameter values of ICs (see Table 4 for all neuron parameters) were chosen such that they exhibited a relatively depolarized threshold of action potentials and bursts of short action potentials lasting hundreds of milliseconds as seen in experimental recordings in embryonic zebrafish (Tong and McDearmid, 2012). They were also modeled to exhibit periodic bursts lasting hundreds of milliseconds in response to a constant tonic drive (Figure 2B). This firing pattern was generated in part by having a low value of a and a relatively depolarized value of c. MNs (Drapeau et al., 1999) and V0ds were modeled to generate tonic repetitive firing in response to a step depolarization (Figure 2B). Finally, to activate the circuit, a constant drive was provided to the left ICs only. Restricting the drive to left ICs prevented the appearance of near-coincident bilateral coils that could be misinterpreted as spinally mediated multiple coils.

Table 4

Model population	a	b	c	d	Vmax	Vr	Vt	k	C	Rostro-caudal position*	Idrive†
MN Single coiling Double coiling Beat-and-glide	0.5 0.5 0.5	0.1 0.1 0.01	−50 −50 −55	0.2 100 100	10 10 10	−60 −60 −65	−45 −50 −58	0.05 0.05 0.5	20 20 20	5.0+1.6*n
IC Single coiling Double coiling	0.0005 0.0002	0.5 0.5	−30 −40	5 5	0 0	−60 −60	−45 −45	0.05 0.03	50 50	1.0	50 35
V0d Single coiling Double coiling	0.5 0.02	0.01 0.1	−50 −30	0.2 3.75	10 10	−60 −60	−45 −45	0.05 0.09	20 6	5.0+1.6*n
dI6 Beat-and-glide (all models)	0.1	0.002	−55	4	10	−60	−54	0.3	10	5.1+1.6*n
V0v Double coiling Beat-and-glide (base) Beat-and-glide (bursting V2a and all tonic models)	0.02 0.01 0.1	0.1 0.002 0.002	−30 −55 −55	11.6 8 4	10 10 10	−60 −60 −60	−45 −54 −54	0.05 0.3 0.3	20 10 10	5.1+1.6*n
V2a Double coiling Beat-and-glide (base and all tonic models) Beat-and-glide (bursting V2a model)	0.5 0.1 0.01	0.1 0.002 0.002	−40 −55 −55	100 4 8	10 10 10	−60 −60 −60	−45 −54 −54	0.05 0.3 0.3	20 10 10	5.1+1.6*n	2.89 3.05
V1 Beat-and-glide (all models)	0.1	0.002	−55	4	10	−60	−54	0.3	10	7.1+1.6*n

1. ^*n=0 to N−1, N being the total number of neurons in that given population.

   ^†Amplitude of tonic motor command drive.

Simulations results

Our simulations show that this model can generate single coils characterized by large body bends to one side of the body lasting approximately 1 s (Figure 2C, Figure 2—video 1). Our base single coiling model generated six evenly interspersed single coils during a 10 s simulation. This 0.6 Hz coiling frequency is within the 0–1.0 Hz range of frequencies observed during zebrafish development (Saint-Amant and Drapeau, 1998; Saint-Amant and Drapeau, 2000). Silencing ICs blocked activity in all spinal neurons (Figure 2—figure supplement 1A), emphasizing the central role of the IC kernel in the generation of single coils.

Previously reported whole-cell patch-clamp recordings of MNs at this developmental stage display two types of events (Saint-Amant and Drapeau, 2000; Saint-Amant and Drapeau, 2001): PDs via electrical synapses and SBs from contralateral spinal glycinergic neurons that are depolarizing at rest due to the depolarized chloride reversal potential observed early in development. These events last hundreds of milliseconds. In our model, SBs were observed in the contralateral ICs and MNs (events during coilings in left neurons seen in the gray traces in Figure 2C). SBs were caused by glycinergic input from V0ds activated during the ipsilateral coilings. As observed experimentally (Saint-Amant and Drapeau, 2001), preventing SBs by silencing glycinergic synapses from V0ds did not preclude the generation of single coiling, nor did it lead to the generation of multiple coilings (Figure 2—figure supplement 1C). PDs can be unmasked by hyperpolarizing MNs sufficiently to prevent the firing of action potentials (Figure 2D). An analysis of the phase delays between ipsilateral neurons during single coils shows that IC neuron firing precedes ipsilateral MN and V0d firing (Figure 2E) and reinforces that ICs drive single coiling events.

To further validate the model, we tested whether the model could still generate single coils with different parameters. First, we tested whether the model could still generate single coils when the number of model somites was increased from 10 to 30 to be closer to the number of biological somites in zebrafish (Stickney et al., 2000). A 30-somite model with IC axons extending to all somites and several modified gap junction weights (Table 1) generated single coils (Figure 2—figure supplement 2).

Next, the base model's sensitivity to within-model parameter variability was tested. Variability in the amplitude of the tonic motor command, the rostrocaudal extent of every axonal projection, every parameter that set the dynamics of the membrane potential of each neuron (a, b, c, d, and V_max, k, C, V_r, and V_t), and all of the weights of gap junction and chemical synapses were modeled by scaling each value by a random number picked for each simulation. The random numbers were derived from a Gaussian distribution with mean =1, and standard deviations, ${\sigma }_d$ (tonic drive), ${\sigma }_l$ (rostrocaudal length of axonal projections), ${\sigma }_p$ (dynamics of membrane potential), and ${\sigma }_w$ (synaptic weights), respectively. Ten 20-s long simulations were run at various values of ${\sigma }_d$, ${\sigma }_l$, ${\sigma }_p$, and ${\sigma }_w$. In each simulation, the variability of only one of the four sets of parameters (amplitude of motor drive, length of axonal projection, membrane potential dynamics, and synaptic weights) was tested, and the standard deviations of the three other sets of parameters were set to 0.

The single coiling model's suitability was assessed by the relative absence of truncated coils, which were movements with only partial contractions restricted to the body’s rostral segments (Figure 2—video 2). We sought to determine the upper limit of variability within which the single coiling model remained suitable. For this reason, the ranges of ${\sigma }_d$, ${\sigma }_l$, ${\sigma }_p$, and ${\sigma }_w$ that were tested differed amongst the four sets of parameters tested (Figure 2F-M). A comparison of the level of variability at which the models start generating more varying frequency of coiling and more truncated coils suggests that the single coiling model is more robust to noise in the amplitude of the tonic motor command (Figure 2F, J) and was most sensitive to variability in the parameters governing the dynamics of the membrane potential (Figure 2H, L). The single coiling model was relatively mildly sensitive to variability in the synaptic weights and the rostrocaudal extent of the axon projections (Figure 2G, I, K, M).

Overall, the model replicated this first locomotor behavior of zebrafish in terms of the duration and frequency of coiling events as well as synaptic events of MNs. We then built upon this model to replicate the next step in the development of locomotion: the appearance of double coiling.

Double coiling (>24 hpf) depends on the timing and strength of contralateral excitation and inhibition

After single coils appear, double coils emerge as a transitory locomotor behavior at around 24 hpf, coexisting with the single coiling behavior (Knogler et al., 2014). Double coiling is characterized by two successive coils, one on each side of the body, and lasts about 1 s (Knogler et al., 2014). Eventually, double coiling becomes the predominant coiling behavior. Double coiling can represent nearly three-quarters of all coiling events at its peak frequency, with the rest mainly being single coils (Knogler et al., 2014).

At the stage when double coiling appears (24 dpf), the previous electrical scaffold for single coils seems to be supplemented with chemical glutamatergic synapses to form a hybrid electrical-chemical circuit (Knogler et al., 2014). Blocking glutamatergic transmission precludes double coils while sparing single coils (Knogler et al., 2014). In contrast, blocking glycinergic synapses led to triple or even quadruple coils (Knogler et al., 2014). These experimental observations suggest that synaptic excitation is required for successive coils after a first coil. Glycinergic transmission seems to prevent the generation of more than two successive coils. Patch-clamp recordings of MNs at this developmental stage exhibit the same isolated PDs and SBs from earlier developmental stages and show mixed events in which a PD event immediately follows an SB or vice-versa (Knogler et al., 2014). Interestingly, the application of CNQX eliminates mixed PD-SB events but not single isolated SBs, suggesting that the coupling of PD and SB in mixed events is glutamatergic (Knogler et al., 2014). Therefore, we aimed to generate a model with the following characteristics: (1) double coils lasting about 1 s in duration accounting for over half of the coiling events, (2) a dependence of double coiling upon excitatory synaptic transmission, (3) an increase in multiple coiling events in the absence of inhibitory synaptic transmission, and (4) the presence of mixed PD-SB events with similar sensitivity to the blockade of excitatory synaptic transmission as double coils.

Network description

(Figure 3A) To implement a model capable of generating double coils that depend upon glutamatergic transmission, we built upon the single coiling model by adding two populations of neurons. We reasoned that if double coiling depended upon excitatory neurotransmission, then a population of commissural excitatory neurons could be necessary to trigger a second contralateral coil in double coils. V0v neurons are a population of glutamatergic commissural INs, some of which may be involved in larger amplitude locomotor movements such as coiling (Jay and McLean, 2019). Thus, we added a chain of V0vs (10 neurons for each side) electrically coupled to the previous scaffold (i.e., the ipsilateral IC-MN-V0d scaffold). To generate the crossing excitation underlying the second coil, all V0vs projected glutamatergic synapses to contralateral ICs. Electrical synapses were formed with neighboring MNs, V0ds, and V0vs (the nearest three of each type of neuron in both the rostral and the caudal directions). Ipsilateral ICs were coupled with V0vs in the first four rostral somites.

Figure 3 with 7 supplements see all

Download asset Open asset

The second population of neurons that we added was ipsilaterally projecting excitatory neurons present at this stage and shown to receive mixed PD-SB events (Knogler et al., 2014). These neurons have been suggested to be circumferential ipsilateral descending neurons that arise from the V2a population. In the model, V2as were electrically coupled to IC neurons and projected glutamatergic synapses to V0vs. This chemical synapse caused a delay after the initiation of the initial coil that facilitates the second contralateral coil (see below). V2as most likely also excite MNs, based on the data from Knogler et al., 2014. For computational simplicity, we omitted this connection as it was unnecessary for double coilings to be generated, though this may reduce the amplitude of the coils. As V2as display SBs at this stage, we modeled glycinergic projections from V0ds to contralateral V2as such that the ith V0d projected to all contralateral V2as between the i−5 and i+5 segments like how V0ds project to contralateral MNs.

Left and right ICs received a tonic motor command though we delayed the activation of the tonic command to the right side by 1500 ms to ensure that double coils were not near-coincident bilateral single coils. We modified several parameters of the ICs, most notably increases in the a and the k parameters, to produce a more extended inter-coiling period (Figure 3B) than seen in single coiling (Knogler et al., 2014). A reminder that the a parameter represents the time scale of the recovery variable u that returns the membrane potential to rest. The k parameter shapes subthreshold dynamics.

Simulations results 

Simulations of the double coiling model frequently generated pairs of successive, left-right alternating coils lasting about 1 s in total (Figure 3C, Figure 3—video 1). In five 100,000-ms long runs of the base model with a minimal amount of variability added to several model parameters (${\sigma }_d$ = 0.5, ${\sigma }_p$ = 0.01, and ${\sigma }_w$ = 0.05), approximately 60% of events were double coils, with the rest mainly being single coils (31%), and very few triple coils or truncated single coils (5% each) (Figure 3I).

The timing of ICs, V2as, and V0vs (Figure 3C) suggests that double coils were generated by ipsilateral recruitment of V2as and V0vs during the first coil, which led to activation of the contralateral ICs to initiate the second coil. This sequence is supported by an analysis of the phase delays (Figure 3D). IC firing precedes the firing of all other ipsilateral spinal neurons suggesting they drive the activity of each coil. V2a activity precedes that of V0vs, which suggests that V2as recruit V0vs. This recruitment of V0vs by V2as is supported by simulations where the V2a to V0v synapse is removed (Figure 3—figure supplement 1A). V0v activity succeeds all other ipsilateral spinal INs, suggesting they are the last INs active during the first coil in a double coiling event. A key to generating double coils in our model was thus to delay the activation of V0vs. This delay enabled the activation of contralateral ICs after the first coiling is completed and when commissural inhibition of the contralateral IC has also terminated. If the activation of ipsilateral V0vs occurred too early during the first coiling, which can be produced by increasing the weight of the V2a to V0v and the V0v to contralateral IC glutamatergic synapses, the occurrence of a second coil is less likely (Figure 3E,I, Figure 3—video 2).

To further underscore the importance of glutamatergic transmission to double coiling as reported experimentally (Knogler et al., 2014), blocking glutamatergic transmission in the model greatly reduced the number of double coils (Figure 3F,I, Figure 3—video 3). On the other hand, blocking glycinergic synapses increased multiple coilings of three or more coils (Figure 3G,I, Figure 3—video 4) as Knogler et al., 2014 reported. This effect presumably occurs due to the unopposed reverberating commissural excitation of ICs by V0vs. Indeed, silencing V0v synapses in a model with no glycinergic synapses blocks double and multiple coils (Figure 3—figure supplement 1B).

The sequencing of commissural excitation and inhibition in the generation of double coils is further underscored by the presence of mixed SB-PD or PD-SB events (Figure 3H) observed experimentally in hyperpolarized MNs (Knogler et al., 2014). In these mixed events, the PDs were generated by gap junction coupled ICs during the ipsilateral coil, whereas contralateral V0ds activated during the contralateral coil generated the SBs in the ipsilateral MNs. Blocking glutamatergic transmission in our model uncoupled PDs and SBs (Figure 3H) as observed experimentally (Knogler et al., 2014).

Just as the robustness of the single coiling model was tested through modifications to the base model and several sensitivity tests, we also tested the robustness of the double coiling model. First, we increased the size of the model from 10 to 30 somites. Modifications of the tonic motor command amplitude, length of IC axons, gap junction coupling from IC to MN, and the synapses from MN to muscle cell, V0v to IC, V2a to V0v enabled the generation of double coils in this model (Figure 3—figure supplement 1C, Figure 3—video 5). In the 10-somite base model, we also tested the role of the glycinergic reversal potential. The value of this parameter was hyperpolarized from −45 mV in the single coiling model to −58 mV in the base double coiling model. This shift was intended to reflect gradual hyperpolarization of the reversal potential of glycine during development (Ben-Ari, 2002; Saint-Amant and Drapeau, 2000, Saint-Amant and Drapeau, 2001). We tested the double coiling model at values ranging between −46 and −70 mV (Figure 3—figure supplement 2A). We found that the proportion of double coils seemed to be higher, and the proportion of multiple coils was increased at more depolarized values of the glycinergic reversal potential. The proportion of double coils was relatively constant at more hyperpolarized values of the glycinergic reversal potential.

To test whether the double coiling model was sensitive to within-model parameter variability, we ran sets of ten 100-s long simulations at various values of ${\sigma }_d$, ${\sigma }_l$, ${\sigma }_p$, and ${\sigma }_w$ (Figure 3J-M). Again, we found that relatively small levels of variability in the parameters governing membrane dynamics (${\sigma }_p$) decreased the proportion of coiling events that were double coils and increased the number of truncated coils (Figure 3L). Moderate levels of variability in the parameters governing axonal length (${\sigma }_l$) or synaptic weight (${\sigma }_w$) decreased the proportion of double coils while increasing single coils and sometimes truncated coils (Figure 3K, M). The proportion of coiling events was largely unaffected by variability in the amplitude of the motor command (${\sigma }_d$, Figure 3J).

Considering that the generation of double coils was sensitive to chemical synaptic activity and gap junctions (Knogler et al., 2014), we tested the sensitivity of the model to variability in the weights of only chemical synapses (${\sigma }_{w,chem}$) and only gap junctions (${\sigma }_{w,gap}$). We found that the proportion of double coils was relatively more sensitive to the variability of gap junctions than chemical synapses (Figure 3—figure supplement 2B, C).

Generation of swimming pattern by spinal network oscillators (>2–3 dpf)

Around 2 or 3 dpf, zebrafish transition from coiling movements to swimming (Drapeau et al., 2002; Saint-Amant and Drapeau, 1998). This transition entails two fundamental changes in locomotor movements: first, long, slow coils are replaced by quick, short tail beats; and second, the number of consecutive tail beats are increased from the two side-to-side coilings seen in double coils to multiple consecutive tail beats that compose each swimming episode. One of the emerging swimming movements is beat-and-glide swimming, characterized by short swimming episodes lasting several hundreds of milliseconds separated by gliding pauses and lasting several hundreds of milliseconds (Budick and O'Malley, 2000; Buss and Drapeau, 2001). Swim episodes consist of repetitive left-right alternating, low-amplitude tail beats that propagate from the rostral toward the caudal end of the fish body and are generated approximately at 20–65 Hz (Budick and O'Malley, 2000; Buss and Drapeau, 2001).

Beat-and-glide swimming can be produced in isolated larval zebrafish spinal cord preparations by NMDA application (Lambert et al., 2012; McDearmid and Drapeau, 2006; Wiggin et al., 2012) or by optogenetic stimulation of excitatory spinal neurons (Wahlstrom-Helgren et al., 2019). This capacity suggests that the transition from coiling to swimming involves a delegation of rhythm generation from ICs to spinal locomotor circuits (Roussel et al., 2020). Therefore, we sought to model a spinal network that could generate beat-and-glide swimming activity hallmarks—swim episodes lasting about 200–300 ms with repeated left-right alternating low-amplitude tail beats at around 20–65 Hz—without relying on pacemaker cells.

Recent experimental studies have also started to delineate the contributions of specific populations of spinal neurons to swimming. Ablation of ipsilaterally projecting, excitatory neurons in the V2a population eliminates swimming activity (Eklöf-Ljunggren et al., 2012). Genetic ablation of ipsilaterally projecting, inhibitory neurons in the V1 population affects swim vigor but has no effects on the patterning of swimming (Kimura and Higashijima, 2019). Genetic ablation of a subset of commissural inhibitory neurons in the dI6 population reduces left-right alternation (Satou et al., 2020). We sought to replicate the role of these neurons in our model.

Network description 

(Figure 4A) Whereas coiling is likely to be generated by primary MNs, swimming is more likely to involve secondary MNs (Ampatzis et al., 2013; Liu and Westerfield, 1988). There are more secondary than primary MNs, and new spinal neurons are born at the same time as secondary MNs (Bernhardt et al., 1990). To emulate the increase in the number of spinal neurons that may underlie swimming, we increased the size of the fish from 10 to 15 segments and accordingly increased the number of MNs, V0vs, and V2as from 10 to 15. Thus, each model somite in our swimming model represented two biological somites instead of three in our coiling models. IC neurons were removed from the model to reduce computational load. We are not aware of any experimental evidence of the involvement of IC neurons in later swimming stages.

Figure 4 with 1 supplement see all

Download asset Open asset

We introduced two populations of neurons for the beat-and-glide model. Commissural inhibitory CoBL neurons (including neurons from the dI6 and V0d populations) are active during swimming (Liao and Fetcho, 2008; Satou et al., 2020). V0d neurons were involved with faster swimming, whereas dI6s were more likely to be active during slower swimming (Satou et al., 2020). Therefore, we modeled CoBL neurons as dI6s. dI6s thus replaced the V0ds in the coiling models as the source of contralateral inhibition in the swimming model. CoBL neurons have been shown to project to MNs, dI6s, and unidentified ipsilateral descending spinal neurons that could be V2as (Satou et al., 2020). We added fifteen dI6s per side. We modeled the projection pattern of dI6s based on the bifurcating axons with short ascending and long descending branches of a subset of neurons in the dI6 subpopulation (Satou et al., 2020). Thus, the ith dI6s projected ascending branches to their rostral targets in the i−1th segment and projected descending branches to their caudal targets between the i+1 and i+3 segments.

A second new population of neurons was the V1 INs that include circumferential ascending INs that emerge during the second wave of neurogenesis in the spinal cord (Bernhardt et al., 1990). While V1 neurons were not included in the coiling models because their role in that form of movement is unclear, experiments in which the genetically identified V1 neurons are ablated suggest a role in controlling swim vigor (Kimura and Higashijima, 2019). We thus modeled V1s as a population of tonic firing neurons with ipsilateral ascending glycinergic projections. We added 15 V1s per side. We distributed V1s from segment 2 to the caudal end because our preliminary simulations suggested that starting the distribution of V1s at segment 1 made the episode duration more variable. In our model, V1s project segmental and ascending ipsilateral glycinergic synapses with rostral V2as (Kimura and Higashijima, 2019), dI6s, and V0vs (Sengupta et al., 2021) such that the V1s in the ith segment project to their rostral targets in the i−1 to i−2 segments. The V1 projections were short based upon recent evidence that their projections to motor circuits are constrained to segmental and immediately neighboring somites (Sengupta et al., 2021). Reciprocally, V2as formed glutamatergic synapses to caudally located V1. The reversal potential of the glycinergic synapses from dI6 and V1 neurons was set at −70 mV. This value is hyperpolarized compared to values in the coiling model considering the known hyperpolarization of the chloride equilibrium potential during development (Ben-Ari, 2002), which has also been observed in the zebrafish nervous system (Zhang et al., 2010).

V2as were considered the primary source of rhythmogenesis in our models of beat-and-glide swimming based on previous studies showing the necessity and sufficiency of V2a neurons to swimming activity (Eklöf-Ljunggren et al., 2012; Ljunggren et al., 2014). In the swimming model, V2as projected segmental and descending projections to dI6s, MNs, V0vs, V1s, and caudal V2as (the ith V2a projected to all caudal V0vs, dI6s, V2as, and MNs between the i+1 and i+6 segments). Connections between V2as and V2as (Ampatzis et al., 2014; Menelaou and McLean, 2019; Song et al., 2020) and from V2a neurons to MNs (Ampatzis et al., 2014; Menelaou and McLean, 2019; Song et al., 2020) have been reported. Connections from V2as to dI6s have not been studied yet, and so we based them on the reported connections from V2a neurons to V0d neurons, another population of commissural inhibitory INs (Menelaou and McLean, 2019). A subtype of V2a neurons that project to MNs was shown to bifurcate and have short ascending branches (Menelaou and McLean, 2019), which we modeled in addition to an ascending V2a to V0v connection (the ith V2a projected to rostral V0vs and MNs in the i−1 and i−2 segments) that remains to be confirmed.

Less is known about the connection patterns of commissural excitatory neurons at larval stages. However, in adult zebrafish, V0v commissural excitatory neurons have been shown to have bifurcating axons with shorter ascending branches and longer descending branches (Björnfors and El Manira, 2016). We modeled V0vs to project only to contralateral V2as. The ith V0vs projected ascending branches to their rostral targets in the i−1th segment and projected descending branches to their caudal targets between the i+1 and i+3 segments.

Whether rhythmic motor commands from supraspinal commands generate rhythmic tail beats at the spinal cord level is unclear. There is evidence for rhythmic and tonic activity in reticulospinal neurons involved in swimming (Kimura et al., 2013). However, the isolated zebrafish spinal cord can generate rhythmic activity similar to swimming (McDearmid and Drapeau, 2006; Wahlstrom-Helgren et al., 2019; Wiggin et al., 2012; Wiggin et al., 2014). Therefore, in our model, V2as received a tonic motor command in the form of a DC current to test whether the rhythm and pattern of swimming could be generated solely from the activity of spinal circuits.

Most spinal neurons at larval stages exhibit either tonic firing or firing with spike rate adaptation (Kimura and Higashijima, 2019; Menelaou and McLean, 2012; Menelaou and McLean, 2019; Satou et al., 2020), while a subset of MNs showing intrinsic burst firing (Menelaou and McLean, 2012). Therefore, we posited that the generation of rhythmic tail beats was not dependent upon the presence of intrinsically bursting neurons. Almost all of the neurons in the beat-and-glide swimming model were modeled to fire tonically (Figure 4B). Our base model was able to generate the beat-and-glide swimming pattern—alternating episodes of tail beats followed by silent inter-episode intervals each lasting hundreds of seconds—if V0vs were modeled to exhibit a more chattering or bursting firing pattern (Figure 4B).

As the model is symmetrical, including the motor command, we found that the model produces synchronous left-right activity unless we introduced some variability in commissural connections. With no a priori knowledge of where such variability could arise from, we chose to introduce a small amount of variability in the contralateral inhibition of dI6s by dI6s. The synaptic weights of this connection were scaled by a random number picked from a Gaussian distribution with mean =1, and standard deviation of 0.1, and this was sufficient to generate alternating left-right alternation. We did not seek to further characterize the variability required to generate left-right alternation.

Simulation results 

The beat-and-glide swimming model exhibited short-duration (hundreds of milliseconds) swimming episodes with left-right alternation and tail beat frequencies between 20 and 60 Hz (Figure 4C–G, Figure 4—video 1). The characteristics of the swimming episodes in our simulations were close to those described for free swimming in larval zebrafish by Buss and Drapeau, 2001, though the swimming output in our simulations had larger episode durations, shorter inter-episode intervals, and lower tail beat frequencies (Table 5).

Table 5

Parameter	Base model	Buss and Drapeau, 2001
Mean episode duration (ms)	234±6	180±20
Mean inter-episode interval (ms)	242±20	390±30
Mean tail beat frequency (Hz)	30.0±0.6	35±2

1. Values in mean± standard error. n=ten 10,000 ms-long simulations for the base model, and n=12 animals for the data from Buss and Drapeau, 2001.

An analysis of the phase delay between neuron populations during beat-and-glide swimming in the base model shows that the activity of ipsilateral, glutamatergic V2as precedes the activity of all ipsilateral neurons (Figure 4H). This earlier firing of V2as suggests that these spinal neurons drive the activity of the ipsilateral spinal swimming circuit. On the other hand, the glycinergic V1 neurons succeed all ipsilateral spinal neurons, suggesting they provide negative feedback to ipsilateral spinal swimming circuits. A hyperpolarized reversal potential of glycinergic synapses at this developmental stage strengthens the negative feedback (see Sensitivity testing to different values of the glycinergic reversal potential below). The contralateral dI6s and V0vs are out-of-phase with contralateral spinal neurons, consistent with their role in mediating left-right coordination. The longest delay between V2as and their ipsilateral counterparts was with the V0vs, which is reminiscent of V2a firing preceding V0v firing in the double coiling model to ensure a sufficient delay for the initiation of the second coil. Thus, the generation of alternating left-right tail beats would seem to require a certain delay in the excitation of contralateral swimming circuits.

To further investigate the role of specific neurons in the model’s swimming activity, we performed simulations composed of three 5000 ms long epochs: epoch 1, where the model was intact; epoch 2, where we silenced the targeted neurons by removing their synaptic inputs; and epoch 3, where the synaptic inputs to the targeted neurons were restored. Silencing V2as abolished the generation of tail beats (Figure 5A–F, Figure 5—figure supplement 1A–F), underscoring their primacy to the generation of tail beats (Eklöf-Ljunggren et al., 2012; Ljunggren et al., 2014). Commissural excitation mediated by V0vs seemed to be very important in maintaining the beat-and-glide pattern. Silencing V0vs diminished but did not eliminate the rhythmic firing of V2as or MNs. During epoch 2, the tonic motor command continued to activate V2as. Pairs of left-right tail beats may result from the commissural inhibition by dI6s that was still present. However, removing the contralateral excitation by V0v prevented the repetitive activation of the silent side after each tail beat, which severely reduced episode duration and the number of tail beats generated in each episode (Figure 5G–L, Figure 5—figure supplement 1G–L, Figure 5—video 1). These simulations suggest that commissural excitation is necessary to repeatedly activate the silent contralateral side during ongoing swimming to ensure successive left-right alternating tail beats and longer swim episodes (Björnfors and El Manira, 2016; Saint-Amant, 2010).

Figure 5 with 3 supplements see all

Download asset Open asset

Figure 5—source data 1

P-values for V2a null and V0v null simulations.

https://cdn.elifesciences.org/articles/67453/elife-67453-fig5-data1-v2.xlsx

Download elife-67453-fig5-data1-v2.xlsx

Genetic ablation of the ipsilaterally ascending inhibitory V1 INs increased swim vigor but produced no overt changes in swimming patterns (Kimura and Higashijima, 2019). Consistent with those experimental results, simulating the removal of ipsilateral ascending inhibition by silencing V1 INs seemed to increase the amplitude of MN activity (Figure 6A,B, Figure 6—figure supplement 1A–F, Figure 6—video 1). While the overall beat-and-glide pattern persisted, the duration of episodes was increased, and inter-episode intervals were shortened between epochs 1 and 2 (Figure 6C,D). Tail beat frequency was increased (Figure 6E). Left-right alternation was reduced in these simulations and during simulations where dI6s were silenced (Figure 6F,G,L, Figure 6—figure supplement 1G–L). The reduction in left-right alternation during simulations where dI6s were silenced was greater in caudal somites than in rostral somites (Figure 6L). Note that while left-right alternation was reduced, this did not preclude left-right alternating tail beats from being generated (Figure 6—video 2). The reduction of left-right coordination seen here was comparable to levels seen after genetic ablation of a commissural inhibitory subpopulation of dI6 INs (Satou et al., 2020) but is not sufficient to prevent left-right alternation. Since swimming is generated by rostrocaudal propagation of contractile waves, any left-right alternation in rostral segments will inevitably sway the rest of the body, as suggested by the musculoskeletal model. The precise kinematics of the tail beats will be affected by the reduction of left-right alternation observed (Figure 6—figure supplement 2, Figure 6—video 2). Finally, silencing dI6s had negligible effects on the episode duration and inter-episode interval but increased tail beat frequency when comparing epochs 1 and 2 and may increase the amplitude of motor activity (Figure 6H–K).

Figure 6 with 4 supplements see all

Download asset Open asset

Figure 6—source data 1

P-values for V1 null and dI6 null simulations.

https://cdn.elifesciences.org/articles/67453/elife-67453-fig6-data1-v2.xlsx

Download elife-67453-fig6-data1-v2.xlsx

Previous experimental results showed that strychnine application disturbed swimming in the 20–40 Hz range at later stages of development (Roussel et al., 2020). Our model’s behavior to loss of glycinergic transmission was tested. Removal of all glycinergic transmission led to continual tail beats with minimal gliding periods (Figure 7). Motor output was increased during the epoch of no glycinergic transmission (Figure 7B), episode duration was considerably lengthened, and the inter-episode interval was shortened (Figure 7C–D). The frequency of tail beats increased (Figure 7E). Left-right alternation was reduced, particularly at caudal somites (Figure 7F), which did not preclude left-right tail beats but altered swimming kinematics (Figure 7—figure supplement 1, Figure 7—video 1). These results indicate that removing glycinergic transmission in the model led to near-continuous swimming activity with altered kinematics and greater frequencies of tail beats.

Figure 7 with 2 supplements see all

Download asset Open asset

Figure 7—source data 1

P-values for strychnine simulations.

https://cdn.elifesciences.org/articles/67453/elife-67453-fig7-data1-v2.xlsx

Download elife-67453-fig7-data1-v2.xlsx

We then proceeded to test the sensitivity of the base model to several model parameters. Since some V2a INs in adult zebrafish have pacemaker capacities (Song et al., 2020), we tested whether we could also generate beat-and-glide swimming in a model with bursting V2a (Figure 8A). A model with bursting V2as where we also decreased the strength of synapses from V2as to other neurons and increased the connection strength of V0vs and dI6s to contralateral V2as (Table 3) generated 100–400 ms swimming episodes of left-right alternating tail beats at frequencies around 20–80 Hz interspersed by 100–400 ms long silent inter-episode intervals (Figure 8B–F, Figure 8—video 1). Surprisingly, a model with only tonic firing neurons (Figure 8G) was able to generate the hallmarks of beat-and-glide swimming as well (Figure 8H–L, Figure 8—video 2). While there were eventually longer episodes with shorter inter-episode intervals after 6000 ms, this simulation suggests that the architecture of the network is sufficient to generate beat-and-glide swimming for long periods despite the absence of any neurons with bursting properties.

Figure 8 with 5 supplements see all

Download asset Open asset

Figure 8—source data 1

P-values for sensitivity testing to values of E_gly.

https://cdn.elifesciences.org/articles/67453/elife-67453-fig8-data1-v2.xlsx

Download elife-67453-fig8-data1-v2.xlsx

We could generate beat-and-glide swimming in a 30-somite and a 10-somite model (Figure 8—figure supplement 1, Figure 8—video 3). We also examined the effects of changing the glycinergic reversal potential that was set at −70 mV for the base beat-and-glide model (Figure 8—figure supplement 2). We tested the model with E_gly values between −56 mV, near the value of E_gly in the double coiling model, and −72 mV (Figure 8—figure supplement 2A,B). The lower bound of −72 mV is close to the value of E_GABA reported in zebrafish retinal ganglion cells around 3–4 dpf (Zhang et al., 2010). Episode duration was increased, and inter-episode intervals decreased at more depolarized values of the glycinergic reversal potential leading to the loss of the beat-and-glide pattern (Figure 8—figure supplement 2C). Tail beat frequency also increased as glycinergic currents decreased at depolarized glycinergic reversal potentials, and left-right alternation was replaced by left-right synchrony (Figure 8—figure supplement 2D).

Finally, the sensitivity of the base model to variability was tested by running sets of ten 10,000-ms long simulations at various values of ${\sigma }_d$, ${\sigma }_l$, ${\sigma }_p$, and ${\sigma }_w$. We also performed ten 10,000-ms long simulations of the base model (a reminder that there is a random scaling factor to the dI6 to contralateral dI6 synapse in the base model). The episode duration, inter-episode interval, average tail beat frequency in each episode, and the minimum coefficient of the cross-correlation of the left and right muscle output were analyzed (Figures 9 and 10). Increases in variability in the motor command drive (${\sigma }_d$) seemed to affect inter-episode intervals and left-right alternation but not episode duration and tail beat frequency (Figure 9A-D). At similar levels of variability in rostrocaudal axonal length (${\sigma }_l$), all four measures of swimming activity were perturbed (Figure 9E-H). On the other hand, variability in synaptic weights (${\sigma }_w)$ affected only episode duration and inter-episode duration (Figure 9I-L).

Figure 9

Download asset Open asset

Figure 9—source data 1

P-values for sensitivity testing to values of sigma D, L and W.

https://cdn.elifesciences.org/articles/67453/elife-67453-fig9-data1-v2.xlsx

Download elife-67453-fig9-data1-v2.xlsx

Figure 10

Download asset Open asset

Figure 10—source data 1

P-values for sensitivity testing to values of Sigma P.

https://cdn.elifesciences.org/articles/67453/elife-67453-fig10-data1-v2.xlsx

Download elife-67453-fig10-data1-v2.xlsx

Smaller variability in the parameters shaping the membrane potential dynamics (${\sigma }_p)$ disrupted the beat-and-glide pattern with no beat-and-glide swimming observed at some values of ${\sigma }_p$ (Figure 10A, B). As ${\sigma }_p$ was increased, the episode duration, inter-episode interval, and tail beat frequency were disrupted (Figure 10C, D). Even slight variability in the parameters shaping membrane potential dynamics (e.g., ${\sigma }_p$=0.01) resulted in changes in membrane excitability and, in some cases, firing patterns as evidenced by the conversion of some V0vs from burst to tonic firing (Figure 10E). Thus, the beat-and-glide model was most susceptible to variations in the parameters determining membrane potential dynamics and similarly sensitive to the other parameters tested.

To our knowledge, this study presents some of the first models of spinal locomotor circuits in developing zebrafish. We have built several spinal locomotor circuit models that generate locomotor movements of the developing zebrafish (Figure 11, Figure 11—video 1). These models support mechanisms of network operation of developing zebrafish spinal locomotor circuits described experimentally. Our models suggest that the circuitry driving locomotor movements could switch from a pacemaker kernel located rostrally during coiling maneuvers to network-based spinal circuits during swimming. Results from simulations where populations of spinal neurons are silenced were consistent with experimental studies. Our sensitivity analysis suggests that the correct operation of spinal circuits for locomotion is not immune to variations in firing behaviors, length of axonal projections, motor command amplitude, and synaptic weighting. However, the sensitivity to these parameters is variable.

Figure 11 with 1 supplement see all

Download asset Open asset

The code for the models and for the figures, as well as the data used to make the figures, can be accessed at https://github.com/bui-lab/code (copy archived at https://archive.softwareheritage.org/swh:1:rev:3268f1f684937f619b6ad87cd0297f8c7ea66db0). Updates and revisions to the models will also be made available at this site.

1. 1. Ampatzis K
   2. Song J
   3. Ausborn J
   4. El Manira A

   (2013) Pattern of innervation and recruitment of different classes of motoneurons in adult zebrafish

   Journal of Neuroscience 33:10875–10886.

   https://doi.org/10.1523/JNEUROSCI.0896-13.2013

   • PubMed
   • Google Scholar
2. 1. Ampatzis K
   2. Song J
   3. Ausborn J
   4. El Manira A

   (2014) Separate microcircuit modules of distinct v2a interneurons and motoneurons control the speed of locomotion

   Neuron 83:934–943.

   https://doi.org/10.1016/j.neuron.2014.07.018

   • PubMed
   • Google Scholar
3. 1. Anderson TM
   2. Abbinanti MD
   3. Peck JH
   4. Gilmour M
   5. Brownstone RM
   6. Masino MA

   (2012) Low-threshold calcium currents contribute to locomotor-like activity in neonatal mice

   Journal of Neurophysiology 107:103–113.

   https://doi.org/10.1152/jn.00583.2011

   • PubMed
   • Google Scholar
4. 1. Ausborn J
   2. Shevtsova NA
   3. Caggiano V
   4. Danner SM
   5. Rybak IA

   (2019) Computational modeling of brainstem circuits controlling locomotor frequency and gait

   eLife 8:e43587.

   https://doi.org/10.7554/eLife.43587

   • PubMed
   • Google Scholar
5. 1. Bagnall MW
   2. McLean DL

   (2014) Modular organization of axial microcircuits in zebrafish

   Science 343:197–200.

   https://doi.org/10.1126/science.1245629

   • PubMed
   • Google Scholar
6. 1. Ben-Ari Y

   (2002) Excitatory actions of gaba during development: the nature of the nurture

   Nature Reviews Neuroscience 3:728–739.

   https://doi.org/10.1038/nrn920

   • PubMed
   • Google Scholar
7. 1. Bennett MV
   2. Zukin RS

   (2004) Electrical coupling and neuronal synchronization in the mammalian brain

   Neuron 41:495–511.

   https://doi.org/10.1016/S0896-6273(04)00043-1

   • PubMed
   • Google Scholar
8. 1. Bernhardt RR
   2. Chitnis AB
   3. Lindamer L
   4. Kuwada JY

   (1990) Identification of spinal neurons in the embryonic and larval zebrafish

   The Journal of Comparative Neurology 302:603–616.

   https://doi.org/10.1002/cne.903020315

   • Google Scholar
9. 1. Bicanski A
   2. Ryczko D
   3. Cabelguen JM
   4. Ijspeert AJ

   (2013) From lamprey to salamander: an exploratory modeling study on the architecture of the spinal locomotor networks in the salamander

   Biological Cybernetics 107:565–587.

   https://doi.org/10.1007/s00422-012-0538-y

   • PubMed
   • Google Scholar
10. 1. Björnfors ER
    2. El Manira A

    (2016) Functional diversity of excitatory commissural interneurons in adult zebrafish

    eLife 5:e18579.

    https://doi.org/10.7554/eLife.18579

    • PubMed
    • Google Scholar
11. 1. Böhm UL
    2. Prendergast A
    3. Djenoune L
    4. Nunes Figueiredo S
    5. Gomez J
    6. Stokes C
    7. Kaiser S
    8. Suster M
    9. Kawakami K
    10. Charpentier M
    11. Concordet J-P
    12. Rio J-P
    13. Del Bene F
    14. Wyart C

    (2016) CSF-contacting neurons regulate locomotion by relaying mechanical stimuli to spinal circuits

    Nature Communications 7:1–8.

    https://doi.org/10.1038/ncomms10866

    • Google Scholar
12. 1. Brustein E
    2. Chong M
    3. Holmqvist B
    4. Drapeau P

    (2003) Serotonin patterns locomotor network activity in the developing zebrafish by modulating quiescent periods

    Journal of Neurobiology 57:303–322.

    https://doi.org/10.1002/neu.10292

    • PubMed
    • Google Scholar
13. 1. Brustein E
    2. Drapeau P

    (2005) Serotoninergic modulation of chloride homeostasis during maturation of the locomotor network in zebrafish

    Journal of Neuroscience 25:10607–10616.

    https://doi.org/10.1523/JNEUROSCI.2017-05.2005

    • PubMed
    • Google Scholar
14. 1. Budick SA
    2. O'Malley DM

    (2000) Locomotor repertoire of the larval zebrafish: swimming, turning and prey capture

    Journal of Experimental Biology 203:2565–2579.

    https://doi.org/10.1242/jeb.203.17.2565

    • Google Scholar
15. 1. Buss RR
    2. Drapeau P

    (2000) Physiological properties of zebrafish embryonic red and white muscle fibers during early development

    Journal of Neurophysiology 84:1545–1557.

    https://doi.org/10.1152/jn.2000.84.3.1545

    • PubMed
    • Google Scholar
16. 1. Buss RR
    2. Drapeau P

    (2001) Synaptic drive to motoneurons during fictive swimming in the developing zebrafish

    Journal of Neurophysiology 86:197–210.

    https://doi.org/10.1152/jn.2001.86.1.197

    • PubMed
    • Google Scholar
17. 1. Callahan RA
    2. Roberts R
    3. Sengupta M
    4. Kimura Y
    5. Higashijima SI
    6. Bagnall MW

    (2019) Spinal V2b neurons reveal a role for ipsilateral inhibition in speed control

    eLife 8:e47837.

    https://doi.org/10.7554/eLife.47837

    • PubMed
    • Google Scholar
18. 1. Danner SM
    2. Zhang H
    3. Shevtsova NA
    4. Borowska-Fielding J
    5. Deska-Gauthier D
    6. Rybak IA
    7. Zhang Y

    (2019) Spinal V3 interneurons and Left-Right coordination in mammalian locomotion

    Frontiers in Cellular Neuroscience 13:516.

    https://doi.org/10.3389/fncel.2019.00516

    • PubMed
    • Google Scholar
19. 1. Del Negro CA
    2. Hayes JA
    3. Pace RW
    4. Brush BR
    5. Teruyama R
    6. Feldman JL

    (2010) Synaptically activated burst-generating conductances may underlie a group-pacemaker mechanism for respiratory rhythm generation in mammals

    Progress in Brain Research 187:111–136.

    https://doi.org/10.1016/B978-0-444-53613-6.00008-3

    • PubMed
    • Google Scholar
20. 1. Domenici P
    2. Hale ME

    (2019) Escape responses of fish: a review of the diversity in motor control, kinematics and behaviour

    Journal of Experimental Biology 222:166009.

    https://doi.org/10.1242/jeb.166009

    • Google Scholar
21. 1. Drapeau P
    2. Ali DW
    3. Buss RR
    4. Saint-Amant L

    (1999) In vivo recording from identifiable neurons of the locomotor network in the developing zebrafish

    Journal of Neuroscience Methods 88:1–13.

    https://doi.org/10.1016/S0165-0270(99)00008-4

    • PubMed
    • Google Scholar
22. 1. Drapeau P
    2. Saint-Amant L
    3. Buss RR
    4. Chong M
    5. McDearmid JR
    6. Brustein E

    (2002) Development of the locomotor network in zebrafish

    Progress in Neurobiology 68:85–111.

    https://doi.org/10.1016/S0301-0082(02)00075-8

    • PubMed
    • Google Scholar
23. 1. Eklöf-Ljunggren E
    2. Haupt S
    3. Ausborn J
    4. Dehnisch I
    5. Uhlén P
    6. Higashijima S
    7. El Manira A

    (2012) Origin of excitation underlying locomotion in the spinal circuit of zebrafish

    PNAS 109:5511–5516.

    https://doi.org/10.1073/pnas.1115377109

    • PubMed
    • Google Scholar
24. 1. England S
    2. Batista MF
    3. Mich JK
    4. Chen JK
    5. Lewis KE

    (2011) Roles of hedgehog pathway components and retinoic acid signalling in specifying zebrafish ventral spinal cord neurons

    Development 138:5121–5134.

    https://doi.org/10.1242/dev.066159

    • PubMed
    • Google Scholar
25. 1. Favero M
    2. Cangiano A
    3. Busetto G

    (2014) The timing of activity is a regulatory signal during development of neural connections

    Journal of Molecular Neuroscience 53:324–329.

    https://doi.org/10.1007/s12031-013-0128-z

    • PubMed
    • Google Scholar
26. 1. Ferrario A
    2. Merrison-Hort R
    3. Soffe SR
    4. Borisyuk R

    (2018) Structural and functional properties of a probabilistic model of neuronal connectivity in a simple locomotor network

    eLife 7:e33281.

    https://doi.org/10.7554/eLife.33281

    • PubMed
    • Google Scholar
27. 1. Gabriel JP
    2. Mahmood R
    3. Kyriakatos A
    4. Söll I
    5. Hauptmann G
    6. Calabrese RL
    7. El Manira A

    (2009) Serotonergic modulation of locomotion in zebrafish: endogenous release and synaptic mechanisms

    Journal of Neuroscience 29:10387–10395.

    https://doi.org/10.1523/JNEUROSCI.1978-09.2009

    • PubMed
    • Google Scholar
28. 1. Goldman MS
    2. Golowasch J
    3. Marder E
    4. Abbott LF

    (2001) Global structure, robustness, and modulation of neuronal models

    The Journal of Neuroscience 21:5229–5238.

    https://doi.org/10.1523/JNEUROSCI.21-14-05229.2001

    • PubMed
    • Google Scholar
29. 1. Golowasch J
    2. Marder E

    (1992) Ionic currents of the lateral pyloric neuron of the stomatogastric ganglion of the crab

    Journal of Neurophysiology 67:318–331.

    https://doi.org/10.1152/jn.1992.67.2.318

    • PubMed
    • Google Scholar
30. 1. Higashijima S-I
    2. Mandel G
    3. Fetcho JR

    (2004) Distribution of prospective glutamatergic, Glycinergic, and GABAergic neurons in embryonic and larval zebrafish

    The Journal of Comparative Neurology 480:1–18.

    https://doi.org/10.1002/cne.20278

    • Google Scholar
31. 1. Hull MJ
    2. Soffe SR
    3. Willshaw DJ
    4. Roberts A

    (2016) Modelling feedback excitation, pacemaker properties and sensory switching of electrically coupled brainstem neurons controlling rhythmic activity

    PLOS Computational Biology 12:e1004702.

    https://doi.org/10.1371/journal.pcbi.1004702

    • PubMed
    • Google Scholar
32. 1. Ijspeert AJ
    2. Crespi A
    3. Ryczko D
    4. Cabelguen JM

    (2007) From swimming to walking with a salamander robot driven by a spinal cord model

    Science 315:1416–1420.

    https://doi.org/10.1126/science.1138353

    • PubMed
    • Google Scholar
33. Book

    1. Izhikevich EM

    (2007)

    Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting

    The MIT Press.

    • Google Scholar
34. 1. Jay M
    2. De Faveri F
    3. McDearmid JR

    (2015) Firing dynamics and modulatory actions of supraspinal dopaminergic neurons during zebrafish locomotor behavior

    Current Biology 25:435–444.

    https://doi.org/10.1016/j.cub.2014.12.033

    • PubMed
    • Google Scholar
35. 1. Jay M
    2. McLean DL

    (2019) Reconciling the functions of even-skipped interneurons during crawling, swimming, and walking

    Current Opinion in Physiology 8:188–192.

    https://doi.org/10.1016/j.cophys.2019.02.003

    • PubMed
    • Google Scholar
36. 1. Jha U
    2. Thirumalai V

    (2020) Neuromodulatory selection of motor neuron recruitment patterns in a visuomotor behavior increases speed

    Current Biology 30:788–801.

    https://doi.org/10.1016/j.cub.2019.12.064

    • PubMed
    • Google Scholar
37. 1. Kimmel CB
    2. Warga RM
    3. Kane DA

    (1994) Cell cycles and clonal strings during formation of the zebrafish central nervous system

    Development 120:265–276.

    https://doi.org/10.1242/dev.120.2.265

    • PubMed
    • Google Scholar
38. 1. Kimmel CB
    2. Ballard WW
    3. Kimmel SR
    4. Ullmann B
    5. Schilling TF

    (1995) Stages of embryonic development of the zebrafish

    Developmental Dynamics 203:253–310.

    https://doi.org/10.1002/aja.1002030302

    • PubMed
    • Google Scholar
39. 1. Kimura Y
    2. Satou C
    3. Fujioka S
    4. Shoji W
    5. Umeda K
    6. Ishizuka T
    7. Yawo H
    8. Higashijima S

    (2013) Hindbrain V2a neurons in the excitation of spinal locomotor circuits during zebrafish swimming

    Current Biology 23:843–849.

    https://doi.org/10.1016/j.cub.2013.03.066

    • PubMed
    • Google Scholar
40. 1. Kimura Y
    2. Higashijima SI

    (2019) Regulation of locomotor speed and selection of active sets of neurons by V1 neurons

    Nature Communications 10:2268.

    https://doi.org/10.1038/s41467-019-09871-x

    • PubMed
    • Google Scholar
41. 1. Knogler LD
    2. Ryan J
    3. Saint-Amant L
    4. Drapeau P

    (2014) A hybrid electrical/chemical circuit in the spinal cord generates a transient embryonic motor behavior

    Journal of Neuroscience 34:9644–9655.

    https://doi.org/10.1523/JNEUROSCI.1225-14.2014

    • PubMed
    • Google Scholar
42. 1. Kozlov A
    2. Huss M
    3. Lansner A
    4. Kotaleski JH
    5. Grillner S

    (2009) Simple cellular and network control principles govern complex patterns of motor behavior

    PNAS 106:20027–20032.

    https://doi.org/10.1073/pnas.0906722106

    • PubMed
    • Google Scholar
43. 1. Kozlov AK
    2. Kardamakis AA
    3. Hellgren Kotaleski J
    4. Grillner S

    (2014) Gating of steering signals through phasic modulation of reticulospinal neurons during locomotion

    PNAS 111:3591–3596.

    https://doi.org/10.1073/pnas.1401459111

    • PubMed
    • Google Scholar
44. 1. Lambert AM
    2. Bonkowsky JL
    3. Masino MA

    (2012) The conserved dopaminergic diencephalospinal tract mediates vertebrate locomotor development in zebrafish larvae

    Journal of Neuroscience 32:13488–13500.

    https://doi.org/10.1523/JNEUROSCI.1638-12.2012

    • PubMed
    • Google Scholar
45. 1. Liao JC
    2. Fetcho JR

    (2008) Shared versus specialized glycinergic spinal interneurons in axial motor circuits of larval zebrafish

    Journal of Neuroscience 28:12982–12992.

    https://doi.org/10.1523/JNEUROSCI.3330-08.2008

    • PubMed
    • Google Scholar
46. 1. Liu YC
    2. Hale ME

    (2017) Local spinal cord circuits and bilateral mauthner cell activity function together to drive alternative startle behaviors

    Current Biology 27:697–704.

    https://doi.org/10.1016/j.cub.2017.01.019

    • PubMed
    • Google Scholar
47. 1. Liu DW
    2. Westerfield M

    (1988) Function of identified motoneurones and co-ordination of primary and secondary motor systems during zebra fish swimming

    The Journal of Physiology 403:73–89.

    https://doi.org/10.1113/jphysiol.1988.sp017239

    • PubMed
    • Google Scholar
48. 1. Ljunggren EE
    2. Haupt S
    3. Ausborn J
    4. Ampatzis K
    5. El Manira A

    (2014) Optogenetic activation of excitatory premotor interneurons is sufficient to generate coordinated locomotor activity in larval zebrafish

    Journal of Neuroscience 34:134–139.

    https://doi.org/10.1523/JNEUROSCI.4087-13.2014

    • PubMed
    • Google Scholar
49. 1. Marder E
    2. Taylor AL

    (2011) Multiple models to capture the variability in biological neurons and networks

    Nature Neuroscience 14:133–138.

    https://doi.org/10.1038/nn.2735

    • PubMed
    • Google Scholar
50. 1. McDearmid JR
    2. Drapeau P

    (2006) Rhythmic motor activity evoked by NMDA in the spinal zebrafish larva

    Journal of Neurophysiology 95:401–417.

    https://doi.org/10.1152/jn.00844.2005

    • PubMed
    • Google Scholar
51. 1. McLean DL
    2. Fan J
    3. Higashijima S
    4. Hale ME
    5. Fetcho JR

    (2007) A topographic map of recruitment in spinal cord

    Nature 446:71–75.

    https://doi.org/10.1038/nature05588

    • PubMed
    • Google Scholar
52. 1. McLean DL
    2. Masino MA
    3. Koh IY
    4. Lindquist WB
    5. Fetcho JR

    (2008) Continuous shifts in the active set of spinal interneurons during changes in locomotor speed

    Nature Neuroscience 11:1419–1429.

    https://doi.org/10.1038/nn.2225

    • PubMed
    • Google Scholar
53. 1. McLean DL
    2. Fetcho JR

    (2009) Spinal interneurons differentiate sequentially from those driving the fastest swimming movements in larval zebrafish to those driving the slowest ones

    Journal of Neuroscience 29:13566–13577.

    https://doi.org/10.1523/JNEUROSCI.3277-09.2009

    • PubMed
    • Google Scholar
54. 1. Menelaou E
    2. McLean DL

    (2012) A gradient in endogenous rhythmicity and oscillatory drive matches recruitment order in an axial motor pool

    Journal of Neuroscience 32:10925–10939.

    https://doi.org/10.1523/JNEUROSCI.1809-12.2012

    • PubMed
    • Google Scholar
55. 1. Menelaou E
    2. McLean DL

    (2019) Hierarchical control of locomotion by distinct types of spinal V2a interneurons in zebrafish

    Nature Communications 10:4197.

    https://doi.org/10.1038/s41467-019-12240-3

    • PubMed
    • Google Scholar
56. 1. Messina JA
    2. St. Paul A
    3. Hargis S
    4. Thompson WE
    5. McClellan AD

    (2017) Elimination of Left-Right reciprocal coupling in the adult lamprey spinal cord abolishes the generation of locomotor activity

    Frontiers in Neural Circuits 11:00089.

    https://doi.org/10.3389/fncir.2017.00089

    • Google Scholar
57. 1. Montgomery JE
    2. Wahlstrom-Helgren S
    3. Wiggin TD
    4. Corwin BM
    5. Lillesaar C
    6. Masino MA

    (2018) Intraspinal serotonergic signaling suppresses locomotor activity in larval zebrafish

    Developmental Neurobiology 78:807–827.

    https://doi.org/10.1002/dneu.22606

    • Google Scholar
58. 1. Myers PZ
    2. Eisen JS
    3. Westerfield M

    (1986) Development and axonal outgrowth of identified motoneurons in the zebrafish

    The Journal of Neuroscience 6:2278–2289.

    https://doi.org/10.1523/JNEUROSCI.06-08-02278.1986

    • PubMed
    • Google Scholar
59. 1. Picton LD
    2. Bertuzzi M
    3. Pallucchi I
    4. Fontanel P
    5. Dahlberg E
    6. Björnfors ER
    7. Iacoviello F
    8. Shearing PR
    9. El Manira A

    (2021) A spinal organ of proprioception for integrated motor action feedback

    Neuron 109:1188–1201.

    https://doi.org/10.1016/j.neuron.2021.01.018

    • PubMed
    • Google Scholar
60. 1. Roussel Y
    2. Paradis M
    3. Gaudreau SF
    4. Lindsey BW
    5. Bui TV

    (2020) Spatiotemporal transition in the role of synaptic inhibition to the tail beat rhythm of developing larval zebrafish

    Eneuro 7:ENEURO.0508-18.2020.

    https://doi.org/10.1523/ENEURO.0508-18.2020

    • PubMed
    • Google Scholar
61. Software

    1. Roussel Y
    2. Gaudreau SF
    3. Kacer ER
    4. Topcu E
    5. Bui TV

    (2021) Zebrafish spinal locomotor circuit, version swh:1:rev:dd36ce928a2775eeed45149444962a422cb16446

    Software Heritage.

    https://archive.softwareheritage.org/swh:1:dir:ae8d2601852f0bc977d717bafa389d53125cdd94;origin=https://github.com/Bui-lab/Code;visit=swh:1:snp:690ef33a4fa168e1527c9ed3f391a6462ca40316;anchor=swh:1:rev:dd36ce928a2775eeed45149444962a422cb16446
62. 1. Ryczko D
    2. Charrier V
    3. Ijspeert A
    4. Cabelguen JM

    (2010) Segmental oscillators in axial motor circuits of the salamander: distribution and bursting mechanisms

    Journal of Neurophysiology 104:2677–2692.

    https://doi.org/10.1152/jn.00479.2010

    • PubMed
    • Google Scholar
63. Book

    1. Saint-Amant L

    (2010)

    Development of Motor Rhythms in Zebrafish Embryos

    Elsevier B.V.

    • Google Scholar
64. 1. Saint-Amant L
    2. Drapeau P

    (1998) Time course of the development of motor behaviors in the zebrafish embryo

    Journal of Neurobiology 37:622–632.

    https://doi.org/10.1002/(SICI)1097-4695(199812)37:4<622::AID-NEU10>3.0.CO;2-S

    • PubMed
    • Google Scholar
65. 1. Saint-Amant L
    2. Drapeau P

    (2000) Motoneuron activity patterns related to the earliest behavior of the zebrafish embryo

    The Journal of Neuroscience 20:3964–3972.

    https://doi.org/10.1523/JNEUROSCI.20-11-03964.2000

    • PubMed
    • Google Scholar
66. 1. Saint-Amant L
    2. Drapeau P

    (2001) Synchronization of an embryonic network of identified spinal interneurons solely by electrical coupling

    Neuron 31:1035–1046.

    https://doi.org/10.1016/S0896-6273(01)00416-0

    • PubMed
    • Google Scholar
67. 1. Satou C
    2. Kimura Y
    3. Higashijima S

    (2012) Generation of multiple classes of V0 neurons in zebrafish spinal cord: progenitor heterogeneity and temporal control of neuronal diversity

    Journal of Neuroscience 32:1771–1783.

    https://doi.org/10.1523/JNEUROSCI.5500-11.2012

    • PubMed
    • Google Scholar
68. 1. Satou C
    2. Sugioka T
    3. Uemura Y
    4. Shimazaki T
    5. Zmarz P
    6. Kimura Y
    7. Higashijima SI

    (2020) Functional diversity of glycinergic commissural inhibitory neurons in larval zebrafish

    Cell Reports 30:3036–3050.

    https://doi.org/10.1016/j.celrep.2020.02.015

    • PubMed
    • Google Scholar
69. 1. Sautois B
    2. Soffe SR
    3. Li WC
    4. Roberts A

    (2007) Role of type-specific neuron properties in a spinal cord motor network

    Journal of Computational Neuroscience 23:59–77.

    https://doi.org/10.1007/s10827-006-0019-1

    • PubMed
    • Google Scholar
70. 1. Sengupta M
    2. Daliparthi V
    3. Roussel Y
    4. Bui TV
    5. Bagnall MW

    (2021) Spinal V1 neurons inhibit motor targets locally and sensory targets distally

    Current Biology 72:053.

    https://doi.org/10.1016/j.cub.2021.06.053

    • Google Scholar
71. 1. Song J
    2. Dahlberg E
    3. El Manira A

    (2018) V2a interneuron diversity tailors spinal circuit organization to control the vigor of locomotor movements

    Nature Communications 9:1–14.

    https://doi.org/10.1038/s41467-018-05827-9

    • Google Scholar
72. 1. Song J
    2. Pallucchi I
    3. Ausborn J
    4. Ampatzis K
    5. Bertuzzi M
    6. Fontanel P
    7. Picton LD
    8. El Manira A

    (2020) Multiple Rhythm-Generating circuits act in tandem with pacemaker properties to control the start and speed of locomotion

    Neuron 105:1048–1061.

    https://doi.org/10.1016/j.neuron.2019.12.030

    • Google Scholar
73. 1. Stickney HL
    2. Barresi MJ
    3. Devoto SH

    (2000) Somite development in zebrafish

    Developmental Dynamics 219:287–303.

    https://doi.org/10.1002/1097-0177(2000)9999:9999<::AID-DVDY1065>3.0.CO;2-A

    • PubMed
    • Google Scholar
74. 1. Taylor AL
    2. Goaillard JM
    3. Marder E

    (2009) How multiple conductances determine electrophysiological properties in a multicompartment model

    Journal of Neuroscience 29:5573–5586.

    https://doi.org/10.1523/JNEUROSCI.4438-08.2009

    • PubMed
    • Google Scholar
75. 1. Tazerart S
    2. Viemari JC
    3. Darbon P
    4. Vinay L
    5. Brocard F

    (2007) Contribution of persistent sodium current to locomotor pattern generation in neonatal rats

    Journal of Neurophysiology 98:613–628.

    https://doi.org/10.1152/jn.00316.2007

    • PubMed
    • Google Scholar
76. 1. Tong H
    2. McDearmid JR

    (2012) Pacemaker and plateau potentials shape output of a developing locomotor network

    Current Biology 22:2285–2293.

    https://doi.org/10.1016/j.cub.2012.10.025

    • PubMed
    • Google Scholar
77. 1. Vogelstein JT
    2. Park Y
    3. Ohyama T
    4. Kerr RA
    5. Truman JW
    6. Priebe CE
    7. Zlatic M

    (2014) Discovery of brainwide neural-behavioral maps via multiscale unsupervised structure learning

    Science 344:386–392.

    https://doi.org/10.1126/science.1250298

    • PubMed
    • Google Scholar
78. 1. Wahlstrom-Helgren S
    2. Montgomery JE
    3. Vanpelt KT
    4. Biltz SL
    5. Peck JH
    6. Masino MA

    (2019) Glutamate receptor subtypes differentially contribute to optogenetically activated swimming in spinally transected zebrafish larvae

    Journal of Neurophysiology 122:2414–2426.

    https://doi.org/10.1152/jn.00337.2019

    • PubMed
    • Google Scholar
79. 1. Warp E
    2. Agarwal G
    3. Wyart C
    4. Friedmann D
    5. Oldfield CS
    6. Conner A
    7. Del Bene F
    8. Arrenberg AB
    9. Baier H
    10. Isacoff EY

    (2012) Emergence of patterned activity in the developing zebrafish spinal cord

    Current Biology 22:93–102.

    https://doi.org/10.1016/j.cub.2011.12.002

    • PubMed
    • Google Scholar
80. 1. Wiggin TD
    2. Anderson TM
    3. Eian J
    4. Peck JH
    5. Masino MA

    (2012) Episodic swimming in the larval zebrafish is generated by a spatially distributed spinal network with modular functional organization

    Journal of Neurophysiology 108:925–934.

    https://doi.org/10.1152/jn.00233.2012

    • PubMed
    • Google Scholar
81. 1. Wiggin TD
    2. Peck JH
    3. Masino MA

    (2014) Coordination of fictive motor activity in the larval zebrafish is generated by non-segmental mechanisms

    PLOS ONE 9:e109117.

    https://doi.org/10.1371/journal.pone.0109117

    • PubMed
    • Google Scholar
82. 1. Yang L
    2. Rastegar S
    3. Strähle U

    (2010) Regulatory interactions specifying Kolmer-Agduhr interneurons

    Development 137:2713–2722.

    https://doi.org/10.1242/dev.048470

    • PubMed
    • Google Scholar
83. 1. Zhang RW
    2. Wei HP
    3. Xia YM
    4. Du JL

    (2010) Development of light response and GABAergic excitation-to-inhibition switch in zebrafish retinal ganglion cells

    The Journal of Physiology 588:2557–2569.

    https://doi.org/10.1113/jphysiol.2010.187088

    • PubMed
    • Google Scholar
84. 1. Zhong G
    2. Masino MA
    3. Harris-Warrick RM

    (2007) Persistent sodium currents participate in fictive locomotion generation in neonatal mouse spinal cord

    Journal of Neuroscience 27:4507–4518.

    https://doi.org/10.1523/JNEUROSCI.0124-07.2007

    • PubMed
    • Google Scholar

Acceptance summary:

As animals mature, their locomotor patterns become varied, more flexible and complex. In this manuscript, Roussel et al., build models of the spinal network in embryonic and larval zebrafish based on experimental data to understand how these animals generate distinct behaviors during development.

Decision letter after peer review:

Thank you for submitting your article "Modelling spinal locomotor circuits for movements in developing zebrafish" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Ronald Calabrese as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Rishikesh Narayanan (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Reviewer #1 (Recommendations for the authors):

For the B panels in figures 2, 3 and 4, are these responses of isolated neurons or neurons in the fully connected network?

Wherever neuron traces are shown, showing rostrocaudal and left-right neuron traces in the same panel makes them too cluttered. In all cases, where the firing starts and stops is difficult to discern. It may suffice to show one rostral and one caudal neuron trace.

Single coiling:

Include sensitivity analysis in graphical form in Figure 2.

SBs are depolarizing and reverse at -40. How critical are they then for preventing simultaneous contralateral coiling? Why an analysis of neuron deletions has not been done?

Was sensitivity different for glycinergic synapses and electrical synapses?

Was bursting in PMs critical for generating single coiling? This can be included in Figure 2. Figure 2 can be improved to test the single coiling model in multiple ways. In its current form, it doesn't provide much data regarding the performance of the model.

Double coiling:

Line 219: what proportion of PMs were electrically coupled to rostral CENs?

The neuron traces are confusing with too many superimposed.

Line 240: Replace 'with' with 'to'.

What accounts for the long delay in firing between PMs and other neurons though they are electrically coupled? These seem much longer than segmental conduction delays. Such long delays seem critical for CENs to initiate the second coil on the contralateral side once the first coil is over. Yet, this seems inexplicable following known monosynaptic/conduction delays. Likewise, in this circuit configuration, what prevents the circuit from going into third or fourth coils? It is not clear which process decays slowly to disallow such multiple coils. Similarly, it is not clear how under conditions of glycinergic blockade in the model, multiple coilings are generated if inhibition is normally over before CENs are even activated.

Figure 3C-E: Why aren't there multiple traces for IEDs?

Figure 3H: Convert to color for better effect.

Beat and glide swimming:

Line 398: Do you mean 'few'?

Line 431: "Silencing…' – this sentence does not make sense.

How is swimming initiated?

Figure 4A: This schematic is too confusing. Do MNs receive only IED inputs? Do they not receive direct inhibition?

Figure 4C, F: Unable to discern any detail of MN activity due to overlap of traces. Show phase relationships between classes of neurons and comment.

Figure 4D: Same as above, cannot see the color coded information clearly as most angles appear red or blue and therefore are at the extremes of their scalebar. Why aren't there any values for body position 1.0 which should be the tail tip?

Figure 5B: Authors state in line 432 that removing CENs does not affect rhythmic firing in IEDs or MNs. When CENs are removed, is there left-right alternation? From the traces shown, there doesn't appear to be and the video shows what appear to be unilateral tail flicks. I couldn't tell for sure as they are at a fast playback speed. If that is the case, authors need to be more convincing of this result in the figure and the video and perhaps include L-R phase relationships to establish that alternation occurs. What do the peaks in Row3 correspond to? Is each cycle a bout? Seems too slow to be a tail beat? Do the dampened oscillations in the middle epoch qualify as bonafide swimming? If yes, what is the justification? Also please explain, 'sum of left and right motoneuron activity was summed'. Why use this parameter as motor output when you have the muscle transform calculated? Why is MN activity shown on a different scale compared to the other rows below?

Figure 7: The importance of bursting IED is not really clarified in this figure. Is there a quantitative difference in durations? Both long and short are produced. Why? Perhaps a phase plane analysis of contribution of IED parameters to durations will help sort out the underlying dynamics? Young prefeeding larvae do exhibit long duration swims, so the conclusion that bursting IEDs are essential for beat and glide swimming appears too strong for the results presented here.

The videos shown for CIN null and strychnine experiments show left-right tail beats while figure panels show different results.

Discussion:

Lines 604-635: Authors discuss differential recruitment of spinal neurons as a function of 'speed'. However, the studies cited looked at the recruitment of neurons as a function of frequency. It has recently become clear that speed can be changed even without changing swim frequency and that MNs are recruited as a function of forward speed and not frequency (Jha and Thirumalai 2020).

Methods:

A more detailed methods section is warranted.

It is not clear what the authors mean by tonic motor command- is it a DC current injection or a regular synaptic input?

Line 705: What is s here?

Line 739: Why is capacitance decreasing with age?

Line 745: How is deflection angle calculated per somite when three somites are collapsed into one muscle?

Line 756: 'muscle output' is vague – do you mean V or theta?

Figure 1D: The transformation from local body angle to an overall body midline position can be better schematized to help the reader.

Figure 1E: More intuitive to have time on x-axis as has been done in Figure 7F. Change elsewhere also.

Reviewer #2 (Recommendations for the authors):

1. The sensitivity analyses involving variability in individual parameters is useful. But, it is not clear from the description (starting line 721) on what specific parameters are governed by \σ_p, \σ_w and \σ_d. Importantly, these sensitivity analyses do not seem to cover cell-to-cell variability in Izhikevich model parameters (spanning a, b, c, d, k, C, V_r, V_t; it is not clear if \σ_p alters any of these parameters?). Specifically, all neurons of the same subtype seem to have the same model parameters. While the firing dynamics of Izhikevich model neuron are critically dependent on the model parameters, there are a range of model parameters (perturbations around the employed values) over which similar firing patterns could be achieved. How sensitive are the conclusions presented here to such variability across neurons of the same subtype in the network (with the variability reflecting electrophysiological recordings from larval neurons)?

What would be the impact of heterogeneities in gap junction connectivity? From the description provided, it seems like the model assumes identical weights for gap junction connectivity. How sensitive would the conclusions presented here to adding variability to gap junctional weights?

Additional sensitivity analyses could focus on the numbers of neurons per chain, number of connections within and across chains and glycinergic synaptic connectivity (there seem to be several unsubstantiated assumptions on these connections; e.g. paragraph spanning lines 126-137).

Apart from these, any parameter that has been fixed to be homogeneous should also be considered in addressing the question on whether the results observed are because of the homogeneous nature of that parameter in the simulations. The authors could consider incorporating some of these simulations involving networks with parametric variability, and mention the others in the discussion.

2. All interpretations of the authors from experimental data are based on the summary statistics, and do not account for heterogeneities across different larvae. For instance, the authors cite experimental data that adding CNQX and APV to block glutamatergic transmission precluded double coils while sparing single coils, and that blocking glycinergic synapses led to triple or even quadruple coils. They use these conclusions to drive their modeling outcomes. However, are the impacts of CNQX or APV or blockade of glycinergic synapses the same across all larvae? Do all larvae behave identically when treated with these pharmacological agents? If not, how do the authors account for these heterogeneities across larvae.

Employing summary statistics to define models or assess outcomes is perilous, and the authors should account for heterogeneities in experimental outcomes to define their models (Marder and Taylor, Nature Neuroscience, 2011). Heterogeneities and variability should be accounted for each measurement at each scale in interpreting modeling and experimental outcomes -- the CNQX/APV/Glycinergic synaptic blockers is just one instance. The ideal way to do this is to generate a population of heterogeneous models and derive conclusions from there (Marder and Taylor, Nature Neuroscience, 2011). I recommend that the authors mention this a future direction, refine their interpretations in the Discussion section accounting for heterogeneities across larvae and present the caveats of using summary statistics to drive experimental/model interpretations.

3. Are there several routes (in terms of neuronal intrinsic properties, neuronal firing patterns, network connectivity, electrical/chemical synapse weights, neuromodulatory tones, etc.) to achieving single coiling or double coiling or other locomotor movements? Or, is there a unique route to achieve these across larvae? Could the authors comment on potential degeneracy across scales in achieving the same behavioral outcomes (see Vogelstein et al., Science, 2014; Edelman and Gally, PNAS, 2001)? How would such degeneracy alter the conclusions presented here? How would such degeneracy relate to heterogeneities across larvae (Pt. 2 above) and variability in parameters (Pt. 1 above)? The authors should provide detailed discussion on these questions. They have considered single hand-tuned models across different developmental stages to make their points and demonstrate behavioral outcomes. But, it is important to ask if this is the only route to achieve these behavioral outcomes, and if zebrafish and their larvae might be using different structural combinations to elicit same functional outcomes?

While the authors have talk about potential ion-channel degeneracy in the discussion, they should also comment on circuit level degeneracy involving different neuronal subtypes and disparate synaptic components (in terms of neuronal intrinsic properties, neuronal firing patterns, network connectivity, electrical/chemical synapse weights, neuromodulatory tones, etc. e.g. Prinz et al., Nat Neuroscience, 2004). The authors could mention in the Discussion section about the possibility of using an unbiased stochastic search involving all parameters to assess such potential degeneracy involving intrinsic properties and synaptic connectivity across all stages of the developing fish. The authors could also mention in the discussion the disadvantages of using a single hand-tuned model and in deriving all interpretations from that single model (Marder and Taylor, Nat Neuroscience., 2011).

4. Please add a separate Discussion section on testable predictions that emerge from the model presented here, which provide clear pointers that would test the overall hypothesis.

Reviewer #3 (Recommendations for the authors):

Strengths and weaknesses

The key strength of this manuscript is the detailing of a set of related models detailing the motor output of the larval zebrafish across key stages of development. The models should form a basis for future research. It also a first of its kind – I don't know of similar models focusing on development of locomotor function. The main weakness is the reliance on assumptions of model connectivity. But I suggest that if the model is treated as a basis for the community to refine and validate it will be incredibly useful.

Reviewer #1 (Recommendations for the authors):

For the B panels in figures 2,3 and 4, are these responses of isolated neurons or neurons in the fully connected network?

The responses are from isolated neurons, and this has been clarified in the figure legends of Figures 2-4 and 8.

Wherever neuron traces are shown, showing rostrocaudal and left-right neuron traces in the same panel makes them too cluttered. In all cases, where the firing starts and stops is difficult to discern. It may suffice to show one rostral and one caudal neuron trace.

We have reduced the number of traces in figures depicting neural activity during simulations to a rostral, middle, and caudal neuron.

Single coiling:

Include sensitivity analysis in graphical form in Figure 2.

We have added the sensitivity analysis for the single coiling model to Figure 2.

SBs are depolarizing and reverse at -40. How critical are they then for preventing simultaneous contralateral coiling? Why an analysis of neuron deletions has not been done?

Due to the lack of contralateral excitation, there are no contralateral coils in the single coiling model. We further emphasize this point in a new simulation where we remove contralateral inhibition by V0d, which removes SBs, and find no emergence of multiple coilings (Figure 2 —figure supplement 1). We have added an analysis of neuron deletions for ICs and MNs in Figure 2 —figure supplement 1.

Was sensitivity different for glycinergic synapses and electrical synapses?

We did not perform this analysis for the single coiling model as glycinergic synapses do not contribute to the generation of coils or the attenuation of double or multiple coils in the single coiling model (see our response to the previous point). For the double coiling model, we did perform a sensitivity analysis of chemical versus electrical synapses (Figure 3 —figure supplement 2, Lines 356-360).

Was bursting in PMs critical for generating single coiling? This can be included in Figure 2. Figure 2 can be improved to test the single coiling model in multiple ways. In its current form, it doesn't provide much data regarding the performance of the model.

The bursting in the ICs (formerly PMs) is not essential to single coiling. In response to this question, we created a model with ICs firing a single spike per burst. With an increase in the synaptic weight from IC to MNs, this generated single coiling as well. However, we modelled ICs as bursting to resemble the recordings of ICs from Tong and McDearmid (2012), which show sustained bursts of short spikelets. To address this question, we have clarified the reason why our IC models display bursting. We can add the result of the one-spike bursting ICs as a figure supplement to Figure 2 if desired.

Double coiling:

Line 219: what proportion of PMs were electrically coupled to rostral CENs?

All ipsilateral ICs (we now refer to PMs as their biological equivalents, the IC neurons) were coupled with the most rostral V0vs (formerly CENs) in the first four rostral somites. We have clarified this point in the text (Line 270).

The neuron traces are confusing with too many superimposed.

As mentioned above, we have reduced the number of traces to a rostral, middle, and caudal neuron.

Line 240: Replace 'with' with 'to'.

This change has been made (Line 268).

What accounts for the long delay in firing between PMs and other neurons though they are electrically coupled? These seem much longer than segmental conduction delays. Such long delays seem critical for CENs to initiate the second coil on the contralateral side once the first coil is over. Yet, this seems inexplicable following known monosynaptic/conduction delays. Likewise, in this circuit configuration, what prevents the circuit from going into third or fourth coils? It is not clear which process decays slowly to disallow such multiple coils. Similarly, it is not clear how under conditions of glycinergic blockade in the model, multiple coilings are generated if inhibition is normally over before CENs are even activated.

For double coiling, the critical delay is the delay between IC firing and the subsequent activation of the contralaterally-projecting glutamatergic V0v (formerly CENs). This delay is a sum of the delay in ICs activating V2a (formerly IEDs) and the latter activating V0v through glutamatergic excitation. The V2a to V0v excitation has to be titrated carefully to produce a delay. Otherwise, there is no second contralateral coil if the first coil is not allowed to be completed. To highlight this delay, we have added a phase delay plot (Figure 3D) and also used simulations where the V2a to V0v is increased (Figure 3E) to support these points further (Lines 305-312).

Figure 3C-E: Why aren't there multiple traces for IEDs?

In our model, the more caudal V2as (formerly IEDs) are not activated. This lack of activity in the caudal V2as does not limit the ability to generate double coils. The rostral V0vs are crucial for this behaviour and are recruited by rostral V2as.

Figure 3H: Convert to color for better effect.

This change has been made (now Figure 3I), and we use this color scheme for new figures depicting additional sensitivity testing with the double coiling model (Figure 3J-M and Figure 3 —figure supplement 2).

Beat and glide swimming:

Line 398: Do you mean 'few'?

No, by ‘new’, we meant that this connection from V2as and V1s was new. However, since V1s were new, any connections to and from them are new. To avoid any confusion, we removed “new”. “Reciprocally, V2as formed glutamatergic synapses to caudally located V1.” (Line 429).

Line 431: "Silencing…' – this sentence does not make sense.

This sentence has been modified to better convey our message:

“Silencing V0vs diminished but did not eliminate the rhythmic firing of V2as or MNs. During Epoch 2, the tonic motor command continued to activate V2as. Pairs of left-right tail beats may result from the commissural inhibition by dI6s that was still present. However, removing the contralateral excitation by V0v prevented the repetitive activation of the silent side after each tail beat, which severely reduced episode duration and the number of tail beats generated in each episode (Figure 5G-L, Figure 5 —figure supplement 1G-L, Figure 5 – video 1).” (Lines 505-510).

How is swimming initiated?

The new phase delay plot that we have added suggests that swimming is initiated by V2as that serve to drive all other spinal neurons in each tail beat (Figure 4H and Lines 486-489).

Figure 4A: This schematic is too confusing. Do MNs receive only IED inputs? Do they not receive direct inhibition?

The schematic for the beat and glide model has been modified. The inhibition of MNs by V1s (formerly IIA) has been added for clarification.

Figure 4C, F: Unable to discern any detail of MN activity due to overlap of traces. Show phase relationships between classes of neurons and comment.

The phase relationships have been added in Figure 4H, and these results are commented on (Lines 486-496).

Figure 4D: Same as above, cannot see the color coded information clearly as most angles appear red or blue and therefore are at the extremes of their scalebar. Why aren't there any values for body position 1.0 which should be the tail tip?

We agree that the colour-coding of the angles does not allow for a precise depiction of the angles. Because of the time range depicted, we use this type of heatmap (also used in Figure 8 and Figure 8 —figure supplement 1) mainly to show which somites are active, depict the swimming pattern and demonstrate that episodes consist of left-right tail beats. The more precise changes in angles across the somites are better viewed in the videos, which have been slowed down to observe the kinematics better.

The omission of the values for body position at the tail tip for the local body angle heat maps has been corrected, and they are now present.

Figure 5B: Authors state in line 432 that removing CENs does not affect rhythmic firing in IEDs or MNs. When CENs are removed, is there left-right alternation? From the traces shown, there doesn't appear to be and the video shows what appear to be unilateral tail flicks. I couldn't tell for sure as they are at a fast playback speed. If that is the case, authors need to be more convincing of this result in the figure and the video and perhaps include L-R phase relationships to establish that alternation occurs. What do the peaks in Row3 correspond to? Is each cycle a bout? Seems too slow to be a tail beat? Do the dampened oscillations in the middle epoch qualify as bonafide swimming? If yes, what is the justification? Also please explain, 'sum of left and right motoneuron activity was summed'. Why use this parameter as motor output when you have the muscle transform calculated?

We have revised the text to clarify that there is no left-right alternation without commissural excitation. Furthermore, we have added a supplementary figure with membrane potential traces of all neurons to clarify that while there is still rhythmic activity in V2a and MNs, this activity is greatly diminished, and swimming episodes are considerably shortened. We have added supplementary figures with membrane potential traces of all neurons for all the knockout simulations (Figure 5 —figure supplement 1, and Figure 6 —figure supplement 2).

We have corrected the statement in the original version of the manuscript regarding using the sum of left and right motoneuron activity. It is the sum of left and right muscle output.

Why is MN activity shown on a different scale compared to the other rows below?

We show MN activity on a different scale to show a magnified view of MN activity during one or two swim episodes. In the rows below, the motor output and parameters of swimming activity are displayed for the entire simulation. This change in time scale is to show fluctuations of this output during the progression of the simulations.

Figure 7: The importance of bursting IED is not really clarified in this figure. Is there a quantitative difference in durations? Both long and short are produced. Why? Perhaps a phase plane analysis of contribution of IED parameters to durations will help sort out the underlying dynamics? Young prefeeding larvae do exhibit long duration swims, so the conclusion that bursting IEDs are essential for beat and glide swimming appears too strong for the results presented here.

As described in the introduction to our Author’s Response, the modified beat-and-glide swimming model does not suggest that bursting IEDs are essential anymore, so this section has been removed.

The videos shown for CIN null and strychnine experiments show left-right tail beats while figure panels show different results.

While the figures show reduced left-right alternation in the dI6 (formerly CIN null) and strychnine simulations, the videos do indeed still show left-right tail beats (Figure 6 – video 2 and Figure 7 – video 1). The conserved left-right alternation in these simulations is explained by the presence of some left-right alternation in rostral segments. With the left-right alternation in rostral segments, the rest of the body will inevitably sway from left to right. We have added the left-right cross-correlation for a rostral somite to show that silencing dI6 and strychnine have a lesser effect on rostral segments, which drives the continued presence of left-right tail beats. We have also added two figure supplements, Figure 6 —figure supplement 2 and Figure 7 —figure supplement 1, to illustrate that the kinematics of the tail beats are modified in the dI6 null and the strychnine simulations by the reduced, but not eliminated, left-right alternation.

To the best of our knowledge, the experimental work with dI6 silencing (Satou et al., 2020) showed similar levels of reduced left-right alternation as we do. Their results were obtained from ex-vivo recordings. They do not report the effect on actual swimming, so we do not know whether their transgenic fish still show some left-right tail beats. As for strychnine, the effects of strychnine on actual fish swimming are complicated by its effects on neuromuscular junction transmission (E.g. García-Colunga and Miledi, PNAS, 1999,https://doi.org/10.1073/pnas.96.7.4113).

Discussion:

Lines 604-635: Authors discuss differential recruitment of spinal neurons as a function of 'speed'. However, the studies cited looked at the recruitment of neurons as a function of frequency. It has recently become clear that speed can be changed even without changing swim frequency and that MNs are recruited as a function of forward speed and not frequency (Jha and Thirumalai 2020).

Thank you for clarifying this important distinction. We now refer specifically to swimming frequency rather than swimming speed. We have modified the text to reflect this distinction.

Methods:

A more detailed methods section is warranted.

We have made the methods more detailed. In particular, we have added details about how the muscle activity is converted to body angles (Lines 934 – 950) and how swimming activity was analyzed (e.g. phase delays, episode durations, inter-episode intervals, etc) at Lines 952-974.

It is not clear what the authors mean by tonic motor command- is it a DC current injection or a regular synaptic input?

We modelled the tonic motor command as a DC current injection. This is clarified at Lines 456-458 While this is less biologically accurate, we made this choice to reduce computational load.

Line 705: What is s here?

The s stands for second, and we have clarified this:

“…the transmission speed in arbitrary distance units per second (a.d.u./s)” (Line 883).

Line 739: Why is capacitance decreasing with age?

These differences in capacitance values are a weakness of our simple musculoskeletal model. We found that the lower capacitance values enabled shorter time constants for the high-frequency tail beats instead of the slower coiling movements' speeds. The lower capacitance may reflect differences in muscle types active during coiling versus swimming or kinetics of neuromuscular transmission at stages of development.

Line 745: How is deflection angle calculated per somite when three somites are collapsed into one muscle?

We have clarified the notion of model muscle cell and biological somite so that the deflection angle is calculated for one muscle cell in which are collapsed three somites (coiling models) or two somites (swimming model) (Lines 125-127, 398-400, 922-924).

Line 756: 'muscle output' is vague – do you mean V or theta?

This issue was clarified to mean V:

“Left-right alternation at a particular somite was analyzed using cross-correlation of VRmuscle,i and VLmuscle,ith left and right muscle output at that somite.” (Lines 946-947).

Figure 1D: The transformation from local body angle to an overall body midline position can be better schematized to help the reader.

We have modified this figure to make the transformation more intuitive.

Figure 1E: More intuitive to have time on x-axis as has been done in Figure 7F. Change elsewhere also.

Time has been placed on the x-axis as suggested for this and other similar figures (Figure 8 and Figure 8 —figure supplement 1).

Reviewer #2 (Recommendations for the authors):

1. The sensitivity analyses involving variability in individual parameters is useful. But, it is not clear from the description (starting line 721) on what specific parameters are governed by \σ_p, \σ_w and \σ_d. Importantly, these sensitivity analyses do not seem to cover cell-to-cell variability in Izhikevich model parameters (spanning a, b, c, d, k, C, V_r, V_t; it is not clear if \σ_p alters any of these parameters?). Specifically, all neurons of the same subtype seem to have the same model parameters. While the firing dynamics of Izhikevich model neuron are critically dependent on the model parameters, there are a range of model parameters (perturbations around the employed values) over which similar firing patterns could be achieved. How sensitive are the conclusions presented here to such variability across neurons of the same subtype in the network (with the variability reflecting electrophysiological recordings from larval neurons)?

What would be the impact of heterogeneities in gap junction connectivity? From the description provided, it seems like the model assumes identical weights for gap junction connectivity. How sensitive would the conclusions presented here to adding variability to gap junctional weights?

Additional sensitivity analyses could focus on the numbers of neurons per chain, number of connections within and across chains and glycinergic synaptic connectivity (there seem to be several unsubstantiated assumptions on these connections; e.g. paragraph spanning lines 126-137).

Apart from these, any parameter that has been fixed to be homogeneous should also be considered in addressing the question on whether the results observed are because of the homogeneous nature of that parameter in the simulations. The authors could consider incorporating some of these simulations involving networks with parametric variability, and mention the others in the discussion.

As suggested by the reviewer, we have added several sensitivity analyses. These analyses include sensitivity to variability in just the chemical synapses as well as just the gap-junction weights (double coiling model), the rostrocaudal extent of axon projections (all models), and glycinergic reversal potential (double coiling and beat-and-glide swimming models). We have also tested our ability to generate single and double coiling with a thirty-somite model and beat-and-glide with a thirty-somite model and a shorter 10-somite model.

We have also clarified how the sensitivity analysis was performed to test variability of the parameters setting membrane potential dynamics (σp) "…the parameters that set the dynamics of the membrane potential of each neuron (a, b, c, d, and Vmax, k, C, Vr and Vt , see Materials and methods)" (Lines 204-207).

With regards to homogeneity of the parameters in our model, it is true that in the base model, the parameters are homogeneous for each cell within a population or each specific synaptic connection. In our sensitivity analysis, however, every parameter (whether setting membrane potential dynamics, synaptic weight) is scaled individually. We try to describe this “Variability in the amplitude of the tonic motor command, the rostrocaudal extent of every axonal projection, every parameter that set the dynamics of the membrane potential of each neuron (a, b, c, d, and Vmax, k, C, Vr and Vt), and all of the weights of gap junction and chemical synapses were modelled by scaling each value by a random number picked for each simulation” (Lines 204-208).

We now discuss the homogeneity of some of the parameters in our model and the implications of this homogeneity and other sensitivity analyses to be considered in the future (Lines 676-689).

2. All interpretations of the authors from experimental data are based on the summary statistics, and do not account for heterogeneities across different larvae. For instance, the authors cite experimental data that adding CNQX and APV to block glutamatergic transmission precluded double coils while sparing single coils, and that blocking glycinergic synapses led to triple or even quadruple coils. They use these conclusions to drive their modeling outcomes. However, are the impacts of CNQX or APV or blockade of glycinergic synapses the same across all larvae? Do all larvae behave identically when treated with these pharmacological agents? If not, how do the authors account for these heterogeneities across larvae.

Employing summary statistics to define models or assess outcomes is perilous, and the authors should account for heterogeneities in experimental outcomes to define their models (Marder and Taylor, Nature Neuroscience, 2011). Heterogeneities and variability should be accounted for each measurement at each scale in interpreting modeling and experimental outcomes -- the CNQX/APV/Glycinergic synaptic blockers is just one instance. The ideal way to do this is to generate a population of heterogeneous models and derive conclusions from there (Marder and Taylor, Nature Neuroscience, 2011). I recommend that the authors mention this a future direction, refine their interpretations in the Discussion section accounting for heterogeneities across larvae and present the caveats of using summary statistics to drive experimental/model interpretations.

Thank you for raising the importance of accounting for inter-animal variability in our modelling work. This issue is now discussed (Lines 683-689).

3. Are there several routes (in terms of neuronal intrinsic properties, neuronal firing patterns, network connectivity, electrical/chemical synapse weights, neuromodulatory tones, etc.) to achieving single coiling or double coiling or other locomotor movements? Or, is there a unique route to achieve these across larvae? Could the authors comment on potential degeneracy across scales in achieving the same behavioral outcomes (see Vogelstein et al., Science, 2014; Edelman and Gally, PNAS, 2001)? How would such degeneracy alter the conclusions presented here? How would such degeneracy relate to heterogeneities across larvae (Pt. 2 above) and variability in parameters (Pt. 1 above)? The authors should provide detailed discussion on these questions. They have considered single hand-tuned models across different developmental stages to make their points and demonstrate behavioral outcomes. But, it is important to ask if this is the only route to achieve these behavioral outcomes, and if zebrafish and their larvae might be using different structural combinations to elicit same functional outcomes?

While the authors have talk about potential ion-channel degeneracy in the discussion, they should also comment on circuit level degeneracy involving different neuronal subtypes and disparate synaptic components (in terms of neuronal intrinsic properties, neuronal firing patterns, network connectivity, electrical/chemical synapse weights, neuromodulatory tones, etc. e.g. Prinz et al., Nat Neuroscience, 2004). The authors could mention in the Discussion section about the possibility of using an unbiased stochastic search involving all parameters to assess such potential degeneracy involving intrinsic properties and synaptic connectivity across all stages of the developing fish. The authors could also mention in the discussion the disadvantages of using a single hand-tuned model and in deriving all interpretations from that single model (Marder and Taylor, Nat Neuroscience., 2011).

We agree that there is likely a substantial amount of degeneracy in networks capable of achieving the motor output desired. We show some of this degeneracy by showing that beat and glide can be generated by models with bursting V2a (formerly IED) or V0v (formerly CENs), as well as just tonic neurons. We now discuss this at greater length, in part in light of inter-animal heterogeneity and variability (Lines 742-757).

4. Please add a separate Discussion section on testable predictions that emerge from the model presented here, which provide clear pointers that would test the overall hypothesis.

As suggested by the reviewer, a section on testable predictions has been added to the Discussion (Lines 706-757).

Author details

1. Yann Roussel

   1. Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Biology, University of Ottawa, Ottawa, Canada
   2. Blue Brain Project, École Polytechnique Fédérale de Lausanne, Genève, Switzerland

   Contribution

   Conceptualization, Resources, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

   Competing interests

   No competing interests declared

2. Stephanie F Gaudreau

   Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Biology, University of Ottawa, Ottawa, Canada

   Contribution

   Data curation, Validation, Investigation, Writing - review and editing

   Competing interests

   No competing interests declared

3. Emily R Kacer

   Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Biology, University of Ottawa, Ottawa, Canada

   Contribution

   Software, Investigation, Writing - review and editing

   Competing interests

   No competing interests declared

4. Mohini Sengupta

   Washington University School of Medicine, Department of Neuroscience, St Louis, United States

   Contribution

   Validation, Methodology

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-5234-8258

5. Tuan V Bui

   Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Biology, University of Ottawa, Ottawa, Canada

   Contribution

   Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Project administration, Writing - review and editing

   For correspondence

   tuan.bui@uottawa.ca

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0003-0024-1544

• 1,782

  Page views

• 255

  Downloads

• 6

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.