1. Neuroscience
Download icon

Computational modeling of spinal circuits controlling limb coordination and gaits in quadrupeds

  1. Simon M Danner  Is a corresponding author
  2. Natalia A Shevtsova
  3. Alain Frigon
  4. Ilya A Rybak
  1. Drexel University College of Medicine, United States
  2. Université de Sherbrooke, Canada
Research Article
  • Cited 5
  • Views 931
  • Annotations
Cite as: eLife 2017;6:e31050 doi: 10.7554/eLife.31050

Abstract

Interactions between cervical and lumbar spinal circuits are mediated by long propriospinal neurons (LPNs). Ablation of descending LPNs in mice disturbs left-right coordination at high speeds without affecting fore-hind alternation. We developed a computational model of spinal circuits consisting of four rhythm generators coupled by commissural interneurons (CINs), providing left-right interactions, and LPNs, mediating homolateral and diagonal interactions. The proposed CIN and diagonal LPN connections contribute to speed-dependent gait transition from walk, to trot, and then to gallop and bound; the homolateral LPN connections ensure fore-hind alternation in all gaits. The model reproduces speed-dependent gait expression in intact and genetically transformed mice and the disruption of hindlimb coordination following ablation of descending LPNs. Inputs to CINs and LPNs can affect interlimb coordination and change gait independent of speed. We suggest that these interneurons represent the main targets for supraspinal and sensory afferent signals adjusting gait.

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

Introduction

To effectively move in the environment, limbed animals use different gaits depending on the desired locomotor speed (Grillner, 1981; Hildebrand, 1989). Mice change gait from walk to trot followed by gallop and bound as locomotor speed increases (Clarke and Still, 1999; Herbin et al., 2004; Herbin et al., 2007; Bellardita and Kiehn, 2015; Lemieux et al., 2016). Each gait is characterized by a distinct pattern of limb movements (Hildebrand, 1989), which in turn is defined by neural interactions between spinal circuits controlling each limb.

It is generally accepted that each limb is controlled by a separate central pattern generator (CPGs), because cats are able to walk on split-belt treadmills with limbs stepping at different speeds (Forssberg et al., 1980; Halbertsma, 1983; Frigon et al., 2013, 2015; Thibaudier et al., 2013; Frigon, 2017). Hence, quadrupedal locomotion is controlled by four CPGs (one per limb) located in separate sections of the spinal cord, specifically on the left and right side in the lumbar and cervical enlargements (Kato, 1990; Ballion et al., 2001; Juvin et al., 2005; Juvin et al., 2012). Therefore, the central control of interlimb coordination and gait expression is defined by neurons and neuronal pathways that mediate interactions between circuits controlling each limb. These interlimb pathways are mediated by (a) short projecting commissural interneurons (CINs) that control left-right interactions at the lumbar and, presumably, cervical levels (Talpalar et al., 2013; Bellardita and Kiehn, 2015; Molkov et al., 2015; Shevtsova et al., 2015; Danner et al., 2016), and (b) descending (cervical-to-lumbar) and ascending (lumbar-to-cervical) long propriospinal neurons (LPNs) that control homolateral and diagonal cervical-lumbar interactions (Matsushita et al., 1979; Skinner et al., 1979, 1980; Menétrey et al., 1985; Nathan et al., 1996; Dutton et al., 2006; Brockett et al., 2013; Ruder et al., 2016; Flynn et al., 2017).

The role and speed-dependent involvement of genetically identified CINs, specifically V0V and V0D, in alternation and expression of different gaits have been experimentally investigated (Talpalar et al., 2013; Bellardita and Kiehn, 2015). Specifically, it was shown that mutants lacking V0V CINs selectively lose trot, and mutants lacking both V0V and V0D CINs can only bound at any speed (Bellardita and Kiehn, 2015; Kiehn, 2016). Using modeling (Danner et al., 2016), we presented a possible organization of the central pathways mediating interlimb coordination that exhibited speed-dependent gait expression and the correct loss of gaits after selective removal of V0V and both V0V and V0D CINs. These V0 interneurons may contain both short projecting CINs as well as diagonally projecting LPNs (Ruder et al., 2016).

LPNs that directly couple the cervical and lumbar segments have been anatomically identified (Edinger, 1896; Nathan and Smith, 1959; Matsushita et al., 1979; Skinner et al., 1979; Skinner et al., 1980; Menétrey et al., 1985; Nathan et al., 1996; Dutton et al., 2006; Brockett et al., 2013; Flynn et al., 2017). Several studies in vitro found that LPNs can be involved in coordination of left-right and cervical-lumbar activities recorded from four ventral roots (Ballion et al., 2001; Juvin et al., 2005; Juvin et al., 2012; Cowley et al., 2010). However, their role in speed-dependent gait expression in vivo, remains poorly understood. Recently, Ruder et al. (2016) described axonal projections of different LPNs and their synaptic and genetic properties. These data differ from the connections proposed in our previous model (Danner et al., 2016). In addition, Ruder et al. (2016) have shown that deletion of descending (cervical-to-lumbar) LPNs causes a speed-dependent distortion of left-right coordination. These new findings provided us with the information needed to revise our model and investigate the role of different LPN pathways in interlimb coordination and gait control during locomotion.

In locomotion induced by stimulation of the mesencephalic locomotor region (MLR), increasing the intensity of stimulation causes an increase of locomotor speed accompanied by sequential gait transitions (Orlovskiĭ et al., 1966; Shik et al., 1966; Shik and Orlovsky, 1976; Skinner and Garcia-Rill, 1984; Grillner, 1985; Atsuta et al., 1990; Nicolopoulos-Stournaras and Iles, 2009). From these observations, we hypothesized that the same brainstem drive that controls speed by exciting the rhythm generators also affects interlimb coordination by acting on CINs or LPNs (Danner et al., 2016). This hypothesis is in accordance with multiple data confirming that different CINs and LPNs receive supraspinal inputs (Lloyd, 1942; Lloyd and McIntyre, 1948; Bannatyne et al., 2003; Jankowska et al., 2003; Cowley et al., 2010; Zaporozhets et al., 2011). However, the involvement of these neurons in mediating brainstem control of gait has not been demonstrated. Therefore, organization and role of brainstem inputs to these neurons can so far only be based on predictions derived from computational models.

Here, we present a computational model of the spinal locomotor circuits and investigate the potential role of CINs and LPNs in the control of and transition between different locomotor gaits. The model reproduces the speed-dependent gait expression in intact and genetically transformed mice, lacking V2a, V0V, and all V0 neurons (Crone et al., 2009; Bellardita and Kiehn, 2015), as well as the disruption of hindlimb coordination following ablation of descending LPNs (Ruder et al., 2016). The model proposes the following roles for LPNs: diagonal V0D LPNs supports walk, diagonal V0V LPNs together with local V0V CINs stabilize trot, and homolateral LPNs ensure fore-hind alternation in all gaits. Finally, we show that additional external inputs to CINs and LPNs may affect interlimb coordination and gait expression independent of speed. We suggest that these spinal interneurons are the main targets of supraspinal and somatosensory afferents to adjust interlimb coordination and gait to different environmental and behavioral conditions.

Results

The model

The model consists of four rhythm generators (RGs), each controlling one limb (Figure 1). These RGs interact via several bidirectional pathways mediated by CINs (left-right) and LPNs (homolateral and diagonal). Neural populations were modeled with a simplified, 'activity-dependent’ population model (see Materials and methods). Each RG consisted of a flexor (F) and extensor center (E) that mutually inhibited each other through intermediate inhibitory interneurons (InF and InE, respectively). Similar to our previous models (Rybak et al., 2006a, 2006b, 2013, 2015; Molkov et al., 2015; Shevtsova et al., 2015; Danner et al., 2016), both centers incorporated a slowly inactivating, persistent sodium current (INaP), allowing each center to intrinsically generate rhythmic bursting under certain conditions defined by external tonic drives or level of excitation, which was consistent with experimental data (Hägglund et al., 2013). However, we followed the concept of asymmetric flexion-dominated CPG organization (Pearson and Duysens, 1976; Duysens, 1977; Zhong et al., 2012; Duysens et al., 2013; Machado et al., 2015). Based on this concept and our previous models (Rybak et al., 2015; Molkov et al., 2015; Shevtsova et al., 2015; Danner et al., 2016; Ausborn et al., 2017), we assumed that under normal conditions extensor centers receive relatively high drive that keeps them in the regime of tonic activity. Hence, they exhibited rhythmic bursting only due to rhythmic inhibition from the corresponding intrinsically oscillating flexor centers. In this case, the frequency of oscillation in the model was defined by the brainstem drive to flexor centers and was almost independent of the phasic interactions between RGs mediated by CINs and LPNs (Rybak et al., 2015).

Model concept.

Each limb is controlled by a separate rhythm generator (RG). The local commissural neurons (CINs) as well as homolaterally and diagonally projecting (descending and ascending) long propriospinal neurons (LPNs) couple the four RGs. Brainstem drive acts on the RGs to control locomotor speed and on CINs and diagonal LPNs to control gaits.

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

Left-right CIN connections at cervical and lumbar levels (Figure 2A and Figure 3) are organized in accordance to our previous models (Rybak et al., 2015; Shevtsova et al., 2015; Danner et al., 2016): the inhibitory V0D CINs provide direct mutual inhibition between the left and right flexor centers, the excitatory V0V CINs also provide mutual inhibition between the flexor centers (receiving inputs from excitatory V2a and acting through inhibitory Ini neurons), the excitatory V3 CINs provide mutual excitation between the flexor centers, and the inhibitory CINi CINs provide inhibition from the extensor to the contralateral flexor centers.

Connections within the spinal cord.

(A) Connections between left and right rhythm generators (RG) within each girdle. (B) Connections from the fore to hind RGs via descending (cervical-to-lumbar) long propriospinal neurons (LPNs). (C) Connections from the hind to fore RGs via ascending (lumbar-to-cervical) LPNs. (D) Brainstem drive connections to the extensor and flexor centers, commissural interneurons, and LPNs. Spheres represent neural populations. Excitatory and inhibitory connections are marked by arrowheads and circles, respectively.

https://doi.org/10.7554/eLife.31050.003
Full model schematic.

Spheres represent neural populations and lines synaptic connections. Excitatory and inhibitory connections are marked by arrowheads and circles, respectively. RG, rhythm generator.

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

Connections between the cervical and lumbar circuits are mediated through mono- and polysynaptic neural pathways. Here, we modeled them as LPN connections. Ruder et al. (2016) reported that there are excitatory and inhibitory descending LPNs, whereas ascending LPNs are only excitatory. Both descending and ascending LPNs project homolaterally and diagonally. There are more diagonal then homolateral excitatory LPNs, while the opposite is the case for inhibitory LPNs. The excitatory homolateral LPNs were identified as Shox2 and the excitatory diagonal LPNs as V0V neurons. We have hypothesized that diagonally projecting inhibitory LPNs are of V0D type.

The cervical-to-lumbar connections in our model (Figure 2B) are organized accordingly: the excitatory Shox2 neurons provide excitation from each cervical extensor to its homolateral lumbar flexor center, inhibitory LPNs of unidentified type (LPNi) provide inhibition from each cervical flexor to its homolateral lumbar flexor center, V0V LPNs provide excitation from each cervical flexor to its diagonal flexor center, and inhibitory V0D LPNs provide inhibition from each cervical flexor to its diagonal lumbar flexor center.

