To model and computationally investigate the brainstem control of locomotion, we built upon our previous model of spinal circuits, consisting of four RGs (each controlling one limb), which interact via local cervical and lumbar CINs and LPNs connecting cervical and lumbar compartments (Danner et al., 2017).

In our previous model, the control of locomotor frequency and gaits was provided by changes of tonic excitatory inputs to the RGs and inhibitory inputs to the particular spinal CIN and LPN populations. These inputs were considered as external brainstem drives. The implemented network architecture allowed the model to reproduce multiple experimental data, including a monotonic increase of locomotor frequency and frequency-dependent sequential gait transitions from walk to trot and then to gallop and bound (Shik et al., 1966; Shik and Orlovsky, 1976; Bellardita and Kiehn, 2015), when external drives progressively increased (Danner et al., 2017). Moreover, the model was able to reproduce changes in frequency-dependent gait expression in mutant mice lacking specific genetically identified CINs and LPNs (Bellardita and Kiehn, 2015). Here we have extended this model by incorporating brainstem compartments that include bilaterally located MLR and RF structures providing descending brainstem drives to the spinal cord.

Figure 1 shows a simplified schematic of the extended model and illustrates our main assumptions concerning circuit organization in the brainstem providing inputs to spinal circuits controlling locomotion. The full schematic of the model is shown in Figure 2. The MLR on each side includes CnF and PPN structures and the RF is represented by the LPGi. The CnF, PPN and LPGi include excitatory glutamatergic (Glu) and inhibitory GABAergic or glycinergic (GABA/Gly) neurons. Synaptic interactions within and between the CnF, PPN and LPGi on each side and their bilateral connections are organized according to the existing direct and indirect experimental data. Glutamatergic neurons in the CnF project to the contralateral CnF and to the ipsi- and contralateral PPN and LPGi (Figures 1 and 2; Caggiano et al., 2018). Glutamatergic neurons in the PPN project to the ipsi- and contralateral LPGi neurons (Figures 1 and 2; Caggiano et al., 2018). Inhibitory neurons in the CnF and PPN have local projections (Caggiano et al., 2018). LPGi inhibitory neurons project locally as well as to the ipsi- and contralateral spinal circuits in the cervical and lumbar compartments (Figure 2; Capelli et al., 2017).

Figure 1

Download asset Open asset

Figure 2

Download asset Open asset