Lumbar-to-cervical connections (Figure 2C) mirror the excitatory cervical-to-lumbar connections but do not include inhibitory LPNs.

Finally, there were two brainstem drives (Figure 2D): one that excited flexor centers and inhibited the local V0 CINs at cervical and lumbar levels and the descending diagonal V0D LPNs, promoting left-right and diagonal alternation between the corresponding RGs, and the other that provided constant excitatory drive to the extensor centers. The drive to the flexor centers, CINs, and LPNs was varied to influence locomotor speed and gait. Specifically, the drive to flexor centers defined the locomotor frequency. Based on our suggestion, this drive also inhibited V0 CINs and diagonal V0D LPNs changing the balance of excitatory and inhibitory interactions between the corresponding RGs. The full model schematic can be seen in Figure 3.

The model was used to analyze drive- and frequency-dependent changes in the locomotor phase durations and the expression of different gaits in an intact system and following the removal of particular neuron types.

Control of locomotor frequency and gait by brainstem drive

We investigated the model performance by changing the parameter α, which defined brainstem drive to flexor centers, CINs, and LPNs (see Materials and methods). The model generated oscillations when parameter α was changed from 0 to 1.05. Within this range, an increase of α led to an increase in locomotor frequency from 2 to 12 Hz (Figure 4A and B, top diagram). This increase in frequency occurred mostly due to shortening of the extension phase, while the flexion phase remained almost constant (Figure 4B, bottom diagram), which is a characteristic property of fictive and real locomotion across species (Halbertsma, 1983; Hildebrand, 1989; Frigon and Gossard, 2009; Danner et al., 2015; Shevtsova et al., 2015). The increased frequency in the model was accompanied by sequential gait changes from walk to trot to gallop and then bound (Figure 4A) consistent with experimental observations (Figure 4C, Bellardita and Kiehn, 2015).

Model performance.

(A) Gait changes during 6 s of simulation with increase of brainstem drive (parameter α) every second (shown at the bottom). Gait changes were subject to a brief transitional period (shown as transitional colors). (B) Dependency of frequency and phase durations on α. Vertical dashed lines (labeled as x1-x4) represent selected values of α related to different gaits. Note that the blue line almost completely covers the red one (indicating increase and decrease of α, respectively) in the top diagram. (C) Frequency ranges in which gaits are observed in mice (created from data from Figure 1E of Bellardita and Kiehn, 2015). (D–G) Phase differences and step cycle compositions of four characteristic gaits produced by the model that correspond to the α-values (x1-x4) in B.

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

The expressed walking gait was a lateral-sequence walk (Figure 4A,D), in which flexion phases proceeded in the following order: hind-right, fore-right, hind-left and fore-left. This walk occurred at frequencies from 2 to 4 Hz (Figure 4B). At low frequencies, the extension phase was three to four times longer than the flexion phase duration, the RGs were in flexion sequentially, and the homolateral and diagonal phase differences were close to 0.25 and 0.75, respectively (Figure 4D).

Trot (Figure 4A,E) was characterized by diagonal synchronization, left-right alternation and fore-hind alternation. It was expressed over a large range of frequencies (from 4 to 10.5 Hz; Figure 4B). At low frequencies, with longer extension than flexion phase durations, the gait constituted a walking trot and otherwise a running trot.

Gallops are characterized by synchronization or quasi-synchronization of the hind RGs and a non-zero phase difference between the fore RGs. In the model (Figure 4A,B,F), gallop was expressed at frequencies between 9 and 11 Hz (between 9 and 10.5 Hz both gallop and trot occurred). The sequence of flexion phases at cervical or lumbar levels could be either left-right or right-left but was always the same for both fore and hind RGs. Thus, transverse but not rotary gallops were expressed. The left-right phase difference of the hind RGs was always closer to synchronization than of the fore RGs. Some of the gallops exhibited an interval between the end of the hind extension and the beginning of the fore extension, in which all RGs were in flexion at the same time (Figure 4F). This corresponds to an aerial phase, characteristic for gallops.

Bound (Figure 4A,B,G) was characterized by left-right synchronization at cervical and lumbar levels, and homolateral and diagonal alternation. The fore extension phases were directly followed by hind extension phases, which were then followed by all RGs being in flexion (corresponding to an aerial phase). Bound was expressed at frequencies between 11 and 12 Hz. The flexion phase duration was always longer than the extension phase duration (Figure 4B, bottom diagram).

The frequency ranges, phase differences, and step cycle compositions (relative phase durations) of all gaits were consistent with mouse locomotion (Figure 4B–G; Bellardita and Kiehn, 2015; see gait definitions in Materials and methods).

Walk and the transition to trot: diagonal V0D LPNs ensure stable walk

To investigate control of locomotion and gait transitions by brainstem drive, we used bifurcation diagrams reflecting normalized phase differences between oscillations in hind and fore left-right, homolateral and diagonal RGs with changing parameter α that represented the brainstem drive and was used as bifurcation parameters (see examples in Figure 5). Normalized phase differences of 0.5 correspond to alternation, whereas phase differences of 0 or 1 correspond to synchronization. Gaits were operationally defined based on these phase differences (see Table 2 in Materials and methods) and are marked in the diagrams by colored areas. Blue and red lines indicate the stable phase differences with stepwise increase and decrease of α, respectively. Any discrepancies between the red and blue lines indicate regions of bi- or multistability. For these α values two or more different stable gaits coexisted and could be expressed depending on the initial conditions. Bifurcations can be seen as abrupt changes of the stable phase differences. Figure 5A shows the bifurcation diagrams of the intact model.

Bifurcation diagrams of the intact model (A), after removal of all V0V neurons (B), of only V0V LPNs (C), and of all V0V and V0D neurons (D).

Normalized phase differences of 0.5 correspond to alternation, whereas phase differences of 0 or 1 correspond to synchronization. Blue and red lines indicate the stable phase differences with stepwise increase and decrease of the bifurcation parameter α, respectively. Colored areas indicate the expressed gait.

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

At the lowest values of brainstem drive, both left-right (local V0D and V0V CINs) and diagonal alternation (diagonal V0D LPNs) promoting pathways dominated. Together with three to four times longer extension than flexion phase durations (Figure 4B bottom diagram) and homolateral alternation (present over all frequencies), this ensured that only one RG was in flexion at any given time (Figure 4A,D and Figure 5A).

With increasing α, inhibitory brainstem drive caused the activity of the diagonally projecting V0D LPNs to decrease and the balance to shift from diagonal alternation (promoted by V0D) to diagonal synchronization (promoted by V0V). Together with a progressive decrease in the relative duration of the extension phase, this caused the flexion phases of the diagonal RGs to progressively overlap and the diagonal phase difference to move from a quarter-phase lag toward synchronization (Figure 5A). In parallel, the homolateral phase difference moved toward alternation and the left-right phase differences remained at a strict out-of-phase alternation (0.5 normalized phase difference). Together this constitutes a gradual transition from walk to trot.

The sequence of the flexion phases during walk was determined by the hypothetical CINi population that was introduced in our previous model (CINi2 in Danner et al., 2016) to ensure 0.5 normalized left-right phase difference during walk.

Trot and the transition to gallop and bound: diagonal V0V LPNs stabilize trot

The transition from trot to gallop and bound was governed by brainstem-drive-controlled speed-dependent changes in the left-right and diagonal pathways. The left-right alternation promoting local CINs (V0D and V0V) received inhibition from the brainstem that increased with speed, while the left-right synchronization promoting local CINs (V3 and CINi) did not receive brainstem inhibition (Figure 2D, Figure 3). Thus, with increasing brainstem drive (and frequency), the pathway promoting left-right synchronization became stronger than those promoting left-right alternation. This resulted in a transition from trot to gallop and then to bound (Figure 5A). The diagonally projecting V0V LPNs, promoting diagonal synchronization, acted synergistically with local V0V CINs.

Two sequential bifurcations occurred: when the gait switched from trot to gallop and then from gallop to bound (Figure 5A). The bifurcation from trot to gallop resembled a subcritical (backwards) pitchfork bifurcation in the hind (and fore) left-right phase differences. The transition exhibited hysteresis: perfect left-right alternation and quarter-phase lags (in either direction) were stable at the same α-values (blue and red lines show different stable left-right phase differences in the orange shaded area in Figure 5A). In other words, when drive was increasing (blue lines), the transition from trot to gallop occurred at a higher α-value than the transition from gallop to trot when drive was decreasing (red lines). Such a hysteresis is a common feature of gait transitions in quadrupeds (Heglund and Taylor, 1988). Stable gallops occurred because the decreasing V0V CINs and LPNs activities were accompanied by a decrease of the burst duration, while V3 and CINi CINs were active during the whole flexion or extension phase, respectively. Note, that the previous model (Danner et al., 2016) could not reproduce stable gallops; they only existed as a transitional regime. Therefore, the addition of diagonal LPNs in the present model was critical for stable gallop to exist. The bifurcation from trot to gallop occurred when V0V CINs and LPNs could not support left-right alternation and diagonal synchronization any more, but V3 and CINi CINs were not strong enough to cause perfect left-right synchronization. Therefore, stable states emerged around 0/1 normalized left-right phase differences. The deviation from synchronization was determined by the duration of the bursts of V0V CINs and LPNs. As the burst durations decreased further, the left-right phase differences transitioned towards synchronization (and the diagonal phase difference towards alternation) via a supercritical pitchfork bifurcation. This bifurcation provided the transition from gallop to bound without hysteresis (blue and red lines of the left-right phase differences in Figure 5A follow the same trajectory).

Gait expression following removal of different spinal interneuron types

Experimental observations showed that mice lacking V0V neurons lose trot, mice lacking V0V and V0D only exhibit bound (Bellardita and Kiehn, 2015), and mice lacking V2a neurons transition to left-right synchronization at lower speeds compared to intact animals (Crone et al., 2009). Thus, to further validate the model, we removed several types of genetically identified neurons and compared the model performance with experimental results.

When all V0V neurons were deleted (Figure 5B) only walk, gallop and bound were expressed; trot was selectively lost. Left-right synchronization promoting pathways became dominant at a lower α-value (lower frequency), and left-right alternation was only supported by V0D CINs. Deletion of V2a had the same effect as deletion of V0V, because V2a relayed all inputs to V0V CINs and LPNs.

Deletion of only diagonally projecting V0V LPNs (and not V0V CINs) did not result in a loss of any gait, but caused the transition from trot to gallop to occur at a lower α-value and frequency (Figure 5C). Furthermore, an additional bifurcation in the left-right and diagonal phase differences occurred during trot. At low trot-frequencies and α-values, a strict out-of-phase left-right alternation (0.5 normalized phase difference) was stable, while at medium trot-frequencies, a pair of stable states near 0.5 normalized phase difference emerged that then transitioned via a region of multistability to a pair of stable states near 0/1 normalized phase difference (expression of gallop). This highlights how local V0V CINs and LPNs synergistically stabilize trot by providing left-right alternation and diagonal synchronization.

Deletion of both V0V and V0D resulted in the removal of all left-right alternation promoting pathways and caused the loss of all alternating gaits (Figure 5D). Bound was stable over all frequencies.

Deletion of descending (cervical-to-lumbar) LPNs affects left-right coordination

Ruder et al. (2016) showed that deletion of descending LPNs influences left-right coordination: at low speeds both the fore and hind limbs were alternating, then at medium speeds the hind limbs, and at high speeds both the fore and hind limbs demonstrated disturbed coordination with both alternating and synchronized activities.

In our model, deletion of the cervical-to-lumbar LPNs resulted in the emergence of new stable states (Figure 6B). At low brainstem drive values walk remained the only stable state and with increasing drive it gradually transitioned towards trot. With further increase of α, trot was first the only stable gait (Figure 6B,C). Then, both trot and gallop were stable at the same α-values (Figure 6B,D1,D2,E1,E2) and finally bound became the single stable state (Figure 6B,F). During the gallop at low α-values, the hind RGs were close to synchronization while the fore RGs remained alternating (Figure 6B,D2). With increasing α, the fore left-right phase differences gradually transitioned towards synchronization, thus gallops at higher α-values exhibited fore and hind left-right phase differences close to synchronization (Figure 6B,E2).

Model performance after deletion of all cervical-to-lumbar long propriospinal neurons (LPNs).

Bifurcation diagrams, indicating the stable phase differences depending on the bifurcation parameter α for the intact model (A) and after deletion of all cervical-to-lumbar LPNs (B). (C–F) Extension phases of stable gaits (identified by vertical dashed lines in B labeled as x1-x4). At α={x2,x3} two different gaits were stable simultaneously (depicted in D1, D2 and E1, E2, respectively).

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

The introduced asymmetry by the deletion of pathways in only one direction (from cervical to lumbar but not from lumbar to cervical) caused bifurcations of higher codimension and resulted in several bi- and multistabilities (Figure 6B). Interestingly, the stable states closely approximate the union of the stable states of the intact model Figure 6A and the model after deletion of all diagonally projecting V0V LPNs (Figure 5C). Thus, the disturbance of the fore and hind left-right phase differences was mainly caused by the deletion of diagonally projecting V0V LPNs. The remaining homolaterally projecting LPNs maintained fore-hind alternation after deletion. Local V0V CINs maintained the stable branch of left-right alternation.

Noise causes high step-to-step variability after deletion of cervical-to-lumbar LPNs

In the above section, we described that deletion of the cervical-to-lumbar LPNs resulted in the emergence of new, additional stable states at medium and high brainstem drive values. This means that two gaits were stable at the same α-values (e.g. see x2 in Figure 6B and Figure 6D1,D2). To consider how these changes in the number of steady states affect system behavior, such as step-to-step variability, we increased the noisy currents in all populations (σNoise from 0.005 pA to 1.75 pA). In the intact model, at medium brainstem drives (related to trot), the increased noise caused phase durations and phase differences to become more variable, but both, the fore and hind RGs remained alternating over all steps (Figure 7A). After deletion of cervical-to-lumbar LPNs, the applied noise caused spontaneous transitions between the steady states: the activity of hind RGs was spontaneously switching between synchronization and alternation while the fore RGs remained alternating (see dashed box in Figure 7B).

Model performance under application of noise before (intact) and after ablation of cervical-to-lumbar long projecting propriospinal neurons (LPNs).

(A,B) Exemplary extension phases of the intact model (A) and after removal of cervical-to-lumbar LPNs (B). The dashed box in B indicates synchronization of hind RGs that transiently occurred after removal of cervical-to-lumbar LPNs. (C–E,C1–E1) Phase differences categorized into three equally sized bins (gray: appropriate for trot, yellow: quarter phase difference off, red: antiphase) before and after ablation of cervical-to-lumbar LPNs for 1000 s simulation in the model (C–E) and for experimental data pooled across animals (C1–E1) at three different speeds and α-values. C1–E1 were created from data extracted from Figure 5B,C and S5A,D of Ruder et al. (2016).

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

A similar disruption of interlimb coordination has been described in mice after deletion of cervical-to-lumbar LPNs (Ruder et al., 2016). They reported random occurrences of clusters of misaligned steps and spontaneous transitions between alternation and synchronization, which were similar to those exhibited by our model. Furthermore, the model also reproduced the speed-dependent effect of changes in the proportion of appropriate and misaligned steps: at low speeds hind and fore limb coordination remained appropriate for trot, at medium speeds hind limb coordination was disturbed, and at high speeds both fore and hind limb coordination was affected (Figure 7C–E). Thus, in the regions where both trot and gallop were stable, the applied noise caused spontaneous transitions between the stable states.

Gait changes by brainstem drive and by drive-independent inputs to CINs and LPNs

Mice usually switch gaits and locomotor speeds abruptly with few or no intermediary steps (Bellardita and Kiehn, 2015). To simulate the dynamics of gait transitions, we initialized the model with a constant α-value that defined a certain gait (and frequency) and then abruptly changed α to a value characteristic of another gait. Transitions from walk to trot, from gallop to trot and from walk to gallop and back to walk (Figure 8A–C) mirrored experimental observations closely (Figure 8D–E; Bellardita and Kiehn, 2015): they usually required only one or two cycles to stabilize the new gait and locomotor frequency.

Gait transitions with changes of brainstem drive (A–C) and independent of brainstem drive (G–I).

(A–C) Instantaneous frequency, extension phases before and after parameter α was abruptly changed (bottom trace). Black vertical lines indicate the time when α was changed and the shaded green areas indicate the transitions periods. (A) Transition from walk to trot when α was changed from 0.02 to 0.4. (B) Transition from gallop to trot when α was changed from 0.85 to 0.6. (C) Transition from walk to gallop and back to walk when α was changed from 0.02 to 0.9 and to 0.02. (D–F) Experimentally observed instantaneous speed during gait changes in mice corresponding to A–C (created from data extracted from Figure 3A–C of Bellardita and Kiehn 2015). (G–I) Instantaneous frequency and extension phases during gait changes caused by additional drives to CINs and LPNs. Bottom trace indicates time course of additional drives. (G) Additional inhibitory drive to all V0V CINs and LPNs (mI=0.0, bI=0.2) at α=0.5 and caused a transition from trot to bound. During the transitional period a gallop occurred. (H) An additional excitatory drive to local V0V CINs (mE=0.0, bE=0.05) at α=0.925 caused a transition from gallop to trot. (I) Additional excitatory drive to cervical, local V0V CINs (mE=0.0, bE=0.1) at α=0.975 caused the left-right phase difference of the fore RGs to change from almost synchronized to a quarter-phase lag during gallop.

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

Finally, gait changes could also be induced by additional inputs to CINs and LPNs without changing the common brainstem drive controlling locomotor speed (frequency). For example, additional inhibition to V0V CINs and LPNs during trot induced a transition to bound without significantly affecting the frequency, resulting in a bound at a frequency characteristic of trot (Figure 8G). Excitatory drive to the local V0V CINs during gallop induced a transition to trot, again without significantly affecting speed (Figure 8H). Excitation of the local cervical V0V CINs during gallop caused the fore left-right phase differences to move farther from synchronization and the left and right extension phases of the fore RGs to overlap less, while the hind RG coordination remained almost unchanged (Figure 8I). Thus, tonic inputs to CINs and LPNs were able to alter interlimb coordination and gait almost independently of frequency.

Discussion

Control of locomotor speed and speed-dependent gait transitions

Central control of locomotor speed and speed-dependent gait expression involves many neural, biomechanical, metabolic, environmental and behavioral factors. Most of these factors are beyond the scope of the present study. Here, we focused on the potential role of central interactions between rhythm-generating circuits controlling each limb and their regulation by brainstem drives.

In our model, two such drives were considered: constant drive to the extensor centers and variable drive to the flexor centers and various CINs and LPNs. The drive to the extensor centers put them in a tonic activity regime if the flexor centers were not active. This drive may represent a separate descending supraspinal drive to extensor circuits operating through the medial reticulospinal or vestibulospinal tracts and/or inputs from cutaneous afferents and load receptors from extensor muscles (Hiebert and Pearson, 1999; Dietz and Duysens, 2000; Bouyer and Rossignol, 2003; Rossignol et al., 2006).

The brainstem drive to the flexor centers put them into a rhythmic bursting regime and strong inhibition between flexor and extensor centers mediated by inhibitory interneurons caused rhythmic flexor-extensor alternation (Zhang et al., 2014; Britz et al., 2015; Rybak et al., 2015; Shevtsova et al., 2015; Danner et al., 2016; Shevtsova and Rybak, 2016). Progressive increase of this drive increased locomotor frequency, in accordance with studies in decerebrate cats and rats, where stimulation of the MLR induced locomotor activity that increased in frequency (speed) with increasing stimulation (Orlovskiĭ et al., 1966; Shik et al., 1966; Shik and Orlovsky, 1976; Skinner and Garcia-Rill, 1984; Grillner, 1985; Atsuta et al., 1990; Nicolopoulos-Stournaras and Iles, 2009). In many of these studies, the increase in speed was accompanied by sequential changes in gait from walk to trot and to gallop and bound. The exact mechanisms for these speed-dependent gait transitions are currently unknown, but it is reasonable to suggest that these transitions are initiated and supported by inputs and drives from different excitatory and inhibitory neurons in the medullary and pontomedullary reticular formation that project to the spinal cord and can directly affect different CINs and LPNs (Alstermark et al., 1987; Jankowska et al., 2003, 2005; Matsuyama et al., 2004; Mitchell et al., 2016; Ruder et al., 2016).

Following our previous model (Danner et al., 2016), we suggest that gait transitions are caused by brainstem drives acting on neurons coupling the RGs and affecting the balance between the activities of pathways promoting synchronization and alternation. In the previous model, only homolateral LPNs and local CINs were considered. The transition from alternating (trot) to synchronized gaits (gallop and bound) resulted from excitatory drive to local V3 CINs, promoting left-right synchronization in both cervical and lumbar levels, so that with increasing drive the left-right interactions through these pathways overcame interactions via pathways supporting left-right alternation.

In the current model, the same conceptual idea was implemented by increasing inhibitory influence of the brainstem drive to local V0 CINs (V0D and V0V) and diagonal V0D LPNs (Figure 2D and Figure 3). This increasing inhibition influenced both the walk to trot and the trot to gallop and bound transitions: the walk to trot transition occurred when diagonal V0V became stronger than diagonal V0D LPNs, and the trot to gallop and bound transition occurred when local V3 became stronger than V0V CINs and LPNs. Additionally, this inhibitory drive allowed for stable solutions for gallop, which in the previous model only occurred as transient solutions (Danner et al., 2016).

Unfortunately, there is no experimental evidence for or against brainstem excitation of V3 or inhibition of V0 CINs and LPNs. Either mechanisms could exist and operate in reality. Moreover, they could even coexist and cooperate in multiple ways. This issue can be resolved in future experiments by recording from identified CINs and LPNs and analyzing their response to MLR stimulation.

Long propriospinal neurons and the organization of interactions between rhythm-generating circuits

The left-right coordination in the present model is mediated by both the local CINs, coupling left and right RGs at cervical and lumbar levels, and the diagonal LPNs, that couple cervical and lumbar RGs. The local CIN interactions (Figure 2A and Figure 3) were implemented based on earlier experimental studies (Bellardita and Kiehn, 2015) and followed our previous computational model (Danner et al., 2016). The homolateral and diagonal LPN-mediated pathways (Figure 2B,C and Figure 3) were incorporated to be consistent with (Ruder et al., 2016). The bidirectional homolateral excitatory (from each extensor to the corresponding flexor center, mediated by Shox2 interneurons) and the descending inhibitory (between the flexor centers, mediated by LPNi interneurons) LPN pathways acted synergistically to ensure fore-hind alternation over all frequencies and gaits. The descending diagonal inhibitory LPNs (V0D) together with the local inhibitory V0D CINs ensured stable quarter-phase lags between homolateral and diagonal RGs during walk. The diagonal ascending and descending excitatory LPNs (V0V) promoted diagonal synchronization and (together with the inhibitory homolateral interactions) supported left-right alternation during trot.