In the present model, we hypothesized the existence of two distinct populations of glutamatergic neurons in each LPGi (LPGi-Glu-1 and LPGi-Glu-2, labeled in Figures 1 and 2 as ‘1’ and ‘2’, respectively, and, correspondingly, two distinct pathways from these populations to the spinal cord circuits. The LPGi-Glu-1 populations receive bilateral excitation from the glutamatergic neurons of both CnFs and PPNs and project to the flexor centers of the homolateral (via excitatory interneuron populations Ine) and the contralateral (via the descending commissural populations CINe) RGs in both cervical and lumbar compartments (Figure 2). This allows these LPGi populations and their descending pathways to control the frequency of locomotor oscillations generated by the RGs.

In contrast, the LPGi-Glu-2 populations receive bilateral excitatory inputs only from the CnFs and project to descending spinal interneurons (dIni) inhibiting the cervical and lumbar V0_V and V0_D CINs and the descending V0_D LPNs (Figure 2). This allows these LPGi populations and corresponding pathways to influence limb coordination and gait.

To keep the extensor centers in a tonic mode (Zhong et al., 2012; Shevtsova et al., 2015; Danner et al., 2016; Shevtsova and Rybak, 2016; Danner et al., 2017; Ausborn et al., 2018) a constant input drive (D_VN,ex = 2.15) was applied to the left and right brainstem VN populations which project to all four RG extensor centers (E, Figure 2, see also Danner et al., 2016). The VN populations in our model represent vestibular nuclei as well as other potential sources of excitatory inputs to the extensor centers of the central pattern generator involved in postural control. Vestibulospinal neurons are known to be tonically active (Orlovsky, 1972), preferentially project to spinal interneurons and slow motoneurons of the extensor pools (Miller et al., 1975; Lemon, 2008; Basaldella et al., 2015) and have been implicated in mediating extensor tone (Fulton et al., 1930; Burke et al., 1972). However, other supraspinal centers might also be involved. All other lumbar and cervical circuits and their local, ascending and descending interconnections shown in Figure 2 are described in Danner et al. (2017). Connection weights are listed in Table 1.

In their recent studies, Caggiano et al. (2018) and Josset et al. (2018) have explored the anatomical and molecular heterogeneity of the MLR, highlighting the differential role of glutamatergic neurons within the CnF and PPN. Our model was implemented and adjusted to reproduce their main results:

1. Unilateral selective activation of glutamatergic CnF neurons in the mouse can drive the full range of locomotor speeds with speed-dependent gait expression, including walk, trot, gallop and bound (Caggiano et al., 2018; Josset et al., 2018). In our model, progressively increasing the excitatory tonic drive (α, see Materials and methods) to the glutamatergic population of the left CnF resulted in an increase in locomotor frequency (upper panel in Figure 3A) up to 10–11 Hz. This increase of frequency occurred because the applied stimulation produced a progressive activation of the glutamatergic population in the contralateral CnF and in both PPNs, and finally a progressive activation of both LPGi-Glu-1 populations. The LPGi-Glu-1 populations, via the intermediate left and right Ine and CINe populations, progressively activated the flexor centers (F) of all four RGs (Figure 2) causing the observed monotonic increase of locomotor frequency. Simultaneously, stimulation of the glutamatergic population of the left CnF resulted in a progressive activation of both LPGi-Glu-2 populations, which, via both descending inhibitory dIni, progressively inhibited cervical and lumbar V0_V and V0_D CINs and descending V0_D LPNs, leading to a stimulation-dependent sequential gait transition from walk to trot and then to gallop and bound. This is illustrated in Figure 3A where left-right, homolateral and diagonal phase differences between the four RGs characterize the different gaits. Figure 3—figure supplement 1 shows examples of extensor activity and polar diagrams for each gait. For definitions and a detailed explanation of gait characteristics see Danner et al. (2016) and Danner et al. (2017). Since the structure of connections within the cord remained identical to our previous study (Danner et al., 2017), both the sequence of the expressed gaits and the intermediate regimes of bistability in Figure 3A reproduce our previous results.

2. Unilateral selective activation of glutamatergic neurons in the dorsal PPN in the mouse can only trigger low speed locomotion and generate alternating gaits (walk and trot), even if stimulated at high intensities (Caggiano et al., 2018). Our model reproduced this experimental observation since glutamatergic PPN neurons, exciting only the LPGi-Glu-1 populations, had relatively weak synaptic influence on these populations not allowing the generation of high locomotor frequencies (see Figure 3B, top diagram). Moreover, the system could not switch from alternating (walk and trot) to synchronous (gallop and bound) gaits, since PPN neurons did not excite the LPGi-Glu-2 populations and, correspondingly, did not activate the transitions to synchronous gaits via the LPGi-Glu-2-spinal cord pathways projecting to spinal CINs and LPNs (see Figure 3B).

3. Unilateral selective activation of glutamatergic CnF neurons while PPN activity is suppressed bilaterally in the mouse leads to a decrease in maximum locomotor frequency with all gaits present (Caggiano et al., 2018). Corresponding to these experimental results, in our model, unilateral activation of CnF glutamatergic neurons generated locomotor oscillations with all gaits expressed when PPN neurons were inactivated bilaterally but with slightly reduced maximum frequencies (Figure 3C). Thus, both CnF and PPN controlled low-frequency alternating-gait locomotion in the model, but only activation of the CnF (with or without the PPN) produced high frequency and left-right synchronous gaits like gallop and bound.

Figure 3 with 1 supplement see all

Download asset Open asset

While stimulation of the MLR can initiate and support locomotion of different frequencies and gaits, it does so - in the biological system (Capelli et al., 2017) as well as in our model - through activation of glutamatergic neurons in the RF, whose axons descend bilaterally to spinal circuits and presumably activate spinal RGs. Specifically, Capelli et al. (2017) have recently shown that unilateral selective stimulation of glutamatergic neurons of the LPGi, a structure within the medullary RF, can initiate locomotion and elicit locomotor oscillations and gaits in a range corresponding to that of selective CnF stimulation. Our model was able to reproduce these findings as well. Unilateral stimulation of glutamatergic LPGi populations in the model produced a similarly wide range of locomotor frequencies and corresponding gaits as was the case for unilateral stimulation of glutamatergic neurons in the CnF (Figure 3D). This suggests that the RF, and particularly the LPGi, is involved in mediating MLR control of locomotion.

Together, the above simulations have demonstrated that the proposed brainstem-spinal cord connectome allowed our model to reproduce the experimentally observed effects of stimulation of glutamatergic populations within the CnF, PPN, and LPGi.

To explicitly examine the dependence of distinct locomotor gaits on locomotor frequency, the bifurcation diagrams of Figure 3A and C were rebuilt to plot changes in phase difference against locomotor frequency (Figure 4A). We then compared locomotor gaits when glutamatergic neurons of the left CnF where stimulated (Figure 4A top diagram, same simulation as in Figure 3A) with the same stimulation while the PPN was bilaterally inactivated (Figure 4A, bottom diagram, same simulation as Figure 3C). The inactivation of the PPN shifted the transition from alternating gaits (walk and trot) to synchronized gaits (gallop and bound) to lower locomotor frequencies.

Figure 4

Download asset Open asset

The mechanism of this shift in our model is the following. As described above, the LPGi-Glu-1 pathways control locomotor frequency and the LPGi-Glu-2 pathways control gait by promoting the transition from alternating to synchronous gaits (Figures 1 and 2). When the PPN was bilaterally inactivated, only the drives to LPGi-Glu-1 were reduced since the PPN on each side only projects to LPGi-Glu-1 while the drive to LPGi-Glu-2 (promoting the transition to synchronous gaits) on each side remained unaffected (Figures 1 and 2). Thus, frequency increased more slowly with increasing CnF stimulation and the transition to synchronous gaits (gallop and bound) occurred at lower locomotor frequencies (Figure 4A,B). This shift was even more pronounced after increasing noisy current in all neurons (${\displaystyle {\sigma }_{Noise}}$ increased from 0.005 pA to 1.0 pA; Figure 4C). Incorporating a moderate noise allowed us to reproduce a natural step-to-step variability and variable frequency-dependent changes similar to those during natural locomotion. Phase differences were evaluated for each step cycle and plotted in equally spaced bins between 0 and 1 over the corresponding locomotor frequency (partitioned in 0.25 Hz bins). All values were then normalized to the total number of bins and can be interpreted as the relative frequency of occurrence of each phase difference-frequency pair (see Materials and methods). Figure 4C shows that inactivation of the PPN, while the left CnF glutamatergic population was stimulated, resulted not only in a shift of the transition from alternating to synchronized gaits to a lower frequency (Figure 4A,B) but also affected the stability of the steady-state gaits. This can be seen by an increased variability of the left-right phase differences around their stable-state solutions when the PPN is inactivated. The variability increased with increasing frequency and resulted in a wider transition period between trot (left-right alternation) and gallop (left-right synchronization).

Incorporation of step-to-step variability, as described at the end of the previous section, also allows the analysis of variable gait expression as observed in natural locomotion. To approximate frequency-dependent gait expression under noise conditions, we calculated frequency-dependent relative probabilities of expression of each gait (see Materials and methods) for three cases: when stimulation was applied to either CnF, PPN, or LPGi glutamatergic neurons (Figure 5). Activation of glutamatergic neurons in the CnF with fixed noise (${\displaystyle {\sigma }_{Noise}}$ = 1 pA) and varying drives (α) produced a wide range of frequencies and expression of left-right alternating (walk and trot) as well as left-right synchronous gaits (gallop and bound), while activation of glutamatergic neurons in the PPN resulted in lower maximum frequencies and only expression of left-right alternating gaits (Caggiano et al., 2018). While gait distributions were not reported for LPGi stimulations our model suggests that LPGi stimulation can induce both alternating and synchronized gaits (Figure 5C). Those gait transitions are expected given the relationship between gait and speed.

Figure 5

Download asset Open asset

Unilateral activation of inhibitory neurons in the CnF, PPN, or LPGi decelerates or stops ongoing locomotor activity (Capelli et al., 2017; Caggiano et al., 2018). In freely behaving mice, optogenetic activation of CnF inhibitory neurons (defined by expression of the Vgat promoter) reduces locomotor speed and, in some trials, halts locomotion completely, while activation of PPN inhibitory neurons reduces locomotor frequency to a much lesser extent with only occasional stopping of locomotion (Extended Data Figure 3 in Caggiano et al., 2018). The optogenetic activation of LPGi inhibitory neurons, like that of CnF inhibitory neurons, reduces locomotor speed down to a complete halt of locomotion at higher stimulation intensities (Figure 2g in Capelli et al., 2017).

These experiments were performed in freely walking mice, thus to simulate the ongoing locomotor activity in the model triggering the whole movement behavior from slow to high speed, we bilaterally activated CnF glutamatergic neurons by applying a constant excitatory tonic drive (D_CnF,ex) to these neurons. In each series of simulations, the inhibitory neurons in one region of interest (CnF-GABA/Gly, PPN-GABA/Gly, or LPGi-GABA/Gly) were stimulated unilaterally by applying monotonically increasing excitatory drive to the neurons (D_{CnF/PPN/LPGi,in}) within the corresponding region. The effects of these stimulations depended on the stimulated region and the initial frequency of the locomotor activity defined by D_CnF,ex.

Progressive unilateral activation of inhibitory neurons within the CnF (by increasing D_CnF,in) resulted in a reduction of locomotor frequency and an orderly progression from bound to gallop, trot and walk, and finally stopped locomotion (Figure 6A and examples for D_CnF,ex = 3.04 in Figure 6D).

Figure 6

Download asset Open asset

The same unilateral progressive activation of inhibitory population within the PPN could only decrease locomotor frequency and was not able to stop locomotion (Figure 6B and E).

Both these simulations were qualitatively consistent with the experimental data of Caggiano et al. (2018), and systematically demonstrate the possible role of inhibitory neurons in the CnF and PPN in the regulation of locomotor frequency and gait expression.

Unilateral stimulation of the LPGi inhibitory population had a similar effect on locomotor activity to that of stimulation of the CnF inhibitory neurons (Figure 6C). Moreover, similar to the experimental studies of Capelli et al. (2017) (their Figure 2f,g), progressive activation of glutamatergic populations in the CnF (increase in D_CnF,ex at fixed D_LPGi,in) in our model increased locomotor frequency, whereas progressive activation of the inhibitory population (increase in D_LPGi,in at fixed D_CnF,ex) decreased locomotor frequency up to termination of locomotor oscillations. Also, as is the case of activation of the inhibitory population in the CnF, our model predicts that the decrease of frequency with progressive activation of the inhibitory LPGi population is accompanied by orderly gait transitions (Figure 6F).

Importantly, despite similarities of the effects, the underlying mechanisms for frequency reduction with activation of LPGi inhibitory neurons was different from that in the CnF. While activation of inhibitory neurons in the CnF reduced CnF and PPN activity locally within the MLR, LPGi inhibitory neurons suppressed the activity of their downstream targets in the spinal cord, in addition to local glutamatergic neurons in the LPGi (see Figures 1 and 2).

In our previous modeling studies (Danner et al., 2016; Danner et al., 2017) we suggested that locomotor speed and gait (limb coordination) are controlled by descending brainstem drives to different targets within the spinal cord. Specifically, locomotor speed (which is dependent on the frequency of the locomotor rhythm) is determined by descending brainstem drives to the locomotor rhythm generating circuits (RGs) controlling the limbs, whereas the phase relationships between these RGs (defining inter-limb coordination and thus gait) is controlled by descending brainstem drives to the specific CIN and LPN populations mediating left-right and fore-hind interactions between the RGs. However, the previous models did not include brainstem centers and the brainstem drives were simply introduced as external inputs. In the present model, we explicitly simulated brainstem circuits and their descending pathways to the spinal cord.

To simulate MLR-controlled locomotion, the brainstem model included two bilaterally located major MLR nuclei (CnF and PPN) whose outputs converged at the left and right LPGi nuclei in the RF, which in turn projected to spinal locomotor circuits (Liang et al., 2011; Liang et al., 2016; Capelli et al., 2017; Brownstone and Chopek, 2018; Caggiano et al., 2018; Josset et al., 2018; Oueghlani et al., 2018). The decision to include only the LPGi nucleus of the RF into our model was based on results by Capelli et al. (2017), who had been able to initiate locomotion by optogenetic activation of glutamatergic neurons only in this nucleus of the RF and reported that the LPGi receives projections from both CnF and PPN. However, locomotion could also be evoked via other RF-mediated pathways not involving the LPGi (reviewed in Brownstone and Chopek, 2018). Indeed, when Capelli et al. (2017) removed LPGi glutamatergic neurons, slow locomotion could still be observed, suggesting that other nuclei can also be involved in mediating MLR pathways to the spinal cord. Hence, what in the model is defined as LPGi may also include other subpopulations of glutamatergic neurons within or even outside the RF. The identity and location of these neurons are not known at this moment.

Finally, based on our simulations we suggest that each (left and right) LPGi has separate glutamatergic populations that give rise to two separate pathways controlling locomotor frequency through activation of spinal rhythm generating circuits and gait via regulation of specific spinal CIN and LPN populations (Figures 1 and 2). We also suggest that the activities of these LPGi populations are mediated and distributed within the spinal cord bilaterally and between cervical and lumbar circuits by ipsi- and contralaterally projecting populations of interneurons (such as the Ine, CINe, and dIni populations in Figure 2). Such spinal interneurons that receive descending inputs from the RF and distribute their activity widely within the spinal cord have been found in cats (Jankowska et al., 2003; Matsuyama et al., 2004) and rats (Mitchell et al., 2016). These suggestions await experimental testing in the future.

Recent studies in mice raised a possibility that slow exploratory-type locomotion and fast escape-like locomotion might be initiated and controlled by distinct brainstem circuits. Specifically, it has been shown that slow locomotion can be evoked by activation of glutamatergic neurons in the PPN (Caggiano et al., 2018), whereas fast locomotion is initiated and controlled by glutamatergic neurons in the CnF (Caggiano et al., 2018; Josset et al., 2018). The role of the PPN in the generation (Caggiano et al., 2018) and control (Josset et al., 2018) of slow locomotion appears to be more complex. In contrast to Caggiano et al. (2018), Josset et al. (2018) conclude that activation of PPN glutamatergic neurons in resting animals can only activate flexor muscles but is not very effective in initiating locomotion per se. Yet, similar to Caggiano et al. (2018), the suppression of the PPN in Josset et al. (2018) also reduces the locomotor speed.

To simulate and propose a mechanistic explanation for the functional difference between the activation of the PPN and CnF we incorporated bilaterally interacting CnF and PPN nuclei (Figures 1 and 2). In our model, in addition to ipsi- and contralateral projections from the glutamatergic CnF populations to the glutamatergic PPN populations, both CnF and PPN glutamatergic populations project bilaterally to the left and right LPGi nuclei. However, we suggest that the projections from the PPN target only LPGi-Glu-1 populations (responsible for control of locomotor frequency, but not gait), whereas the projections from the CnF affect both LPGi populations (LPGi-Glu-1 and LPGi-Glu-2) and hence can control locomotor frequency and provide frequency-dependent control of locomotor gait (Figures 1 and 2). In addition, we assume that synaptic weights of excitatory projections from the CnF to LPGi-Glu-1 are much stronger than those from the PPN. Therefore, in our model, a unilateral progressive activation of the glutamatergic population within the CnF results in monotonic increase of locomotor frequency up to 10–11 Hz accompanied by frequency-dependent gait transitions from walk to trot to gallop and bound (Figures 3A and 5A). These results are consistent with experimental studies of Caggiano et al. (2018) (their Figure 1h–j) and Josset et al. (2018) (their Figure 6c1), providing indirect validation of the network architecture implemented in our model. In contrast, a unilateral progressive activation of the glutamatergic population within the PPN can only produce an increase of frequency up to ~8 Hz accompanied by only left-right asymmetric gaits, such as walk and trot (Figures 3B and 5B). These results are in agreement with experimental studies of Caggiano et al. (2018) (their Figure 1h–j) and support the observations of Josset et al. (2018) (their Figure 6c2 and e4) of a maintenance of asymmetrical walking at subthreshold or threshold PPN activations.

A unilateral progressive stimulation of glutamatergic populations in the LPGi produces effects - similar to those during unilateral monotonic activation of the glutamatergic population in the CnF - a monotonic increase of locomotor frequency and frequency-dependent gait transitions (Figures 3D and 5C). While the increase of frequency during unilateral monotonic activation of glutamatergic neurons in the LPGi is consistent with the results of Capelli et al. (2017), the authors of this study did not analyze gaits during their experiments, but our results suggest sequential gait changes with increasing frequencies similar to that of CnF stimulation.

Inactivation of glutamatergic neurons in the PPN leads to an overall reduction of locomotor speed (Caggiano et al., 2018; Josset et al., 2018). More specifically, Caggiano et al. (2018) studied the effect of pharmacological inactivation of the PPN on locomotor characteristics. They found that inactivation of the PPN led to a reduction of locomotor speed, but the animals were still able to express the full spectrum of locomotor gaits. Our simulations were qualitatively consistent with this conclusion (see Figure 3C). Moreover, our simulations have shown that the suppression of the PPN shifts the transition from alternating to synchronized gaits towards lower locomotor frequencies (Figure 4A,B). This shift is even more pronounced after adding noise to the model (Figure 4C).

It has been shown that unilateral optogenetic stimulation of Vgat inhibitory neurons within the CnF in freely moving mice reduces locomotor speed and can fully stop locomotion, while the same stimulation of inhibitory neurons in the PPN affects locomotor speed to a lesser extent and usually does not stop locomotion (Caggiano et al., 2018, their Extended Data Figure 3a–c). The fact that optogenetic activation of inhibitory neurons in the CnF can fully suppress locomotor activity suggests that CnF inhibitory neurons inhibit the PPN, which was implemented in our model (Figures 1 and 2). These inhibitory connections, however have not been shown experimentally. Optogenetic activation of inhibitory neurons in the RF including the LPGi progressively decreases locomotor speed and can fully stop locomotion (Capelli et al., 2017, their Figure 2f,g and Extended Data Figure 6d,h).

Our simulations were qualitatively consistent with these experimental findings (Figure 6). Notably, in our simulations, the decrease of locomotor frequency during unilateral activation of inhibitory populations in the CnF, PPN and LPGi was always accompanied by orderly progressive changes of locomotor gaits toward slower gaits like trot and walk (Figure 6D–F). Although Caggiano et al. (2018) and Capelli et al. (2017) did not analyze gait transitions during progressive increase of stimulation intensity to inhibitory neurons in these areas, our findings suggest that frequency/gait relations are not affected by unilateral activation of inhibitory neurons in these areas.

In this study we aimed to reproduce new experimental data on the brainstem control of locomotion and the specific role of CnF and PPN nuclei and the RF (Capelli et al., 2017; Caggiano et al., 2018; Josset et al., 2018) within the framework of our previous model (Danner et al., 2017). Correspondingly, we incorporated MLR and RF circuits into the model with only minimal changes to the original network structure. This ensures that our previous model assumptions and conclusions hold true for the extended model as well, but at the same time increases the complexity of the model to a point where thorough mathematical and performance analyses are not feasible. Therefore, in future studies it will be beneficial, in parallel with detail models, to generate related, but simplified network models allowing systems level mathematical analysis and more thorough investigation of the dynamic mechanisms underlying the brainstem control of locomotor frequency and gait. An example of such simplifications and systems level analyses has been performed by others for our previous model (Lodi et al., 2017; Lodi et al., 2018),

Similar to our previous modeling investigations (Danner et al., 2016; Danner et al., 2017), in this study we focused only on central neural interactions, without considering other brainstem circuits beyond MLR and RF and circuits in the spinal cord operating below the RGs, such as reflex circuits and motoneurons. Neither have we developed a full biomechanical model to investigate the role of biomechanics and sensory feedback in the control of locomotor speed and gait. We also did not take into account heterogeneous activity profiles of individual neurons (Caggiano et al., 2018; Oueghlani et al., 2018) as populations are represented by an activity based model. All the above will be the focus of future investigations.

The model represents a bilateral network of interconnected populations of neurons and includes the simulated brainstem and spinal cord compartments. The present model was based on, and represents an extension of, our previous model (Danner et al., 2017). While keeping the same spinal cord circuitry we added brainstem compartments including bilaterally interacting CnF and PPN compartments as well as LPGi compartments mediating descending brainstem drives to the spinal cord (Figures 1 and 2). Several additional relay neuron populations (Ine, CINe, and dIni) were incorporated to mediate and distribute brainstem signals to the spinal cord bilaterally and between the cervical and lumbar compartments of the cord (Figure 2, more detailed description of new network structures in Results).

Each population in the model was represented by a non-spiking, ‘activity-based’ model (Ermentrout, 1994). The flexor and extensor RG centers (F and E populations, see Figure 2) incorporated a persistent sodium current and had intrinsic oscillating properties. The average membrane potential, V, in these populations was described as:

In all other populations, the average membrane potential obeyed the following equation:

where C is the membrane capacitance, I_NaP the persistent sodium current, I_L the leak current, I_SynE and I_SynI excitatory and inhibitory synaptic currents, respectively, and I_Noise a noisy current. The output function f(V) translates V into the integrated population activity representing population output as defined by the linear piecewise function:

A complete description of the population model and parameters are presented in Danner et al. (2017). In all brainstem (CnF, PPN, and LPGi) and relay (Ine, CINe, and Ini) neuron populations the conductance variable, g_L, was equal to 5 nS.

Synaptic connection weights were adapted from our previous model (Danner et al., 2017) and weights for newly introduced connections were selected within their operating ranges and tuned to produce gait transitions similar to those of our model from Danner et al. (2017) and to reproduce the relevant experimental data (Capelli et al., 2017; Caggiano et al., 2018). Connection weights are listed in Table 1.

To simulate the effect of activation of a neuron population i (i $\in$ [CnF, PPN, LPGi, VN]) in the brainstem, we applied a tonic excitatory drive, D_i,j to this population. The stimulation strength D_i,j was given by the following equation:

where m_i is the slope and b_i the intercept. The scalar α $\in$ [0, 1.05] characterizes the variable stimulation strength. The index j (j $\in$ [ex, in]) indicates if excitatory (ex, Glu) or inhibitory (in, GABA/Gly) populations are stimulated.

The set of differential equations was solved with the same custom C++ code using odeint of the boost library used in Danner et al. (2017). The C++ code was compiled as a python module and python 3.6 was used to interface with the simulation and to analyze the results. Source code and python scripts to create all simulations presented here are available on GitHub at https://github.com/SimonDanner/CPGNetworkSimulator (Danner, 2019; copy archived at https://github.com/elifesciences-publications/CPGNetworkSimulator). Data analysis procedures are described in Danner et al. (2017).

Similar to our previous models (Shevtsova et al., 2014; Shevtsova et al., 2015; Danner et al., 2016; Danner et al., 2017), the extensor RG centers in the current model were in a tonic mode while the flexor RG centers were oscillating.

To produce locomotor activity by activation of glutamatergic populations within the CnF, PPN, or LPGi, the excitatory drive was unilaterally applied to the population of interest. For each of these simulations, values for m_i and b_i (Equation 4) are indicated in the corresponding figure legends. To simulate inactivation of the PPN region, all weights of connections originating from this region were set to 0.

The bifurcation diagrams (Figure 3) were built for four normalized phase differences [hind left-right, fore left-right, homolateral (left fore – left hind) and diagonal (right fore – left hind)]. To this end, α was increased from 0.0 to 1.05 and then decreased back to 0.0 in 1000 equally spaced steps (for details see Danner et al., 2017). At each step, simulations were performed in 10 s intervals until the standard deviation of each phase-difference measured over five locomotor cycles was less than 0.001 or 200 s passed. The frequency of oscillations as the reciprocal of the period was calculated and the hind left-right bifurcation and phase-transition diagrams were built versus frequency to compare model performance in the cases when the CnF was activated with and without PPN inactivation (Figure 4).

To consider step-by-step variability (Figure 4C) and relative probabilities of frequency-dependent gait expression (Figure 5), simulations were performed with increased noisy currents (σ_Noise = 1 pA, see Equation (14) in Danner et al., 2017). To this end, the free parameter α was increased from 0.0 to 1.05 in steps of 0.01. At each step, the simulation was run for 100 s. For the left-right hindlimb phase difference and for the gait, bivariate histograms were created with cycle frequency as the second variable. The phase difference was partitioned into 65 equally spaced bins between 0 and 1 and the frequency was partitioned into 0.25 Hz wide bins from 0 to 14 Hz (Figure 4C). Gaits were evaluated at each step cycle based on the definition in Table 2 of Danner et al. (2017) and gallop and bound were grouped together (Figure 5). The counts per 2D-bin were then divided by the total number of locomotor cycles. Thus, these numbers represent the relative frequency of occurrence of each phase difference-frequency or gait-frequency pair and can be interpreted as a probability.

To simulate the effect of activation of inhibitory (GABA/Gly) neurons in the CnF, PPN, or LPGi, the locomotor-like activity was initially evoked by bilateral application of excitatory drive to glutamatergic neurons in the CnF (D_CnF,ex $\in$ [2.78, 3.06] in 0.02 steps). Then, to simulate activation of inhibitory (Gly/GABA) populations in the CnF, PPN, or LPGi, for each value of D_CnF,ex, the excitatory drive was unilaterally applied to the corresponding inhibitory population (D_{CnF/PPN/LPGi,in} $\in$ [1.15, 2.85] in 0.07 steps).

To test the robustness of the model, we simultaneously varied all connection weights by multiplying each weight by a normally distributed random number with a mean of 1 and standard deviation σ_p between 0.02 and 0.2 in steps of 0.02. For each σ_p, 100 random models were built and bifurcation diagrams were calculated. With σ_p ≤ 0.04 all randomized models retained all stable regimes and their sequential transitions with changes of α. With increasing σ_p an increasing number of models lost some stable solutions (gaits such as bound or trot) and 50% of the models were unstable at σ_p = 0.2. Thus, the final model represents a coarse system allowing parameter variations without dramatic (qualitative) changes in behavior.

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

1. Received: November 15, 2018
2. Accepted: January 19, 2019
3. Accepted Manuscript published: January 21, 2019 (version 1)
4. Version of Record published: January 31, 2019 (version 2)

© 2019, Ausborn 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.