The presence of only cervical-to-lumbar but not lumbar-to-cervical inhibitory LPNs (Ruder et al., 2016) and stronger lumbar-to-cervical excitatory influence than cervical-to-lumbar (Juvin et al., 2005; Brockett et al., 2013) introduced a fore-hind asymmetry. Thus, the dynamics of our model were consistent with left-right but not fore-hind interchange (Schöner et al., 1990; Golubitsky et al., 1999). The left-right symmetry (because of symmetric LPN and CIN interactions) resulted in stable normalized left-right phase differences of 0.5 (in walk and trot), or 0/1 (in bound), or pairs of steady states symmetric around 0.5 (in gallops). These symmetry properties resulted in the stability of both left and right leading gallops. The fore-hind asymmetry caused the lateral-sequence walk to be the only stable walk, which is the preferred walk of several quadrupeds, including mice, dogs, cats, and humans crawling (Patrick et al., 2009; Bellardita and Kiehn, 2015; Righetti et al., 2015; Frigon, 2017), as well as the correct extension phase sequences during gallop and bound. Afferent inputs to LPNs could potentially affect the symmetry properties and result in different stable solutions. For example, cats exhibit diagonal-sequence walk when walking on a split belt treadmill with the fore and hind limbs at a different speed (Thibaudier et al., 2013).

Role of genetically identified spinal interneurons in interlimb coordination

We used our model to simulate the effects of ablating some genetically identified spinal interneurons on the speed-dependent expression and transition of different gaits. The removal of V0V neurons abolished trot (Figure 5B) and only bound could be expressed after the removal of both V0V and V0D (Figure 5D), which is fully consistent with the experimental data of Bellardita and Kiehn (2015). Also, because in our model V2a neurons mediate input to V0V neurons (Figure 2 and Figure 3; see Crone et al., 2008; Talpalar et al., 2013; Rybak et al., 2015; Shevtsova et al., 2015), removal of V2a had the same effect as removing V0V: transition to gallop and bound at a lower frequency, consistent with Crone et al. (2009).

It is necessary to note that the role of V0V neurons in the expression and support of trot in the present model significantly differs from our previous model (Danner et al., 2016). The previous model had been developed before the new data on LPN organization (Ruder et al., 2016) became available and only included homolateral LPNs and not diagonal ones. Therefore, in that model, the role of V0V neurons in supporting trot was entirely based on their involvement in support of left-right alternation. In the present model, both V0V CINs and LPNs support trot, by providing left-right alternation (CINs) and diagonal synchronization (LPNs). The model suggests that removal of one type of these pathways (e.g. the diagonal LPNs, Figure 5C) does not fully eliminate trot, but shifts transition to gallop to a lower locomotor frequency. Transition to gallop results from the increasing inhibitory influence of the brainstem drive to V0V neurons with increasing frequency. The predicted inhibitory influence of brainstem drive (e.g. from MLR) to identified V0V neurons (both local CINs and LPNs) awaits experimental testing.

Effect of ablation of the descending cervical-to-lumbar LPN pathways

Simulating deletion of cervical-to-lumbar connections was performed for additional model validation in relation to the recent experimental data: Ruder et al. (2016) showed that deletion of these connections, did not affect interlimb coordination at low speeds, disturbed left-right coordination between hind limbs at medium speeds, exhibiting spontaneous switching to left-right synchronization, and caused the same disturbances in left-right coordination between both the fore and hind limbs at high speed. In the model, this deletion resulted in additional steady states with almost synchronous left-right activities of hind RGs at the medium drive values (frequencies) and slow migration of such additional steady states to left-right synchronization of fore RGs (Figure 6B).

Introduction of noise led to random transitions between these coexisting stable states, and resulted in spontaneous changes in phase relationships between hind RGs (including occasional synchronization) at medium frequencies and both hind and fore RGs at higher frequencies (Figure 7). This qualitatively corresponded to the experimental results on ablation of cervical-to-lumber connections (Ruder et al., 2016) and hence provided additional validation to our model.

Commissural and long propriospinal neurons can serve as main targets for supraspinal and sensory afferent signals controlling limb coordination

Although limbed animals, including quadrupeds express gaits according to their speed, there are many other factors affecting selection of the appropriate gaits, including energetic/metabolic (Hoyt and Taylor, 1981), biomechanical (Alexander and Jayes, 1983; Biewener, 1990; Farley and Taylor, 1991; Fukuoka et al., 2015), environmental and behavioral factors. For example, animals need to maintain balance when maneuvering in dynamic environments, or when chasing prey or escaping from a predator. As locomotor activities are controlled by networks within the spinal cord, signals affecting interlimb coordination to accommodate these factors can be (a) common or specific supraspinal inputs or drives, reflecting additional visual, auditory, vestibular and different goal-directed factors, and/or (b) proprioceptive and cutaneous afferent feedbacks from the periphery that informs the spinal controller about various biomechanical constraints. The supraspinal control of gait transitions has been confirmed by multiple experiments during fictive or real locomotion controlled by brainstem stimulation (Orlovskiĭ et al., 1966; Shik et al., 1966; Shik and Orlovsky, 1976; Skinner and Garcia-Rill, 1984; Grillner, 1985; Atsuta et al., 1990; Nicolopoulos-Stournaras and Iles, 2009). At the same time, there are experimental data confirming the important role that proprioceptive and exteroceptive feedback plays in the control of interlimb coordination and gait (reviewed by Pearson, 1995; Pearson, 2004; Frigon, 2017).

The locomotor gait of quadrupeds is defined by the phase relationships of the movements of the different limbs. Assuming that each limb is controlled by a separate RG, each gait can be represented as a set of phase differences between the rhythmic activities generated by the RGs. In quadrupeds, all left-right interactions within the cervical and lumbar spinal cord, including interactions between left and right RG circuits, are mediated by local CINs, and all diagonal and homolateral interactions between the cervical and lumbar circuits, including interactions between the corresponding RGs, are mediated by LPNs with descending and ascending, diagonal and homolateral projections, respectively. Therefore, the activities of CINs and LPNs located within the four sections of the cord (left and right hemicords at cervical and lumbar levels) explicitly define interactions between the four rhythm-generating circuits (or RGs; Figure 9). An integrated pattern of activity of these neurons defines interlimb coordination and gait. Multiple studies confirm that spinal CINs and LPNs receive excitatory and inhibitory inputs from supraspinal structures (Skinner et al., 1979; Skinner et al., 1980; Alstermark et al., 1987; Bannatyne et al., 2003; Jankowska et al., 2003; Jankowska et al., 2005; Jankowska et al., 2006; Matsuyama et al., 2004; Cowley et al., 2008; Cowley et al., 2010; Jankowska, 2008; Szokol et al., 2011; Zaporozhets et al., 2011; Juvin et al., 2012; Mitchell et al., 2016; Ruder et al., 2016), which could explicitly provide speed-dependent gait control as in our model. At the same time, as we showed in our simulations (Figure 8G–I), inputs to CINs and LPNs other than supraspinal drives controlling locomotor frequency can also induce gait changes independent of speed. In particular, these inputs can come from sensory afferents. This is consistent with multiple experimental data showing that CINs and LPNs receive direct and indirect input from sensory afferents (Lloyd and McIntyre, 1948; Edgley et al., 2003; Bannatyne et al., 2006; Jankowska et al., 2009; Flynn et al., 2017). In accord with this is the notion that afferent feedback from each limb projects mostly to the spinal cord section controlling the same limb (Frigon, 2017) and hence its influence on spinal circuits and RGs controlling other limbs should be mostly mediated by CINs and LPNs located in this section.

Schematic representation of control of limb coordination and gait by different local commissural interneurons (CINs) and long propriospinal neurons (LPNs).

(A) Circuits controlling four limbs. Each limb is controlled by its own rhythm generator (RG). Local CINs and homolateral and diagonal LPNs couple the RGs and coordinate their activities. Supraspinal signals and sensory feedbacks (directly or also indirectly through dorsal horn interneurons; not shown) provide inputs to the RGs, CINs and LPNs. Inputs to the RGs affect the locomotor frequency (speed), while inputs to CINs and LPNs affect interlimb coordination and gait expression. (B) A more detailed representation of local CINs and LPNs that are located in one cord section (here left cervical) and integrate the corresponding intraspinal and supraspinal inputs and sensory information (from the fore left limb) to mediate the effect of this limb’s activity on spinal circuits controlling the other three limbs.

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

Therefore, we suggest that CINs and LPNs represent the main neural targets for different local/intraspinal, supraspinal, and sensory inputs to control interlimb coordination and adjust gait to various internal and external conditions. Inputs to other neural structures, such as the RG centers, could also alter interlimb coordination by transsynaptically influencing the balance of the CIN and LPN pathways, but would also affect other parameters, such as locomotor frequency. Separate (independent of brainstem drives) signals to CINs and LPNs only (e.g. from sensory afferents) can affect limb coordination and change gait independent of speed (as shown in Figure 8G–I). Figure 9 summarizes the neural circuitry within the spinal cord responsible for interlimb coordination and possible inputs from supraspinal and peripheral afferents. For simplicity, in this schematic, the same CIN and LPN populations receive all intraspinal, supraspinal and afferent inputs. In real cord, the effects of these inputs are rather mediated by separate or overlapping populations of CINs and LPNs, which may interact with each other (e.g. via presynaptic inhibition) competing for their influence on gait expression.

Model limitations

In this study, we focused only on central interactions, without considering biomechanical constraints and afferent feedbacks from the limbs. Mechanical coupling between the limbs and body can also influence interlimb coordination and interact with the neural control (Nishikawa et al., 2007). Moreover, gaits are not only characterized by the phase differences between the limbs but also by a change in kinematics and muscle activation patterns (Bellardita and Kiehn, 2015). Investigating these issues requires simulation of limb and body biomechanics in addition to the spinal circuits, which is outside the scope of the current study. We also did not consider spinal circuits operating below the RGs, such as pattern formation networks, circuits related to muscle afferent inputs, reflex circuits including Ia, and Ib interneurons, Renshaw cells and motoneurons (see Rybak et al., 2006a, 2006b; McCrea and Rybak, 2007, 2008; Zhong et al., 2012), which play an important role in controlling the limbs. Moreover, how inputs from various supraspinal structures interact with different spinal mechanisms and afferent feedback at different spinal levels in the context of controlling interlimb coordination remains largely unknown. All the above will be the focus of our future investigations.

Materials and methods

Neuron model

The model consists of a network of interconnected populations of neurons. Each population was represented by a non-spiking, ‘activity-based’ model (Ermentrout, 1994) and described by

(1) CdV/dt=-INaP-IL-ISynE-ISynI-INoise

for flexor and extensor centers, and by

(2) CdV/dt=-IL-ISynE-ISynI-INoise

for all other populations. In Equations 1 and 2, V represents the average membrane potential, C the membrane capacitance, INaP the persistent sodium current, IL the leakage current, ISynE and ISynI excitatory and inhibitory synaptic currents, respectively, and INoise a noisy current.

The leakage current was described by

(3) IL=gL(V-EL),

where gL is the leakage conductance and EL the leakage reversal potential.

Excitatory and inhibitory synaptic currents (ISynE and ISynI, respectively) of population i were described by:

(4)ISynE,i=gSynE{j[S(wji)f(Vj)]+DE,i}(ViESynE);(5)ISynI,i=gSynI{j[S(wji)f(Vj)]+DI,i}(ViESynI);

where gSynE and gSynI are synaptic conductances, and ESynE and ESynI the reversal potentials of the excitatory and inhibitory synapses; wji is the synaptic weight from population j to i (wji>0 for excitatory connections and wji<0 for inhibitory connections); and function S is defined as

(6) S(x)={x,ifx00,otherwise.

The output function f(V) translates V into the integrated population activity or neural output and was defined by the linear piecewise function: 

(7) f(V)={0,ifV<Vthr(VVthr)/(VmaxVthr),ifVthrV<Vmax1,ifVVmax.

Excitatory DE,i and inhibitory DI,i drives to population i were modeled as a linear function of the free parameter α:

(8) D{E,I},i(α)=m{E,I},iα+b{E,I},i,

where m{E,I},i is the slope, and b{E,I},i the intercept. α represents the value of the variable brainstem drive. If not otherwise specified, m{E,I},i and b{E,I},i were set to 0 (lack of drive input to the corresponding neuron population).

The persistent sodium current in the flexor and extensor centers of the RG was described by

(9) INaP=g¯NaPmh(V-ENa),

where g¯NaP is the maximal conductance, m and h activation and inactivation variables of INaP, respectively, and ENa the sodium reversal potential. Activation of the persistent sodium current was considered instant and was modeled by 

(10) m(V)={1+exp[(V-V1/2,m)/km]}-1

and slow inactivation obeyed the following differential equation

(11) τh(V)dh/dt=h(V)-h

with

(12)h(V)={1+exp[(VV1/2,h)/kh]}1;(13)τh(V)=τ0+(τmaxτo)/cosh[(VV1/2,τ)/kτ].

In Equations 10–13, V1/2 and k represent half-voltage and slope of the corresponding variable, respectively, and τ0 and τmax are the baseline and maximum of inactivation time constant τh, respectively.

The noisy current INoise was described as an Ornstein-Uhlenbeck process

(14) dINoise/dt=-INoise/τNoise+σNoise2/τNoiseξi(t)

where τNoise is the time constant, σNoise the standard deviation or strength of the noise, and ξi(t) (in 1/s) was the population-specific normalized Gaussian noise.

Model parameters

The model parameters were adapted from our previous model (Danner et al., 2016): C=10 pF; gL=4.5 nS for RG centers and gL=2.8 nS for all other neurons; g¯NaP=4.5 nS; gSynE=gSynI=10.0 nS; EL=-62.5 mV for RG centers and EL=-60 mV for all other neurons; ENa=50 mV; ESynE=-10 mV; ESynI=-75 mV; Vthr=-50 mV; Vmax=0 mV; V1/2,m=-40 mV; V1/2,h=-45 mV; km=-6 mV; kh=4 mV; τ0=80 ms; τmax=160 ms; V1/2,τ=-35 mV; kτ=15 mV; and τNoise=10 ms.

Following drive parameters were used: mE=0.0 and bE=0.1 for the extensor centers; mE=0.1 and bE=0.0 for the flexor centers; mI=0.75 and bI=0.0 for all V0D CINs; and mI=0.15 and bI=0.0 for the local homologous V0V CINs. Thus, the extensor centers received constant drive (independent of α) and the flexor centers and V0 CINs received variable drive whose changes were proportional to α. The connection weights are listed in Table 1. Weights of connections within the RGs were adapted from our previous model (Danner et al., 2016) and other weights were selected within their operating ranges and tuned to allow gait transitions to occur at experimentally observed locomotor frequencies (Bellardita and Kiehn, 2015). To test the robustness of the model, we simultaneously varied all connection weights; each weight was multiplied by a normally distributed random number with a mean of 1 and standard deviation σp of 0.01, 0.02, 0.05, and 0.10. For each σp, 100 random models were created and bifurcation diagrams were calculated. With σp0.02 all randomized models retained all stable regimes and their sequential transitions with changes of α, at σp=0.05 17% and at σp=0.10 40% of the models lost some stable solutions (gaits such as bound or trot). Thus, the final model represents a coarse system allowing parameter variations without dramatic (qualitative) changes in behavior. To simulate deletion of particular populations, their outputs were set to 0.

Table 1
Connection weights.
https://doi.org/10.7554/eLife.31050.011
SourceTarget (wij)
Within girdle and side of the cord
RG-Fi-Ini-F (0.40), i-V0D (0.70), i-V2a (1.00), i-V3 (0.35), i-V2a-diag (0.50)
f-RG-Fi-Ini-Hom (0.70), i-V0D-diag (0.50)
RG-Ei-Ini-E (0.40), i-CINi (0.40), i-Sh2-Hom (0.50)
Ini-Fi-RG-E (–1.00)
Ini-Ei-RG-F (–0.08)
V2ai-V0V (1.00)
V2a-diagi-V0V-diag (0.90)
IniV0Vi-RG-F (–0.07)
Between left and right homologous circuits
V0Dc-RG-F (–0.07)
V0Vc-IniV0V (0.60)
V3c-RG-F (0.03)
CINic-RG-F (–0.03)
Between fore and hind homolateral circuits
f-Ini-Homh-RG-F (–0.01)
f-Sh2-Homh-RG-F (0.01)
h-Sh2-Homf-RG-F (0.125)
Between diagonal circuits
f-V0D-diagdh-RG-F (–0.075)
f-V0V-diagdh-RG-F (0.02)
h-V0V-diagdf-RG-F (0.065)
  1. i-, ipsilateral; c-, contralateral; f-, fore; h-, hind; d, diagonal. CINi, inhibitory commissural interneurons. Ini, regular inhibitory interneurons. RG-F, flexor center; RG-E, extensor center.

Computer simulations

The set of differential equations were solved using the odeint (Ahnert et al., 2011) implementation of the Runge-Kutta-Fehlberg 7–8 variable step-size solver of the boost C++ library (version 1.63.0). INoise was calculated before the simulation with the Forward Euler method and for a fixed step size of 1 ms. The custom C++ code was compiled with clang 800.0.42.1 (Apple LLVM 8.0.0) for macOS 10.12.3. Simulation results were analyzed using Matlab 2016b. Source code and Matlab scripts for all simulations are available in ModelDB (McDougal et al., 2017) at http://modeldb.yale.edu/234101.

Data analysis

Activities of the flexor and extensor center were used to assess the model behavior. A RG was considered as being in flexion when f(V) of its flexor center was greater or equal to 0.1, otherwise it was considered as being in extension. The locomotor period was defined as the duration between two consecutive onsets of the left hind flexor center; the frequency as the reciprocal of the period; extension and flexion phase durations as the duration between onset and offset of extension and flexion, respectively. Normalized phase differences were calculated as the delay between the onsets of the extension phases of a pair of rhythm generators divided by the locomotor period. Four normalized phase differences were calculated: hind left-right (right hind – left hind), fore left-right (right fore – left fore), homolateral (left fore – left hind) and diagonal (right fore – left hind). Gaits were operationally defined based on the phase differences (see Table 2).

Table 2
Operational definition of gaits.
https://doi.org/10.7554/eLife.31050.012
Normalized phase differences
GaitLeft-right hindHomolateralDiagonal
Walk*[0.25,0.75][0.1,0.4) and (0.6,0.9](0.1,0.4] and [0.6,0.9)
Trot[0.25,0.75][0.25,0.75][0.0,0.1] and [0.9,1.0)
Gallop(0.025,0.25] and [0.75,0.975)[0.25,0.75][0.25,0.75]
Bound[0.0,0.025] and [0.975,1.0)[0.25,0.75][0.25,0.75]
  1. *Classification of walk additionally required longer extension than flexion phase durations.

Analysis of model performance

To identify stable solutions of the model, bifurcation diagrams were built for the four normalized phase differences. To this end, α was increased from 0.0 to 1.05 and then decreased back to 0.0 in 1000 equally spaced steps (step size of 0.00105). At each step, the simulation was run for 10 s and the final state was used as the initial condition for the next step. Circular statistics were evaluated for the last five locomotor cycles. At each step, the circular standard deviation of the normalized phase differences of the last five steps was smaller than 0.001. Thus, the mean phase differences can be regarded as the stable solutions. The sequential increase and decrease of α was performed to uncover regions where multiple states are stable. An additional simulation run with stepwise change of α in the opposite direction was initialized when a discrete change of phase differences between two α-values occurred, to uncover potentially missed stable trajectories. Initial conditions were randomized and multiple runs were evaluated. σNoise was set to 0.005 pA, making noisy currents several orders of magnitude smaller than any other current. The noise did not affect the model behavior other than to ensure that the system did not remain on an unstable trajectory.

To model step-to-step variability as presented in decabper, simulations with increased noisy currents (σNoise=1.75 pA) were performed. The free parameter α was set to 0.3, 0.6, and 0.75. At each value, the simulation was run for 1000 s. Phase differences were calculated and partitioned into three bins (cf. Figure 7), depending on their appropriateness for trot and to ensure comparability with experimental data (Ruder et al., 2016). The counts per bin were then divided by the total number of locomotor cycles and multiplied by 100.

To model dynamics of gait transitions, extension phases and the instantaneous frequency are illustrate before and after an abrupt change in α-value or additionally introduced drives. To this end, the model was initialized with a predefined α-value and ten seconds of settling period were allowed before the abrupt parameter change. Parameters used are specified in the legend of Figure 8.

References

  1. 1
  2. 2
  3. 3
  4. 4
    Characteristics of electrically induced locomotion in rat in vitro brain stem-spinal cord preparation
    1. Y Atsuta
    2. E Garcia-Rill
    3. RD Skinner
    (1990)
    Journal of Neurophysiology 64:727–735.
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
    Reflex control of locomotion as revealed by stimulation of cutaneous afferents in spontaneously walking premammillary cats
    1. J Duysens
    (1977)
    Journal of Neurophysiology 40:737–751.
  24. 24
  25. 25
  26. 26
    Verlag Von F.C.W. Vogel
    1. L Edinger
    (1896)
    Vorlesungen über den Bau der nervösen Centralorgane des Menschen und der Thiere, Verlag Von F.C.W. Vogel, 5 edn, Leipzig, Deutschland.
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
    Control of locomotion in bipeds, tetrapods, and fish
    1. S Grillner
    (1981)
    In: J. M Brookhart, V. B Mountcastle, editors. Handbook of Physiology—The Nervous System II. American Physiological Society. John Wiley and Sons, Inc. pp. 1179–1236.
  38. 38
  39. 39
    The stride cycle of the cat: the modelling of locomotion by computerized analysis of automatic recordings
    1. JM Halbertsma
    (1983)
    Acta Physiologica Scandinavica. Supplementum 521:1–75.
  40. 40
    Speed, stride frequency and energy cost per stride: how do they change with body size and gait?
    1. NC Heglund
    2. CR Taylor
    (1988)
    The Journal of Experimental Biology 138:301–318.
  41. 41
  42. 42
  43. 43
    Contribution of sensory feedback to the generation of extensor activity during walking in the decerebrate Cat
    1. GW Hiebert
    2. KG Pearson
    (1999)
    Journal of Neurophysiology 81:758–770.
  44. 44
  45. 45
  46. 46
  47. 47
    Neuronal basis of crossed actions from the reticular formation on feline hindlimb motoneurons
    1. E Jankowska
    2. I Hammar
    3. U Slawinska
    4. K Maleszak
    5. SA Edgley
    (2003)
    Journal of Neuroscience 23:1867–1878.
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
    Analysis of forelimb-hindlimb reflex activity in acutely decapitate cats
    1. DP Lloyd
    2. AK McIntyre
    (1948)
    Journal of Neurophysiology 11:455–470.
  57. 57
  58. 58
    Mediation of descending long spinal reflex activity
    1. DPC Lloyd
    (1942)
    Journal of Neurophysiology 5:435–458.
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71
  72. 72
    [Locomotion induced by stimulation of the mesencephalon]
    1. GN Orlovskiĭ
    2. FV Severin
    3. ML Shik
    (1966)
    Doklady Akademii nauk SSSR 169:1223–1226.
  73. 73
  74. 74
    Neural control of locomotion
    1. KG Pearson
    2. J Duysens
    (1976)
    In: R. M Herman, S Grillner, P. S. G Stein, D. G Stuart, editors. Function of Segmental Reflexes in the Control of Stepping in Cockroaches and Cats. Boston: Springer. pp. 519–537.
  75. 75
  76. 76
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81
  82. 82
  83. 83
  84. 84
  85. 85
  86. 86
  87. 87
    [Control of walking and running by means of electric stimulation of the midbrain]
    1. ML Shik
    2. FV Severin
    3. GN Orlovskiĭ
    (1966)
    Biofizika 11:659–666.
  88. 88
    Neurophysiology of locomotor automatism
    1. ML Shik
    2. GN Orlovsky
    (1976)
    Physiological reviews 56:465–501.
  89. 89
  90. 90
  91. 91
  92. 92
  93. 93
  94. 94
  95. 95
  96. 96
  97. 97

Decision letter

  1. Ronald L Calabrese
    Reviewing Editor; Emory University, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Long propriospinal neurons and gait expression in quadrupeds" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom, Ronald L Calabrese (Reviewer #1), is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Eve Marder as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Stephen Soffe (Reviewer #2); Rob Brownstone (Reviewer #3).

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

Summary:

This modeling study uses the latest available data on spinal locomotor networks and behavior in mice and constructs a model to explore possible mechanisms of gait transitions. The model focuses on interneurons that interconnect limb CPGs both locally within the fore or hind limb enlargement and between the limb centers. It emphasizes the introduction of contralaterally projecting long propriospinal neurons (LPNs) as previous models from the same group did not include these recently described LPNs, but only short projecting commissural interneurons (CINs) and ipsilaterally projecting LPNs. Certain populations of LPNs and CINs are subject to changing descending excitatory drive that targets the flexor half-center of the limb CPGs. The model captures well speed-dependent gait expression based on this descending drive. The model is further constrained by testing its ability to capture changes in gait caused by genetic ablation of specific groups of spinal interneurons. In particular, it captures newly discovered disruption of hindlimb coordination following ablation of descending LPNs. They further show that inputs to CINs and LPNs can affect interlimb coordination and change gait independent of speed. The model makes the strong prediction that LPN and CINs interneurons represent the main targets for supraspinal and sensory afferent signals adjusting gait.

Spinal networks controlling locomotion are a very active area of research and this model organizes the vast amount of available data on genetic identification, connectivity, genetic ablation, and behavior and presents testable scheme for the organization of interlimb coordination. We are impressed by the model and the way it has evolved from a previous version, and found the paper a very interesting read. It is an advance on the previous published version because it incorporates new data that endows the model with new explanatory and predictive power. As such, it should be of wide interest.

The work is carefully done and the illustrations are beautiful and helpful as well as presenting the necessary data.

Essential revisions:

1) Rather than being a fundamentally new approach, this paper presents an advance in a model that has been evolving: for example recently presented in their 2016 (J. Physiol) paper. This new version of the model presents a more sophisticated connectivity based on new physiological/behavioral data, i.e. it follow the availability of more detailed evidence on homolateral and crossed pathways between fore and hind-limb centers (Ruder et al., 2016). Interestingly, the 2016 version already presented a model within which control of descending drive to specific points in the spinal network produced changes in locomotor speed and concomitant gait changes. So what does the addition of further connectivity add? It’s necessary to make clearer the relationship between the two models. Do the diagonal connections add something fundamentally different or are they there to fine tune coordination? It seems that they 'improve' some transitions between gaits (e.g. discussed briefly around paragraph four of the Discussion) but do they fill an important functional gap in the 2016 model? The Introduction should be a more thorough job of putting the model in the context of what is known in the field and what the major issues are.

Does the similarity in the capability of the 2016 and new model show a weakness of relatively unconstrained modelling if you can get similar functionality out of quite different levels of coupling? Some discussion along these lines would seem appropriate.

2) The paper provides a clear story of how coordination could work through a balance of different influences. However we are concerned about how robust the model is in terms of choices of connection weights. The authors acknowledge, quite rightly, that direct experimental evidence is seriously lacking. How were the various strengths determined and how critical are they? A further addition to the previous model (as far as we can tell) is the element of descending inhibitory drive to V0Vs and V0Ds; is there experimental evidence to support this? It would be helpful to provide a clearer picture in this paper of some of the key history and limitations of the model.

3) A major concern relates to the descending drive. The authors have justified why they have set it up as is (because it works). Whether there is a necessity for projection to the extensor centers (which receive "relatively high drive") is not clear (as it produces only tonic activity in them). It's not until the second paragraph of the Discussion that it becomes clear that there are 2 different (independent?) descending drives – clarity about this is important. Is there a physiological basis to this? It gets less clear in the third paragraph of the Discussion. Could the "drive" to extensors not be accomplished by simply an increase in excitability of the target neurons? Modulators, for example? Most importantly, though: α is poorly defined. Furthermore, at times α and locomotor frequency are used interchangeably. Is the relationship between the two the same in all conditions (the various knockouts)? As the drive forms the crux of the study, it is important to be crystal clear about it, how the values were assigned, what they mean, what is the physiological basis for these values, etc. Can Figure 2 be used to illustrate what happens when α is changed? Finally, given the relative paucity in knowledge about descending drive, explicitly stating the predictions of its nature in the discussion would be helpful. A corollary concern relates to sensitivity of the model to changes in α. For example, in Figure 8 how dependent are the transitions in G-I on the exact α value chosen? Why are 3 significant figures used to specify this parameter; is extreme precision needed?

4) Variability, noise, and Figure 7. This is very difficult to follow. It is not clear what this means in terms of physiology – where is the increasing noise coming from? And given that the noise is increased to the degree that it is, can the reader assume that there is a high degree of tolerance of these circuits to noise? Or do they normally operate near this high limit and thus need to be very well tuned to prevent susceptibility to noise? Figure 7 is difficult to interpret.

5) Afferent feedback is first mentioned early in the Introduction (in an odd sentence stating that the feedback projects to the limbs?). It then arises later in the Discussion. This topic appears tangential to the manuscript and not fully explored. The authors might want to consider leaving this for a subsequent manuscript devoted to the topic, exploring how sensory stimulation can alter circuit function. In any case, they might want to consider removing this topic from the current manuscript to gain on focus and clarity. (This would also mean Figure 9 isn't needed.) (Furthermore, even sensory input confined to fibres terminating in the dorsal horn (not included in this model) can affect gait. Presumably this would be via indirect pathways.) In fact, the authors say in subsection “Model limitations” that they aren't considering afferent input, so why are they?

6) It might be interesting from a comparative biology point of view to consider scaling in the discussion. What happens as you move from mice to cats to horses, for example? Could there be a "simple" evolutionary process (in terms of these circuits) that changes the relationship between α and locomotor frequency? And speaking of horses, there has been recent work on Dmrt3, which is expressed in a subset of spinal interneurons. Dmrt3 mutations lead to alterations in interlimb coordination (studies replicated in mice). How do these neurons (dI6?) fit into the picture? I came to this by looking at Figure 6D2, which is very close to an even, 4-beat gait as seen in the Icelandic horse tölt. What would be needed for this circuit to produce a tölt? Or to produce a pace, in which the two left limbs alternate with the two right limbs?

7) The model. How were the values in Table 1 determined? Is there a physiological basis for them, or were they derived iteratively? Either way is okay, but we should know where the numbers come from. Furthermore, how robust is the model to changes in these values? Would the degree of robustness tell us anything about the circuits or their development?

8) Several qualitative comparisons to experimental data are made. These should to the extent that is possible be made quantitative, e.g., the data of Figure 7 C-E vs. C1-E1.

9) Not all readers interested in locomotor CPGs are well versed in bifurcation diagrams. While we concur with the authors that such diagrams convey a lot of information and insight succinctly, the crux is that they must be understood. The bifurcation diagrams are not explained adequately in the text. Help lead the reader through them by explicit reference to specific curves. In particular, describe what is multistability and how it appears in the diagrams. Be clear about the hysteresis and explain how it is uncovered and visualized in the diagrams. Help the reader through these diagrams (one example thoroughly explained will do); your work deserves to be understood.

10) Materials and methods are not fully adequate. Give more details about the model, e.g., how many cells are in a population? Did all cells in a population have exactly the same parameters, inputs, noise? Beef up the sections on Analysis of model performance and Data analysis. Make the complete model and all the simulations available on ModelDB or similar site and on DRYAD or similar sites respectively. The model code etc. MUST be posted in a way that others have direct and complete access to it.

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

Author response

Essential revisions:

1) Rather than being a fundamentally new approach, this paper presents an advance in a model that has been evolving: for example recently presented in their 2016 (J. Physiol) paper. This new version of the model presents a more sophisticated connectivity based on new physiological/behavioral data, i.e. it follow the availability of more detailed evidence on homolateral and crossed pathways between fore and hind-limb centers (Ruder et al., 2016). Interestingly, the 2016 version already presented a model within which control of descending drive to specific points in the spinal network produced changes in locomotor speed and concomitant gait changes. So what does the addition of further connectivity add? It’s necessary to make clearer the relationship between the two models. Do the diagonal connections add something fundamentally different or are they there to fine tune coordination? It seems that they 'improve' some transitions between gaits (e.g. discussed briefly around paragraph four of the Discussion) but do they fill an important functional gap in the 2016 model? The Introduction should be a more thorough job of putting the model in the context of what is known in the field and what the major issues are.

Does the similarity in the capability of the 2016 and new model show a weakness of relatively unconstrained modelling if you can get similar functionality out of quite different levels of coupling? Some discussion along these lines would seem appropriate.

We substantially revised the Introduction to clarify the state of knowledge in the field in general and in connection with differences between our previous (Danner et al., 2016) and present model.

Concerning the differences between the previous (Danner et al., 2016) and the present models, there are the two major differences. The first difference concerned the organization of descending (cervical-to-lumbar) and ascending (lumbar-to-cervical) connections mediated by LPNs. The previous model had been developed before the data of Ruder et al. became available. Therefore, at that time, we just tried to find and a minimal solution for LPN pathways organization that would allow the model to reproduce the major gaits and transitions between them. The LPN connectome proposed in that model included only homolateral inhibitory interactions between the cervical and lumbar RGs and did not include any diagonal interactions. The study of Ruder et al., 2016 did not confirm the proposed connectome. Hence, our task in the present study was to find another solution that, on one hand, would allow reproducing the major gaits and their transitions, and, on the other hand, would be consistent with the current knowledge on the organization of LPN pathways, particularly with the recent Ruder et al., 2016 study. Such a solution has been found and, we believe, represents a step toward better understanding of spinal circuit organization. Incorporating diagonal LPN connections (and removal of ascending inhibitory connections) in the present model was primarily made to fit Ruder et al. data, and may not represent a fundamental difference except the fact that these revised connections allowed the existence of stable gallop, which in the previous model could only occur as a transient regime. We have highlighted the difference between the models concerning gallops in the last paragraph of the section “Trot and the transition to gallop and bound: diagonal V0V LPNs stabilize trot” of the Results.

The second difference between the previous (Danner et al., 2016) and the present model was in the organization of supraspinal inputs to CINs and LPNs provided the speed-dependent gait transitions. Importantly, in both models we have suggested that gait transitions are caused by brainstem drives acting on these neurons and changing the balance between the pathways promoting left-right synchronization and left-right alternation. In the previous model, the transition from alternating (trot) to synchronized gaits (gallop and bound) resulted from excitatory drive to local V3 CINs, promoting left-right synchronization at both cervical and lumbar levels. In the current model, the same conceptual idea was implemented by increasing inhibitory influence of the brainstem drive on local V0 CINs (V0D and V0V) and diagonal V0D LPNs involved in left-right alternation. Both implementations of the same idea can exist and operate in reality. Moreover, they can coexist and cooperate in multiple ways. Unfortunately, at the current insufficiency of experimental data, we cannot say which of these solutions is more realistic. However, we do not think that suggesting alternative solutions represents a weakness of the model; we would rather consider this as an advantage. Moreover, we now not only suggest the possible alternative versions, but also propose experiments (e.g., recording the responses of different identified CINs and LPNs to MLR stimulation) that can resolve this issue and show which of the suggested versions is more realistic and closer to reality. This issue is now additionally discussed at the very end of the section “Control of locomotor speed and speed-dependent gait transitions” of the Discussion.

2) The paper provides a clear story of how coordination could work through a balance of different influences. However we are concerned about how robust the model is in terms of choices of connection weights. The authors acknowledge, quite rightly, that direct experimental evidence is seriously lacking. How were the various strengths determined and how critical are they? A further addition to the previous model (as far as we can tell) is the element of descending inhibitory drive to V0Vs and V0Ds; is there experimental evidence to support this? It would be helpful to provide a clearer picture in this paper of some of the key history and limitations of the model.

The weights of connections were selected (within their operating ranges) to allow the changes in the balance (and hence in the corresponding gait transitions) to occur at the values of α that corresponded to the experimentally measured values of frequency at which these gait transitions occurred. Hence, relative changes of these weights can change/shift the corresponding points of gait transitions (values of α and transition frequencies). Yet, the model represents a coarse system and allows parameter variations without dramatic (qualitative) changes in behavior. This is now mentioned at the end of section “Model parameters” of Materials and methods.

Inhibition of V0 neurons by brainstem drive as a mechanism for gait transitions does not have experimental support so far. This is our suggestion and model prediction. In the previous model (Danner et al., 2016), we had a different idea that the increasing brainstem drive activates the excitatory V3 CINs. Unfortunately, there is no experimental evidence for or against brainstem excitation of V3 or inhibition of V0 CINs and LPNs. Either, or both, mechanism can operate in reality. This issue can be resolved in future experiments by recording from genetically identified CINs and LPNs and analyzing their response to MLR stimulation. We state this now in the Discussion section (last paragraph of “Control of locomotor speed and speed-dependent gait transitions”).

3) A major concern relates to the descending drive. The authors have justified why they have set it up as is (because it works). Whether there is a necessity for projection to the extensor centers (which receive "relatively high drive") is not clear (as it produces only tonic activity in them). It's not until the second paragraph of the Discussion that it becomes clear that there are 2 different (independent?) descending drives – clarity about this is important. Is there a physiological basis to this? It gets less clear in the third paragraph of the Discussion.

This comment concerns several different issues and we respond to the individual issues separately. First of all, the idea that increasing brainstem drive causes increasing locomotor frequency (speed) and speed-dependent gait transitions was based on multiple experimental data showing that MLR stimulation in decerebrate cats and rats induce locomotor activity, whose frequency increased with increasing intensity of stimulation (Orlovskii et al., 1966; Shik et al., 1966; Shik and Orlovsky, 1976; Skinner and Garcia-Rill, 1984; Grillner, 1985; Atsuta et al., 1990; Nicolopoulos-Stournaras and Iles, 2009). In many of these studies, such an increase in locomotor frequency was accompanied by sequential changes in gait from walk to trot, gallop and bound. This is described in our paper in the third paragraph of the section “Control of locomotor speed and speed-dependent gait transition” in the Discussion.

An independent high drive to the extensor centers represents a very important issue. The possible origins of this drive had been only shortly discussed in the Discussion (second paragraph; i.e., separate descending drive through reticulospinal or vestibulospinal tracts and/or inputs from cutaneous afferents and load receptors from extensor muscles). This issue is related to the flexion-dominated concept of locomotor CPG organization (Pearson and Duysens, 1976; Duysens, 1977; Zhong et al., 2012; Duysens et al., 2013; Machado et al., 2015), which was used in our previous models (Molkov et al., 2015; Shevtsova et al., 2015; Danner et al., 2016; reviewed in detail in Rybak et al., 2015, and Ausborn et al., 2017). The idea is that the extensor centers can independently generate intrinsic rhythmic bursting under certain conditions (similar to the flexor centers, e.g., see Hägglund et al., 2013). However, as in the previous models, we suggested that under normal conditions, they receive high excitation, keeping them out of intrinsic bursting mode, and hence they normally operate in a regime of tonic activity and exhibit rhythmic bursting due to the rhythmic inhibition from the intrinsically oscillating flexor centers. This was a basis for considering separate drives to flexor and extensor centers.

We agree that we have not adequately explained the two different drives. We have now included an explanation and clarification of the extensor and flexor drives with the corresponding references to the first paragraph of the section "The model" in Results and made corresponding changes in second paragraph from the end of this section.

Could the "drive" to extensors not be accomplished by simply an increase in excitability of the target neurons? Modulators, for example?

It could, and as soon as extensors are excited enough to operate in the regime of tonic activity (if isolated). Therefore, it does not matter whether their excitation is produced by a separate excitatory drive (through synaptic activation) or by modulators. However, the former scenario looks more realistic, since it is more difficult to imagine that modulators act selectively only on extensor-related RG neurons.

Most importantly, though: α is poorly defined. Furthermore, at times α and locomotor frequency are used interchangeably. Is the relationship between the two the same in all conditions (the various knockouts)? As the drive forms the crux of the study, it is important to be crystal clear about it, how the values were assigned, what they mean, what is the physiological basis for these values, etc.

In the first sentence of the section “Control of locomotor frequency and gait by brainstem drive” in the Results we stated that parameter α “defined brainstem drive to flexor centers, CINs, and LPNs” with reference to Materials and methods (where the dependence of brainstem drive(s) on α is clearly defied by equation 8 with all parameters specified in this chapter). Then, we stated “the model generated oscillations when parameter α was changed from 0 to 1.05. Within this range an increase of α led to an increase in locomotor frequency from 2 to 12 Hz (Figure 4A and B, top diagram) […].” Also, Figure 4B (top diagram) clearly shows the direct relationship between α and locomotor frequency generated by the model, including the relationships for values of α when gait changes. The frequency of oscillation in this and previous models was fully defined by the brainstem drive to flexor centers and were almost independent of the phasic interactions between RGs mediated by CINs and LPNs, which were much weaker and could only affect phase relationships between the RGs, but not the oscillation frequency or amplitude. This was previously discussed by Rybak et al., 2015 and has been now mentioned in first paragraph of the section “The Model” in Results. Because of this, the dependences of frequency and phase durations on α, shown in Figure 4B, were almost unchanged with removal of any CIN and LPN types.

The only physiological basis for the relationship between the α (conditionally representing brainstem or MLR drive) and the locomotor frequency is the data showing that MLR stimulation in decerebrate cats and rats induce locomotor activity whose frequency increased with increasing intensity of stimulation (Orlovskii et al., 1966; Shik et al., 1966; Shik and Orlovsky, 1976; Skinner and Garcia-Rill, 1984; Grillner, 1985; Atsuta et al., 1990; Nicolopoulos-Stournaras and Iles, 2009), which we mentioned in several places through the paper, but the exact dependence has never been studied experimentally, neither in intact, nor in genetically transformed animals.

Can Figure 2 be used to illustrate what happens when α is changed? Finally, given the relative paucity in knowledge about descending drive, explicitly stating the predictions of its nature in the discussion would be helpful. A corollary concern relates to sensitivity of the model to changes in α. For example, in Figure 8 how dependent are the transitions in G-I on the exact α value chosen? Why are 3 significant figures used to specify this parameter; is extreme precision needed?

Figure 2 shows (physical) connectivity in the model. We are not sure how we can illustrate what happens when α changes. We have changed the color of the brainstem drive to the extensor centers in Figures 2 and 3 to make clear that this drive is separate from the drive to flexor centers, and that the drive to the flexors and to CINs/LPNs have the same origin. The current knowledge about the brainstem drive is now described in the revised Introduction (second paragraph from the end of Introduction) and in the second paragraph of the Discussion. In the Discussion (including the revised last paragraph) we describe our predictions concerning the effects of brainstem drive to CINs and LPNs. Based on the above, any change in α (drive to the flexor centers) will result in the corresponding change in the locomotor frequency and, in some critical points, can also change gait as shown in Figure 4A,B.

The three examples shown in Figure 8G-I were chosen to show three different gait transitions. For example, Figure 8G shows a transition from trot to gallop and bound at an α value of 0.5 when additional inhibitory drives were applied to V0V CINs and LPNs. This α value was chosen because it is in the middle of the range of α where trot is expressed (see Figure 5A). The α-values in Figures 8H,I have more significant digits because gallop is expressed in a much narrower α-range than trot (see Figure 5A).

We were specifically interested in showing the dynamics gait transitions so they can be related to the examples of Bellardita and Kiehn, 2015 and to the model when α was changed. The three figures (diagrams) in Figure 8A-C were selected to reproduce similar changes as exemplified in the paper of Bellardita and Kiehn (2015, Figure 3A-C). The correspondence here was only qualitative, the exact tuning of α was not necessary.

4) Variability, noise, and Figure 7. This is very difficult to follow. It is not clear what this means in terms of physiology – where is the increasing noise coming from? And given that the noise is increased to the degree that it is, can the reader assume that there is a high degree of tolerance of these circuits to noise? Or do they normally operate near this high limit and thus need to be very well tuned to prevent susceptibility to noise? Figure 7 is difficult to interpret.

This issue requires explanation. It does not matter where this noise came from and what is the level of noise (if the latter is within a reasonable range). The noise (representing natural stochastic processes or multiple signals of unknown origin) is always present in any real system. What is important here is how many steady states (regimes) simultaneously exist (co-exist) in the system. If the system has only one steady regime, it operates close to this regime, and noise can produce a small deviation around it. However, if a nonlinear system has two (or more) steady states (regimes) that coexist (which indicates bi-/multistability), then noise allows the system to spontaneously “jump” from one steady regime to another. In such case, during some time intervals, the system operates in, or close to, one steady regime and during some other time intervals, operates in or around another steady regime(s). This is a general property of multistable dynamic systems. In this connection, the level of the noise just defines how often such jumps happen. Therefore, the noise (including natural noise or uncorrelated inputs in real system) allows distinguishing the above two situations.

In the section “Deletion of descending (cervical-to-lumbar) LPNs affects left-right coordination” of the Results, with reference to Figure 6, we described that deletion of the cervical-to-lumbar LPNs in our model resulted in the emergence of new stable states when α (and frequency) was increasing. Specifically, at low brainstem drive (α) values (or low locomotor frequencies) walk and then trot remained the only stable state, but at medium drive (α) values (medium frequencies) both steady trot and steady gallop coexisted at the same drive (α) values. These changes in steady regimes become obvious with including noise as shown in Figure 7. An increased noise resulted in spontaneous transitions between two steady regimes (seen in Figure 6 and described above): one characterized by left-right alternation (phase differences around 0.5, specific for trot, and the other characterized by left-right relationships corresponding to gallop (see Figure 7B, indicated by dashed square, and Figure 7D,E, orange and red shaded areas). These modeling results were qualitatively very similar to the experimental data of Ruder et al., 2016 and hence provide further validation of our model. We have significantly revised the whole section “Noise causes high step-to-step variability after deletion of cervical-to-lumbar LPNs”, including descriptions related to Figure 7, to make them easier to follow for the reader.

5) Afferent feedback is first mentioned early in the Introduction (in an odd sentence stating that the feedback projects to the limbs?). It then arises later in the Discussion. This topic appears tangential to the manuscript and not fully explored. The authors might want to consider leaving this for a subsequent manuscript devoted to the topic, exploring how sensory stimulation can alter circuit function. In any case, they might want to consider removing this topic from the current manuscript to gain on focus and clarity. (This would also mean Figure 9 isn't needed.) (Furthermore, even sensory input confined to fibres terminating in the dorsal horn (not included in this model) can affect gait. Presumably this would be via indirect pathways.) In fact, the authors say in subsection “Model limitations” that they aren't considering afferent input, so why are they?

We agree that the topic of afferent inputs has not been explicitly explored and this is now mentioned in the section “Model limitations” of the Discussion. Our main conclusion of the paper is that we believe that CINs and LPNs are the main target for signals that control or influence interlimb coordination. In Figure 8G–I we showed that additional signals (other than descending brainstem drive that also controls frequency) to CINs and LPNs can cause changes in interlimb coordination and gait. Such control, could originate from some supraspinal centers. However, there is extensive evidence of projections from sensory afferents to CINs and LPNs. Therefore, we suggest that sensory feedback from each limb can also contribute to limb coordination and gait expression through CINs and LPNs. Based on this, we would prefer to keep the discussion of afferent feedback and Figure 9 in the paper. We updated the text to make clear that afferent feedback has not been explicitly simulated in our model. We also clarified the relationship between this conclusion and the results in Figure 8G–I. Appropriate changes have been made in the last two paragraphs of the section “Commissural and long propriospinal neurons can serve as main targets for supraspinal and sensory afferent signals controlling limb coordination” in Discussion.

The potential role of dorsal horn interneurons in mediating sensory inputs has been mentioned in the legend of Figure 9.

6) It might be interesting from a comparative biology point of view to consider scaling in the discussion. What happens as you move from mice to cats to horses, for example? Could there be a "simple" evolutionary process (in terms of these circuits) that changes the relationship between α and locomotor frequency? And speaking of horses, there has been recent work on Dmrt3, which is expressed in a subset of spinal interneurons. Dmrt3 mutations lead to alterations in interlimb coordination (studies replicated in mice). How do these neurons (dI6?) fit into the picture? I came to this by looking at Figure 6D2, which is very close to an even, 4-beat gait as seen in the Icelandic horse tölt. What would be needed for this circuit to produce a tölt? Or to produce a pace, in which the two left limbs alternate with the two right limbs?

In this study, we tried to be careful and to use and reproduce data that are mostly related to mice. We think that scaling or extending our model to some other animals, such as horses, would be beyond the reasonable degree of speculation. Moreover, concerning locomotion, the specific differences may be connected not with evolution per se, but rather with the size, weight, kinematics, and the relationships of these characteristics to the limited maximal muscle force and other muscle characteristics. That can explain why, for example, not all results obtained in mice are reproducible in primates.

Concerning the inclusion of dI6 neurons in the model circuitry: in our previous paper (Danner et al., 2016), we suggested that the inhibitory commissural neurons CINi (see Figure 3), promoting left-right synchronization, could be a subset of dI6 neurons, but this suggestion has no support so far. In fact, we could not find a clear, specific formulation of the role of dI6 neurons in locomotion/gait control in mice, or of the effect of their selective ablation, besides the general statement that they project to both sides of the cord and “play multiple roles during locomotor activity” (Andersson et al., 2012; Dyck et al., 2012; Greiner et al., 2017). Therefore, we have been so far unable to formulate task for their modeling.

Furthermore, the presence or absence of the DMRT3 mutation influences the horse’s ability to pace (Andersson et al., 2012). However, as far as we know mice do not express pace (Bellardita and Kiehn 2015; Lemieux et al., 2015). This further complicates modeling of dI6 neurons and their involvement in the control of locomotion. Pace could potentially be reproduced in our model by including additional circuits that promote synchronization of homolateral rhythm generators (e.g., mutual excitation of the flexor centers) or alternation of diagonal rhythm generators (e.g., mutual inhibition of diagonal flexor centers) together with a mechanism that suppresses the present pathways promoting homolateral alternation and diagonal synchronization.

7) The model. How were the values in Table 1 determined? Is there a physiological basis for them, or were they derived iteratively? Either way is okay, but we should know where the numbers come from. Furthermore, how robust is the model to changes in these values? Would the degree of robustness tell us anything about the circuits or their development?

Table 1 specifies the weights of connections between neural populations in the model. The connection weights for the RG circuits were adapted from the previous model (Danner et al., 2016). All other weights have been selected within their operating ranges and tuned to reproduce the necessary behaviors of the model under considered conditions (experimentally described changes in locomotion/gait after genetic removal of particular neuron types). There is no specific (explicit) experimental data about connectivity patterns of these neurons or strength of their connections.

In the previous model (Danner et al., 2016), we tuned the parameters of RG units and their interactions to allow them to generate locomotor oscillation within the physiological range of locomotor frequencies with realistic phase durations under control of external (brainstem drive). Changes in these parameters can affect values and range of generated frequencies and phase durations. These parameters were adapted in the present model. Connection weights of the LPN and CIN pathways were tuned to make sure that gait transitions occur at the proper values of α (frequency) in order to fit corresponding experimental data. Changes in these weights can lead to gait transition at wrong frequencies or (with more dramatic changes in weights) suppress the expression of particular gaits. However, in general, the model represents acoarse system allowing parameter variations without dramatic (qualitative) changes in behavior. The process of parameter selection has been described with more details in the “Model parameters” section in Materials and methods.

8) Several qualitative comparisons to experimental data are made. These should to the extent that is possible be made quantitative, e.g., the data of Figure 7 C-E vs. C1-E1.

With the account of multiple simplifications made in the model, we believe, that only qualitative comparisons make sense here. This specifically concerns Figure 7, when our goal was only to show the speed-dependent increase of errors in limb coordination after ablation of cervical-to-lumbar LPN connections. The quantitative matching was not expected there and, in our opinion, was not really necessary.

9) Not all readers interested in locomotor CPGs are well versed in bifurcation diagrams. While we concur with the authors that such diagrams convey a lot of information and insight succinctly, the crux is that they must be understood. The bifurcation diagrams are not explained adequately in the text. Help lead the reader through them by explicit reference to specific curves. In particular, describe what is multistability and how it appears in the diagrams. Be clear about the hysteresis and explain how it is uncovered and visualized in the diagrams. Help the reader through these diagrams (one example thoroughly explained will do); your work deserves to be understood.

We thank the reviewers for this comment. We have inserted detailed explanations on how to read bifurcation diagrams in the very beginning of the section “Walk and the transition to trot: diagonal V0D LPNs ensure stable walk”. Furthermore, we now explicitly reference to the specific curves when we explain the bifurcation diagrams of the intact model.

10) Materials and methods are not fully adequate. Give more details about the model, e.g., how many cells are in a population? Did all cells in a population have exactly the same parameters, inputs, noise? Beef up the sections on Analysis of model performance and Data analysis. Make the complete model and all the simulations available on ModelDB or similar site and on DRYAD or similar sites respectively. The model code etc. MUST be posted in a way that others have direct and complete access to it.

This is probably misunderstanding. As we described in Materials and methods, “each population was represented by a non-spiking, ‘activity-based’ model (Ermentrout, 1994).” So in contrast to many our previous models (e.g. Rybak et al., 2006a,b; 2013; Zhong et al., 2012; Shevtsova et al., 2015; Shevtsova and Rybak, 2015), where we explicitly model populations of neurons described in the Hodgkin-Huxley style, we here followed our other models (Molkov et al., 2015; Danner et al., 2016), in which each population was modeled as a single unit using a non-spiking, ‘activity-based’ model (Ermentrout, 1994). The model equations are clearly described in Materials and methods (Equations 1-7) and the use of non-spiking population models had been additionally mentioned in the very beginning of the Results.

We deposited the full source code and model files on ModelDB. This code can be used to run all simulations presented in the paper. Additional Matlab scripts are provided that run these simulations, process the data, and create plots. Following line has been added to the “Computer simulations” section in the Materials and methods: “Source code and Matlab scripts for all simulations are available in ModelDB (McDougal et al., 2017) at http://modeldb.yale.edu/234101.” For now, the source code can be accessed with “1234” as a passcode (without quotation marks). ModelDB only allows publication of the model files after publication. This passcode won’t be necessary once the paper has been published.

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

Article and author information

Author details

  1. Simon M Danner

    Department of Neurobiology and Anatomy, Drexel University College of Medicine, Philadelphia, United States
    Contribution
    Conceptualization, Software, Formal analysis, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    simon.danner@gmail.com
    Competing interests
    No competing interests declared
    ORCID icon 0000-0002-4642-7064
  2. Natalia A Shevtsova

    Department of Neurobiology and Anatomy, Drexel University College of Medicine, Philadelphia, United States
    Contribution
    Investigation, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  3. Alain Frigon

    Department of Pharmacology-Physiology, Université de Sherbrooke, Sherbrooke, Canada
    Contribution
    Conceptualization, Writing—review and editing
    Competing interests
    No competing interests declared
  4. Ilya A Rybak

    Department of Neurobiology and Anatomy, Drexel University College of Medicine, Philadelphia, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition, Investigation, Visualization, Methodology, Project administration, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon 0000-0003-3461-349X

Funding

National Institutes of Health (R01 NS081713)

  • Ilya A Rybak

National Institutes of Health (R01 NS090919)

  • Ilya A Rybak

National Institutes of Health (R01 NS095366)

  • Natalia A Shevtsova

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Reviewing Editor

  1. Ronald L Calabrese, Reviewing Editor, Emory University, United States

Publication history

  1. Received: August 5, 2017
  2. Accepted: November 21, 2017
  3. Accepted Manuscript published: November 22, 2017 (version 1)
  4. Version of Record published: December 12, 2017 (version 2)

Copyright

© 2017, Danner et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 931
    Page views
  • 168
    Downloads
  • 5
    Citations

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

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)