Introduction

Although incomplete spinal cord injuries partially disrupt communication across the lesion they are often followed by substantial recovery of locomotor function (Burns and Ditunno, 2001; Kuerzi et al., 2010). These improvements are accompanied by reorganization of spared and severed descending pathways and spinal circuitry (Raineteau and Schwab, 2001; Fouad and Tse, 2008; Takeoka et al., 2014; Filli and Schwab, 2015; Engmann et al., 2020; Lemieux et al., 2024). However, it remains unclear which plastic changes of neural connectivity underlie functional recovery and if the symmetry of the lesion influences those changes.

Locomotion is primarily controlled by spinal circuitry (Grillner, 1981; Orlovsky et al., 1999; Kiehn, 2006, 2016; McCrea and Rybak, 2008; Rybak et al., 2015), which produces the locomotor rhythm and pattern, and controls interlimb coordination in response to supraspinal inputs. Flexor-extensor alternation in each limb is regulated by rhythm generators located bilaterally within the lumbar and cervical enlargements (Kato, 1990; Ballion et al., 2001; Juvin et al., 2012; Rybak et al., 2015; Danner et al., 2016, 2017; Frigon, 2017; Fig. 1A). These rhythm generators receive descending inputs from the brainstem, which modulate locomotor frequency and speed (Danner et al., 2016; Capelli et al., 2017; Caggiano et al., 2018; Josset et al., 2018; Ausborn et al., 2019), and are interconnected by commissural (left–right) and long propriospinal (lumbar–cervical) neuronal pathways, which define interlimb coordination (Talpalar et al., 2013; Danner et al., 2016, 2017; Ruder et al., 2016). Thoracic spinal cord injuries disrupt these descending pathways and inter-enlargement long propriospinal neuronal connections, while leaving the sublesional lumbar and supralesional cervical circuitry largely intact (Fig. 1B,C). Thus, the animals’ ability for speed-dependent gait expression following post-injury recovery can provide insights into the state of reorganization of descending and inter-enlargement pathways (Shepard et al., 2021, 2023; Danner et al., 2023).

Figure 1.

Figure 2.

Indeed, following thoracic lateral hemisection or moderate midline contusion injuries, rats recover weight-bearing overground locomotion, though with a reduced maximal speed and altered speed-dependent gait patterns (Danner et al., 2023). Specifically, post-hemisection recovery allows expression of a subset of pre-injury gaits — walk, trot, canter, and gallop — in a speed-dependent manner, albeit with reliance on the contralesional limb as lead in non-alternating gaits. In contrast, mild-to-moderate midline contusion injuries result in a more pronounced speed reduction, loss of non-alternating gaits, and emergence of novel alternating gait patterns (Danner et al., 2023).

We developed a series of computational models of spinal circuitry that controls speed-dependent interlimb coordination and gait expression (Danner et al., 2016, 2017; Ausborn et al., 2019, 2021; Zhang et al., 2022). By reproducing the speed-dependent gait patterns of intact mice (Bellardita and Kiehn, 2015; Lemieux et al., 2016) and simulating the effects of silencing or ablating various classes of commissural (Talpalar et al., 2013; Bellardita and Kiehn, 2015) and long propriospinal interneurons (Ruder et al., 2016; Zhang et al., 2022), these models offer insights into the organization and function of spinal locomotor circuitry, including the roles of commissural and long propriospinal interneurons as well as the broader interactions between brainstem drive and the spinal circuitry that regulate gait expression. Given the similarity in gait patterns between rats (Danner et al., 2023) and mice (Bellardita and Kiehn, 2015; Lemieux et al., 2016), we adapted these models to rats as a foundation for studying post-injury plasticity.

Here, we use a computational model of spinal circuitry to investigate potential reorganization of neural connectivity underlying locomotor control following recovery from symmetrical thoracic contusion and asymmetrical (lateral) hemisection injuries. The model reproduces speed-dependent gait expression in rats before injury and after recovery from hemisection and contusion. In our model, reproduction of altered gait expression observed after recovery from the fully asymmetrical hemisection injury required substantial functional restoration of severed long propriospinal connections and descending drive. This suggests the involvement of compensatory mechanisms, such as local axon sprouting and rewiring, resulting in formation of detour pathways via the intact hemicord. In contrast, after a symmetrical contusion injury, the model suggests that the observed post-recovery gait changes rely on strengthened intra-enlargement connectivity, increased afferent input, and/or changes in intrinsic excitability of the sublesional circuitry that compensate for weakened long propriospinal connections and descending drive.

Results

Model of the rat spinal locomotor circuitry

The intact (pre-injury) model (Figs 1A and 2) was based on our previous models of mouse speed-dependent gait expression (Danner et al., 2017; Zhang et al., 2022) and was adapted to reproduce rat locomotion. Like the previous models, it consists of four rhythm generators (RGs), located on the left and right sides of the cervical and lumbar compartments (left and right fore and left and right hind RGs). Each RG includes a flexor and extensor center (F and E) that mutually inhibit each other via inhibitory interneuron populations (InF and InE). Both flexor and extensor-centers incorporate a persistent (slowly-inactivating) sodium current (I_NaP, see Methods), allowing them to intrinsically generate rhythmic activity under certain conditions defined by the level of excitation, which is regulated by descending brainstem drive. The flexor centers operate in a bursting mode, while the extensor centers receive relatively high drive that keeps them in the mode of tonic activity and exhibit rhythmic activity only due to rhythmic inhibition from the corresponding flexor centers. Thus, each RG produced alternating flexor-extensor activity defining two major locomotor phases (flexor and extensor). The frequency of RG oscillations was defined by brainstem drive to each flexor center; it increased with increasing drive, primarily due to shortening of the extensor phase durations. This reflects the asymmetric modulation of stance and swing phase durations observed in limbed locomotion across species (Halbertsma, 1983; Hildebrand, 1989; Bertram, 2016), including rats (Cohen and Gans, 1975; Danner et al., 2023), where increases in locomotor frequency are primarily achieved by shortening stance (extensor) phases while swing (flexor) phases remain relatively stable — a pattern also generated intrinsically by spinal locomotor circuits across species (Frigon and Gossard, 2012; Danner et al., 2015; Shevtsova et al., 2015).

The left and right RGs within the cervical and lumbar compartments interact via a series of homologous commissural pathways mediated by commissural (V0_D, V0_V, and V3) and local (V2a, Ini and Ini1) interneurons, coordinating the activities of the left and right homologous RGs. The cervical and lumbar RGs interact via a series of homolateral (Sh2 and LPNi) and diagonal (dV0_D, dV0_V and aV3) long propriospinal pathways, coordinating the fore-hind RG activities. Homologous (intra-enlargement) V0_D and V0_V and descending diagonal V0_D interneurons receive the inhibitory brainstem drive, allowing them to coordinate left–right and fore–hind activities in a speed-dependent manner. The organization of spinal locomotor circuitry was based on multiple sets of experimental findings on speed-dependent gait expression of intact mice and after manipulating various classes of local, commissural, and long propriospinal interneurons (Talpalar et al., 2013; Bellardita and Kiehn, 2015; Rybak et al., 2015; Shevtsova et al., 2015; Danner et al., 2016, 2017, 2019; Ruder et al., 2016; Ausborn et al., 2021; Zhang et al., 2022). The detailed description of the proposed connectome of spinal circuitry can be found in Danner et al. (2017) and Zhang et al. (2022). The connection weights and neuron parameters in the current model were adjusted to reproduce characteristics of locomotor behavior in rats (see Methods).

We quantified frequency-dependent changes in inter-RG coordination and flexor/extensor phase durations in response to variations in brainstem drive (parameter α, see Methods) by calculating normalized phase differences between activities of pairs of RGs and compared the resulting gait expression to rat locomotion. Specifically, inter-RG coordination in the model was compared to interlimb coordination in rats, flexor and extensor phase durations were compared to swing and stance phase durations, and the period of oscillation was compared to the step-cycle duration in rats. These parameters were used to calculate gait in both the model and in rats to assess frequency-dependent gait expression (see Methods).

Model reproduces gait expression of intact rats

To ensure suitability of the model to study potential circuit reorganization following recovery from spinal cord injury, we first compared speed-dependent gait expression of the model in the intact case (intact or pre-injury model) with that of intact rats while traversing a 3-m long tank (Danner et al., 2023). Locomotor bouts consisting of periods of acceleration and deceleration were simulated by successively and repetitively ramping the brainstem drive (parameter α, see Methods) from a value at which trot is usually expressed (α=0.55) to one at which bound is expressed (α=1.05) and back down; step-to-step variability was insured by including a noisy current (σ_Noise = 1.1 pA; see Methods).

Representative examples of bouts of locomotor activities of the intact model and rat are shown in Fig. 3A1 and 3A2, respectively. Step frequencies increased progressively in both the model and rats, driven primarily by shortening of the extensor and stance phases (Fig. 3A1 and 3A2, upper panels). As frequency increased, the hind left–right normalized phase differences in the model and rats (orange lines in Fig. 3A1 and 3A2, bottom panels) gradually shifted from alternation (∼0.5 phase difference) to synchronization (∼0/1 phase difference). Meanwhile, homolateral phase differences remained close to alternation (green), and diagonal phase differences transitioned from synchronization to alternation (blue). As a result, both the model and rat transitioned from trot to gallop and then half-bound gallop with increasing frequency. As shown previously (Danner et al., 2023), speed-dependent gait changes were gradual, rather than abrupt, resulting in the expression of transitional gaits (Fig. 3A).

Figure 3.

The average gaits identified in the intact model (Fig. 3B1) — trot, canter, transverse and half-bound gallops (with left and right lead), and bound — closely matched those observed in rats (Fig. 3B2). In the model, trot was characterized by synchronized diagonal RG activity and alternating left– right activation of homologues fore and hind RGs. Canter exhibited synchronization of only one pair of diagonal RGs. Gallops showed quasi-synchronization of the hind RGs and a non-zero phase difference between the fore RGs. Bound was characterized by left–right synchronization of the fore and hind RGs and alternation of the homolateral and diagonal RGs. Note the striking similarity of average extensor/stance phases and normalized phase differences in modeling and experimental results (Fig. 3B1 and 3B2).

To compare frequency-dependent gait expression between the model and rats across bouts, normalized phase differences between pairs of RGs (Fig. 4A1) and corresponding pairs of limbs (Fig. 4A2) were plotted against locomotor frequency. In both cases, the frequency-dependent gait transitions from trot through gallop to half-bound gallop and bound (Fig. 4B1,B2) can be seen as progressive changes in the normalized phase differences (Fig. 4A1,A2). At low frequencies (∼3–4 Hz), left–right alternation of the hind RGs/hindlimbs and the fore RGs/forelimbs dominated. At medium frequencies (∼4–5.5 Hz), left–right phase differences showed progressive change from alternation to synchronization in both directions up to essential synchronization, while some steps still maintained left–right alternation. Finally, at the highest frequencies (>6.5 Hz in model and >6 Hz in rats), most steps exhibited left–right fore and hind phase differences close to 0 or 1 (synchronized activity). The normalized homolateral phase differences remained alternating, close to 0.5 across the frequency range. The normalized diagonal phase differences gradually changed from synchronized to alternating. In parallel, phase durations showed asymmetric frequency-dependent shortening, with a more pronounced reduction in extensor (stance) phase duration relative to flexor (swing) phase duration, consistent between the model (Fig. 4C1) and rats (Fig. 4C2).

Figure 4.

Trot in the intact model was expressed at oscillation frequencies from about 3 to 6 Hz (Fig. 4B1) and was the most prevalent gait in our simulations (∼45%, Fig. 5A1), which is similar to the intact rats (Figs 4B2 and 5A2). The transverse gallop, both in the model and in rats, was expressed at slightly higher frequencies (∼4.5 to 6.5 Hz) and accounted for ∼30% of all steps, making it the second most prevalent gait (Figs 4B1,B2 and 5A1,A2). Half-bound gallop was observed in the model at ∼5 to ∼7 Hz (Fig. 4B1), similar to rats (Fig. 4B2), and was observed in ∼10% of all steps (Fig. 5A1,A2). Bound was expressed at ∼5 to 7 Hz in both the model and rats and accounted for ∼10% of steps (Figs 4B1,B2 and 5A1,A2). Canter was the least prevalent gait in the intact model (∼0.5% of steps; not shown in Fig. 4B1), and also occurred rarely, mostly as a transitional gait, in rats (Danner et al., 2023). Although the exact proportions of the expressed gaits in the intact model and rats differ (Fig. 5A1,A2), their frequency-dependent distribution is qualitatively similar (compare Fig. 4B1 and 4B2). Importantly, both in the model and intact rats, the ranges of locomotor frequency where particular gaits were expressed overlapped (Fig. 4B1,B2) and transitions between gaits occurred most prevalently between gaits in the sequence of their frequency-dependent expression (Figs 4B1 and 5C1,D1), similar to the gait transitions in the intact rats (Fig. 4B2 and Fig. 2 in Danner et al., 2023).

Figure 5.

Finally, to assess the step-to-step variability while accounting for variance introduced by gait changes, we measured the mean deviations from the circular exponential moving averages of each phase difference per bout (Danner et al. 2023). Both the model and rats exhibited a similar pattern: the mean deviations of the moving averages of the two left–right and diagonal phase differences were comparable to each other and consistently higher than those of the two homolateral phase differences (Fig. 5B1,B2).

In summary, the model reproduces the experimental results on locomotion of intact rats well: both the model and rats show progressive gait transitions with frequency increase (from trot to gallop and then to half-bound gallop and to bound) on a continuum in phase-space, and gaits were expressed in overlapping frequency ranges.

Modelling recovery from lateral thoracic hemisection injury

Following recovery from the asymmetrical lateral thoracic hemisection, rats regain the ability to locomote over a wide range of speeds but lose the ability to perform the fastest gaits (half-bound gallop and bound) and exhibit clear left–right asymmetry in interlimb coordination, predominantly using the contralesional limb as the lead during asymmetric gaits such as canter and gallop (Danner et al., 2023). Nevertheless, the expressed gaits represent a subset of intact gaits with appropriate interlimb coordination. Here, we used the model to investigate potential reorganization of spinal locomotor circuitry that could underlie the observed speed-dependent gait expression following recovery from hemisection.

To do this, we first modeled the acute effect of the injury by eliminating severed connections from the model (Figs 1B and 6A). Specifically, right thoracic hemisection was modeled by eliminating brainstem drive to the ipsilesional lumbar compartment and all long propriospinal neuron (LPN) connections on the side of the injury. Note that axons of contralaterally projecting LPNs generally cross the midline at a spinal level close to the cell body (Reed et al., 2006; Pocratsky et al., 2020), thus we assumed that contralaterally projecting LPNs with contralesional cell bodies and homolaterally projecting LPNs with ipsilesional cell bodies are axotomized by the thoracic hemisection.

Then, by manually tuning model parameters, we sought to identify a small set of connectivity changes sufficient to reproduce speed-dependent gait expression observed experimentally in rats following recovery from hemisection injury (Fig. 6B). Specifically, we assumed that the injury-affected LPN and descending pathways, along with local cervical and lumbar circuitry, undergo plasticity associated with functional recovery, e.g., through regrowth (sprouting), detour pathways, synapse loss/formation, or upregulation of sensory feedback, and that descending control of lumbar and cervical circuitry through uninjured pathways is altered.

Figure 6.

These assumptions were implemented in the model (Fig. 6B) by:

1. (functional) recovery of the injury-affected LPN connections to 40% of their pre-injury weights,

2. recovery of drive to the ipsilesional hind rhythm generator to 90% of its pre-injury strength (this might be partially due to concomitant upregulation of sensory input),

3. recovery of inhibitory drive to the ipsilesional lumbar V0_V commissural interneuron (CIN) to 50% of its pre-injury strength, and

4. reduction of inhibitory drive to cervical and contralesional lumbar V0_V CINs to the same strength (50%), which through disinhibition strengthened left–right alternation.

The post-hemisection model exhibited frequency-dependent expression of interlimb coordination and gait qualitatively similar to those observed in rats following recovery from hemisection (Figs 5, 7, and 8 and Fig. 2 in Danner et al., 2023). Exemplary bouts of locomotor activities from both the model and animals showed progressive transitions from trot to canter, followed by transverse gallop and back to trot (Fig. 7A1,A1, upper panels). In both the model and in rats, hind left–right normalized phase differences during gallop and canter were above 0.5, signifying use of the left-lead variants of both gaits (Fig. 7A1,A2, bottom panels). Across all steps, average extensor phase timing and the normalized phase differences for each gait were qualitatively similar between the modeling and experimental results (Fig. 7B1,B2), with the most noticeable difference being the left–right hindlimb phase difference during trot, which in the model showed a stronger asymmetry, i.e., was shifted farther from alternation.

Figure 7.

Figure 8.

In both the model and in rats, the maximal locomotor frequencies were slightly reduced (compared to the pre-injury case shown in Fig. 4) to approximately 6-6.5 Hz (Fig. 8), and the high-frequency, synchronized gaits (half-bound gallop and bound) were expressed in very few steps (see also Figs 5A1,A2 and 8). On the other hand, the occurrence of canter notably increased compared to the intact case (Fig. 5A1,A2) and was distributed across a wide range of generated frequencies (Fig. 8). In both the model and rats, trot (62% of all steps in the model and ∼32% in rats) and gallop (∼25% in model and ∼36% in rats; Fig. 5A1,A2) remained prevalent and were also more stable (high probability of subsequent steps being of the same gait) than canter (Fig. 5C2,D2 and Fig. 2H2, I2 in Danner et al. 2023). Canter often occurred as a transitional gait (high probability that next step is of a different gait) and was frequently skipped at transitions from trot to gallop or gallop to trot (Fig. 5C2,D2 and Fig. 2H2, I2 in Danner et al. 2023).

The frequency-dependent distribution of phase differences across bouts (Fig. 8A1,B1) exhibits several features consistent with experimental observations in rats following recovery from hemisection injury (Fig. 8A2,B2). All steps were expressed on a continuum of phase differences that appears to be a subset of the phase differences of their pre-injury counterparts (Fig. 4A1,B1). Hind left–right normalized phase differences showed a progressive deviation from alternation (∼0.5 phase difference) to quasi-synchronization (phase difference close to 0/1) that occurred only in the direction of the left-lead canter and gallop; right-lead canter and gallop were almost completely lost (Fig. 7C1,C2). Hind left–right synchronization (phase difference close to 0/1; characteristic for half-bound gallop and bound) was almost never reached. Furthermore, there was a left–right asymmetry of the homolateral phase differences: the ipsilesional homolateral normalized phase differences exhibited lower values compared to pre-injury and to the contralesional homolateral phase differences.

Finally, the model exhibited changes in variability of phase differences relative to the pre-injury condition that were similar to those observed in rats (Fig. 5B1,B2). Specifically, compared to pre-injury, the variability (mean deviation from their moving average) of fore left–right and the contralesional (left) homolateral phase differences decreased, whereas variability of the ipsilesional (right) homolateral phase difference increased. The variability of the other phase differences remained at a similar level to their pre-injury values, with a slight increase in the hind left–right phase difference and a slight decrease of the diagonal ones. Although the hind left–right phase difference increased more strongly in rats (Danner et al. 2023), the model results overall are in qualitative agreement with the experiments. It is interesting to note here that the variability of the contralesional (left) homolateral phase difference decreased in the model even though the contralesional LPN connection weights remained at pre-injury levels.

Altogether, these results show that the suggested reorganization of long propriospinal connectivity and descending drive allowed the model to reproduce key features of rat gait expression following recovery from thoracic hemisection injury: expression of a subset of pre-injury gaits, loss of the highest-speed gaits, and pronounced left–right bias evident in the preferential use of contralesional RG/limb as lead during left–right asymmetric gaits.

Mechanism underlying left–right bias after hemisection

Recovered locomotion following hemisection injury was modeled as a partial (functional) recovery of the pathways severed by the injury and thus, the model included asymmetries of LPN connectivity as well as descending drive (Figs. 1B and 6). To investigate the contributions of these structural asymmetries on the emergence of the left–right bias and loss of high-speed gaits in the post-hemisection model, we created bifurcation diagrams of the normalized phase differences for the brainstem drive parameter α (Fig. 9). These simulations were performed at a low level of noise (σ_Noise = 5 fA, see Methods) to ensure identification of stable trajectories.

Figure 9.

In the intact case, the model generated stable rhythmic locomotor-like activity when parameter α was changed in a stepwise manner from ∼0.1 to 1.1 (Fig. 9A1), which led to an increase in locomotor frequency from ∼1 to ∼7 Hz (Fig. 9B1). Blue and red lines in Fig. 9A1 indicate the stable phase differences with increase or decrease of α, respectively. Discrepancies between the red and blue lines (when ∼0.8<α<0.92) indicate regions of hysteresis and multi-stability. For these α values, multiple stable solutions coexisted and could be expressed depending on the initial conditions. The increased frequency was accompanied by sequential gait changes from lateral-sequence steps to trot to gallop and then half-bound gallop and bound (Fig. 9A1). In the region of multi-stability, both trot and gallop were stable. Importantly, for the left–right asymmetric gaits (gallop and half-bound gallop), where the hind and fore left–right normalized phase differences are not at 0.5 or 0/1, both left–right symmetric solutions were stable. This corresponds to the stability of gaits with left and right lead RGs. The gait transitions and the corresponding bifurcations in the intact model are similar to our previous mouse models and are described in detail in Danner et. all (2017). The frequency range and step cycle compositions (relative phase durations) for all gaits as well as gait expression were consistent with rat locomotion (Danner et al. 2023).

In the case of simulated hemisection, the bifurcation diagrams show (Fig. 9A2) that the range of brainstem drive parameter α, in which the model is stable, was reduced compared to the intact case. This resulted in a reduction of the maximal frequency to ∼6 Hz (Fig. 9B2). With increasing drive, the model transitioned from a gait in-between lateral-sequence and canter to trot and then gallop (Fig. 9A2). Half-bound gallop and gallop were lost as stable solutions, and all phase differences lost their second branches. The hind left–right phase difference was shifted above 0.5 (perfect alternation) and progressively increased with increasing α, resulting in continuous transition from trot to gallop only in the direction of left (contralesional) lead sequences. The fore left–right, the contralesional (left) homolateral, and the diagonal phase difference involving the contralesional hind RG all maintained a similar trajectory to their pre-injury counterparts at low α-values before the trot-to-gallop bifurcation. On the other hand, the ipsilesional (right) homolateral and diagonal phase difference involving the ipsilesional hind RG were shifted compared to their pre-injury trajectories.

Next, we investigated how injury-affected LPN connectivity and drive to the ipsilesional hind RG contributed to altered gait expression following hemisection injury. We created bifurcation diagrams for simulations in which only the injury-affected LPN connections (Fig. 9A3) or the drive to the ipsilesional hind RG (Fig. 9A4) were set to their post-injury recovered values (40% and 90% respectively), while all other connection and drive weights were left at their intact values.

Interestingly, the lowered connection weights of the injury-affected and partially recovered LPN weights caused a progressively increasing bias of the left–right phase differences with increasing speed towards ipsilesional (right) lead gallops (Fig. 9A3) — the opposite direction as observed after hemisection. This is followed by a bistable bifurcation (at ∼0.85<α<0.95) towards left-lead sequence gallops.

On the other hand, the slightly (90% of pre-injury value) lowered drive to the ipsilesional hind RG (Fig. 9A4) caused a left–right asymmetry with a progressive shift away from trot towards gallop and half-bound gallop in the same direction as observed in rats after recovery from hemisection — towards left–lead (contralesional) sequence gaits. Indeed, all six phase differences exhibit trajectories similar to the simulated hemisection (Fig. 9A2), although with stronger deviations of the left–right fore and hind phase differences from alternation at higher frequencies.

Thus, in the model, injury-induced and partially recovered left–right asymmetries of the LPN connections and descending drive to ipsilesional RG have opposing effects on left–right asymmetry of interlimb coordination and gait expression. The hemisection model is a result of the interaction of both of these influences, but the effect of the lower drive to ipsilesional hind RG compared to the other three RGs dominates and largely causes the observed left–right bias.

Modelling recovery from thoracic contusion injury

A midline (left–right) symmetrical thoracic contusion injury (Figs 1C and 10A) partially damages the ascending and descending pathways on both sides of the spinal cord, including supraspinal inputs and long propriospinal connections. However, similar to the hemisection injury, the cervical and lumbar circuitry including their homologues pathways remain intact. After a moderate thoracic midline contusion injury, rats regain functional locomotion dominated by left–right alternating gaits (Danner et al. 2023), and even slower maximum speed than after hemisection. They lose the ability to perform non-alternating gaits (gallops and bounds) and adopt atypical gaits not observed in uninjured animals, such as diagonal-sequence and pace. These two novel gaits are characterized by left–right alternation of both the forelimbs and hindlimbs. The diagonal sequence is a four-beat gait pattern where a hindlimb step is followed by a step of the contralateral (diagonal) forelimb, rather than the ipsilateral forelimb, as is the case in the lateral sequence. Pace is a two-beat gait in which the homolateral limbs move in synchrony. Previously, we hypothesized that the post-contusion gait changes are a result of weakened fore–hind coupling together with maintenance of left–right alternation (Danner et al. 2023).

Figure 10.

We assumed that after symmetrical thoracic contusion injury the brainstem drive and ascending and descending LPN connections on both sides of the model are significantly damaged (Figs 1C and 10A). Thus, to model circuit reorganization underlying recovered post-contusion locomotion, we assumed that spared sublesional lumbar circuitry undergoes plastic reorganization and that supraspinal control over the intact supralesional cervical circuitry is altered. Then, by manually tuning model parameters, we looked for possible changes of connection weights in the model that result in reproduction of locomotor function of rats following recovery from contusion, while also keeping the connection weights of the injured pathways low. Note that we did not make any assumption about the proportion of the translesional fibers that remain functionally intact.

This process resulted in the following specific suggestions that we implemented in the model:

1. reduced brainstem drive to the sublesional lumbar RGs was substituted with enhanced afferent input and/or changes in their intrinsic excitability (simulated as additional drives to the hind RGs);

2. descending control over the supralesional cervical circuitry via (tonic) brainstem drive was adjusted so that the oscillation frequency of the cervical RGs matched those of the lumbar RGs (simulated by equal drives to cervical and lumbar RGs) and to secure fore left–right alternation (reduction of inhibitory drives to cervical V0 CINs);

3. commissural pathways in the lumbar compartment were reorganized to secure hind left–right alternation (descending drive to V0_V CINs was removed, and their excitability was altered by an additional constant drive) (Fig. 10B); and

4. finally, connection weights of all inter-enlargement LPN connections were set to 5% of their pre-injury values.

An exemplary bout produced by the post-contusion model is shown in Fig. 11A1. The model progressively transitioned from a diagonal sequence gait to trot then across lateral sequence to pace and finally back to the diagonal-sequence gait. The example bout of rat post-contusion recovery (Fig. 11A2) exhibited the same sequence of gait transitions except that it starts and ends with a lateral sequence gait pattern. Note that fore and hind left–right alternation was maintained across all steps of the bouts, while the homolateral and diagonal phase differences progressively drift in parallel across the full cycle reaching their approximate initial values at the end of the bouts (Fig. 11A1,A2, bottom panels). In both cases, the fore limbs/RGs performed one more step cycle than the hind limbs/RGs (Fig. 11A1,A2, upper panels).

Figure 11.

Across simulated bouts, the post-contusion model expressed lateral and diagonal-sequence gait patterns, as well as trot and pace (Fig. 5A1). Like contused rats (Fig. 5A2; Danner et al. 2023), the model did not express non-alternating gaits (gallop, half-bound gallop and bound) but produced left– right alternating gaits that the intact model did not (including pace and the diagonal-sequence pattern; Fig. 11B1). The maximal locomotor frequencies were reduced up to approximately ∼6 Hz (Fig. 12) and average gaits of the model exhibited similar extension phase patterns and phase differences to those observed in rats (Fig. 11B1,B2). Across the expressed range of frequencies, both fore and hind left– right normalized phase differences remained alternating (around 0.5) — which is in contrast to the intact case where fore and hind left–right phase differences progressively transition from alternating to in-phase synchronization with increasing frequency (Fig. 12A1,A2). On the other hand, all fore–hind phase differences (both homolateral and diagonal phase differences) were highly variable and distributed across the whole cycle. Frequency-dependence of gait expression was largely lost, with all gaits expressed across the full range of frequencies (Fig. 12A1,A2). Yet, in the model diagonal-sequence steps were more frequent at low frequencies and lateral-sequence steps at high frequencies. This is different to contused rats, which expressed these gaits on average at comparable frequencies (Fig. 12B1,B2).

Figure 12.

Interestingly, gait transition probabilities in the post-contusion model showed a clear pattern (Fig. 5C3,D3) and apart from the not-defined category (capturing all other gaits), the expressed gaits were relatively stable: a step of one gait was likely to be followed by another of the same gait (all probabilities 0.74–0.85). Transitions predominantly occurred between neighboring gaits along the sequence: lateral-sequence, trot, diagonal-sequence, pace, and back to lateral-sequence (Fig. 5C3,D3).

The idealized versions of these gaits are evenly distributed along a line (or rather circle) in phase space (defined by the three orthogonal phase-differences: hind left–right, homolateral, and diagonal). The line connects trot and pace in both directions (through lateral-sequence in one and diagonal-sequence in the other direction) and represents the only remaining degree of freedom, assuming that fore and hind left–right alternation was maintained. This organization also explains the experimentally observed distribution of phase differences in contused rats (Danner et al. 2023).

These results are further complemented by the analysis of phase difference variability (Fig. 5B1,B2). Specifically, the mean variability (mean deviation from the moving average) of both fore and hind left–right phase differences were lower compared to the variability of all fore–hind phase differences in both the model and in rats. Comparing variability post-contusion to the pre-injury condition revealed some discrepancies between the model and experimental results. In the model, variability of fore and hind left–right phase differences decreased after contusion, whereas the experimental results showed no significant change. This suggests that the model may have overestimated the strengthening of left–right coupling due to disinhibition of cervical V0_V CINs and reorganization of the lumbar circuitry. However, the mean variability of all fore–hind phase differences (left and right homolateral and diagonal) increased in both the model and experimental results.

Overall, the model supports our previous hypothesis that, following a symmetrical moderate thoracic contusion injury coupling between the left and right RGs within each enlargement remained at (or recovered back to) pre-injury levels, while the coupling between the lumbar and cervical circuits was significantly reduced. These results suggest that plasticity with the local circuitry linking the left and right RGs within each enlargement are critical for recovery following a symmetrical thoracic contusion injury.

Discussion

The brain and spinal cord have a high capacity to adapt and reorganize. Indeed, following incomplete spinal cord injuries, neural plasticity occurs across the breadth of the nervous system (Filli and Schwab, 2015). This includes reorganization of supraspinal motor centers, compensatory structural and functional plasticity of spared descending, propriospinal, and sensory neurons above and below the lesion (Ballermann and Fouad, 2006; Rosenzweig et al., 2010; Filli et al., 2014; Takeoka et al., 2014; Zörner et al., 2014), local regrowth/sprouting of transected fibers, as well as intrinsic changes to sublesional spinal interneurons and motoneurons (Barriere et al., 2008; Cohen-Adad et al., 2014; Takeoka and Arber, 2019). In this computational modeling study, we have investigated potential reorganization strategies of neural connectivity underlying recovery of locomotor function following thoracic lateral hemisection and moderate midline contusion injuries in rats, representing asymmetrical and symmetrical injuries, respectively, that exhibit roughly comparable functional recovery (Danner et al., 2023). The model itself represents neuronal populations resident in the cervical and lumbar enlargements with axons of ascending and descending inter-enlargement propriospinal neurons along with reticulospinal drive to the individual hind and fore RGs. The model reproduced key features of speed-dependent gait expression of rats before and after spinal cord injury. It suggests that recovery after lateral hemisection involves partial functional restoration of descending drive and long propriospinal pathways, whereas recovery following midline contusion relies on reorganization of sublesional lumbar circuitry combined with altered descending control of cervical networks.

Descending drive from the reticular formation controls locomotor speed and gait in intact animals (Capelli et al., 2017; Ausborn et al., 2019; Hsu et al., 2023) and substantial structural and functional plasticity of reticulospinal neurons occurs in parallel with recovery of locomotor function following hemisection injury (Ballermann and Fouad, 2006; Filli et al., 2014; Engmann et al., 2020; Lemieux et al., 2024). In our model, drive from the reticular formation controls the oscillation frequency of each rhythm generator and modulates the activity of commissural interneurons controlling interlimb coordination. The model suggests that following the asymmetric lateral hemisection injury, the drive to the ipsilesional hind rhythm generator recovers to around 90% of its pre-injury strength (Fig. 6B).

The strength of this drive reflects the net excitatory input to the population of flexor rhythm generator neurons and cannot be easily interpreted in terms of the number of regenerated axons or synaptic terminals. Thus, while these results suggest substantial reorganization of reticulospinal interactions with the lumbar circuitry, other factors, such as changes in intrinsic neuronal properties, increased prevalence of excitatory connections between constituent neurons, and strengthened afferent input (Barriere et al., 2008; Cohen-Adad et al., 2014; Takeoka et al., 2014; Takeoka and Arber, 2019; indicated by the pink arrow to the ipsilesional hind RG in Fig. 6B) can also affect the excitability of the sublesional rhythm generator neurons and help compensate for the reduced descending drive.

Inter-enlargement long propriospinal neurons couple the cervical circuits controlling the forelimbs with the lumbar circuits controlling the hindlimbs and are essential for interlimb coordination (Figs 1A and 2; Edinger, 1896; Nathan and Smith, 1959; Dutton et al., 2006; Ruder et al., 2016; Flynn et al., 2017; Frigon, 2017). These neuronal populations are roughly 50% ipsilaterally and 50% contralaterally projecting with the axons of the latter crossing the midline at the level of the cell bodies (Reed et al., 2006; Pocratsky et al., 2020; Shepard et al., 2023). Lateral thoracic hemisection injury severs all these axons unilaterally (Fig. 6A). Similar to the predicted recovery in descending drive, the model also predicts that partial (functional) recovery of the severed long propriospinal pathways is helpful and possibly necessary for recovery (Fig. 6B). This recovery was simulated by setting the connection weights of the injured long propriospinal neuron populations and drive connections to a proportion of their pre-injury strengths. Long-distance regeneration of central axons is unlikely to occur in adult animals (Filli and Schwab, 2015), thus reorganization of these pathways likely depends on local axon growth and synaptic rewiring onto propriospinal neurons with intact axons, resulting in formation of detour pathways through the intact hemicord (Courtine et al., 2008), as well as compensatory mechanisms of spared reticulospinal and long propriospinal neurons, including increased innervation of the ipsilesional lumbar hemicord via sprouting (Ballermann and Fouad, 2006; Courtine et al., 2008; Takeoka et al., 2014; Filli and Schwab, 2015). Indeed, recovery of locomotor function after hemisection injury has been shown to result in an increased prevalence of long propriospinal neurons bypassing the injury by crossing the midline twice and of reticulospinal connections to the ipsilesional lumbar cord (Takeoka et al., 2014). The work by Reed (Reed et al., 2006) uncovered evidence for a small number of descending inter-enlargement neurons with axons that crossed the midline close to the cell body and then again in the lumbar enlargement to terminate ipsilaterally in uninjured animals suggesting that a descending detour pathway may already exist, even if small. Furthermore, potentiation of ipsilesional muscle responses to optogenetic stimulation of the reticular formation has also been linked to locomotor recovery (Lemieux et al., 2024).

Recovery from lateral thoracic hemisection results in distinct left–right asymmetries of locomotor behavior (Danner et al., 2023). These asymmetries are most obvious at higher speeds and for non-alternating gaits — animals almost exclusively use the contralesional limb as the lead limb. We investigated the individual contribution of the asymmetries of long propriospinal connections and of the descending drive of the simulated hemisection (Fig. 9A3,A4). Interestingly, in the case where only the asymmetrical drive to the lumbar rhythm generators (90% of pre-injury strength to ipsilesional RG) was included and all other connections were set to their pre-injury values, the model (Fig. 9A4) very closely reproduced frequency-dependent changes of interlimb coordination and gait expression of hemisected animals. On the other hand, the asymmetry introduced by the partially recovered long propriospinal connectivity resulted in gait asymmetries in the opposite direction (ipsilesional limb used as lead limb for non-alternating gaits; Fig. 9A3). Thus, the asymmetry of the drive to the ipsilesional hindlimb rhythm generator compensated for the asymmetry of the long propriospinal connections and was the main contributor to the observed gait asymmetries (Figs 7,8 and 9A2).

In the case of the midline, symmetrical contusion injury, our model suggests that connection weights of 5% of their pre-injury values for the long propriospinal neurons are sufficient to reproduce the experimentally observed changes in interlimb coordination and gait expression. These values are much lower than anticipated based on the proportion of axons (white matter) spared by this type of contusion injury, which affects approximately 70–80% of the overall white matter (Basso et al., 1996) but presumably a much smaller proportion of the axons in the outermost layers of the ventrolateral funiculus where the long propriospinal tracts are located (Basso et al., 1996; Reed et al., 2006; Brown et al., 2024). Note that we assumed that post-contusion recovery, descending drive to forelimb rhythm generators is tightly controlled to match level of excitability of the lumber rhythm generators. This assumption was modeled by adjusting the drive so that oscillation frequency of cervical and lumbar rhythm generators match. In reality, the excitability of the lumbar rhythm generators would be affected by somatosensory afferent input and body biomechanics (Rossignol et al., 2006; Nishikawa et al., 2007; Akay et al., 2014; Takeoka et al., 2014; Takeoka and Arber, 2019; Ausborn et al., 2021; Kim et al., 2022). We introduced noise to all neurons to simulate step-to-step variability, but we did not consider afferent feedback or limb/body biomechanics in the model. Thus, we cannot exclude that the long propriospinal connection strengths necessary to reproduce experimental interlimb coordination were underestimated.

Yet, there is experimental support for plasticity resulting in decreased long propriospinal connectivity. Specifically, using viral-based tracing, we recently found that the number of axons/terminals in the cervical enlargement arising from ascending long propriospinal neuron was reduced after recovery from a mild contusion injury (Brown et al., 2024), and further that silencing of a proportion of either ascending or descending long propriospinal neurons after contusion injury actually improves (normalizes) interlimb coordination and other key locomotor indices (Shepard et al., 2021, 2023). Thus, long propriospinal neurons seem to have a disruptive effect on interlimb coordination after injury. One explanation might be that the injury affects different classes of long propriospinal neurons (ascending vs. descending, homolateral vs. diagonal, left vs. right) differently, leading some to make new connections that are maladaptive to recovery and that (further) reduction of their connectivity via silencing reduces their influence and thereby improves interlimb coordination. In support of this idea, in Brown et al. (2024), we also found that ascending propriospinal innervation of segments just caudal to the thoracic injury was increased in the same animals where input to C6 was decreased, all relative to uninjured animals. Thus, unguided plasticity of long propriospinal neurons post-injury seems to have a disruptive effect on interlimb coordination.

Silencing of ascending or descending long propriospinal neurons in intact mice or rats has the opposite effect compared to following contusion; it disrupts interlimb coordination in a speed and context dependent manner, causing transient periods of stepping with disordered left–right coordination (Ruder et al., 2016; Pocratsky et al., 2020). Our previous models attributed this effect to the excitatory diagonally projecting long propriospinal neurons (ascending V3 LPNs, descending V0_V LPNs), which ensure diagonal synchronization for trot but also promote left–right alternation because of the tightly controlled homolateral fore–hind alternation by ipsilaterally projecting long propriospinal neurons (Danner et al., 2017; Ausborn et al., 2021; Zhang et al., 2022). Thus, the influence of long propriospinal interneurons promoting left–right alternation of cervical and lumbar rhythm generators is reduced following contusion injury. Indeed, somewhat more severe contusion injuries that remove more than 80% of white matter result in an increased prevalence of disrupted left–right alternation, and other persistent problems such as dorsal stepping and the loss of the 1:1 forelimb/hindlimb relationship during stepping (Shepard et al., 2021, 2023). Yet, recovery from moderate contusion injuries — where plantar stepping and a 1:1 forelimb-to-hindlimb ratio are recovered, as modeled here (Danner et al., 2023) — results in well-controlled left–right alternation of both fore and hindlimbs. Even animals that show the poorest recovery regain stepping, and the number of steps with appropriate left–right alternation significantly exceed the number of disrupted steps (Shepard et al., 2021, 2023). Together, these findings suggest that reorganization and/or altered modulation of intra-enlargement commissural interneurons occurs during recovery to strengthen left–right alternating influences. In the model this was implemented by reduced descending inhibition of cervical V0 commissural interneurons and altered excitability of lumbar V0 commissural interneurons. This specific implementation reflects one possible scenario, as assumed by the model, where V0 commissural interneuron excitability is regulated by descending inhibition. Yet, descending excitation (Danner et al., 2016) or descending neuromodulatory influences — such as serotonergic, noradrenergic, or dopaminergic input — may also contribute (Butt and Kiehn, 2003; Jordan et al., 2008; Abbinanti and Harris-Warrick, 2012; Abbinanti et al., 2012; Perrier and Cotel, 2015). In addition, local excitatory V2a (Crone et al., 2009) and local inhibitory interneurons contribute to left–right coordination (Kjaerulff and Kiehn, 1997; Butt and Kiehn, 2003; Shevtsova et al., 2015; Danner et al., 2016, 2017; Kiehn, 2016) and may be modulated by descending input, making them potential targets for reorganization following injury.

In conclusion, challenging animals to perform overground locomotion across the full range of speeds before and after two distinct spinal cord injury models revealed injury-specific deficits in interlimb coordination and gait expression (Danner et al., 2023). These deficits differed between the two injury models and extended beyond those observed at walk–trot speeds only. By reproducing these behavioral changes in speed-dependent gait expression in a computational model of spinal locomotor circuitry, we identified potential principles of circuit reorganization underlying post-injury locomotor recovery following asymmetrical lateral hemisection and symmetrical contusion injuries at the low-thoracic level.

The model qualitatively recapitulated the in vivo locomotor characteristics observed after recovery from lateral hemisection when the ipsilesional descending drive and severed long propriospinal neurons were partially recovered, suggesting that the plastic changes responsible for recovery involved new and/or strengthened pathways that bypass the lesion to drive the hindlimb rhythm generators. In contrast, the model changes that effectively recapitulated locomotor recovery after a symmetrical midline contusion injury involved strengthening of the local left–right circuitry within each enlargement.

These findings suggest that injury symmetry partly dictates the location and type of plasticity supporting recovery. Accordingly, clinical approaches — including neuromodulation and strategies to enhance regeneration or sprouting — may be more effective, if injury symmetry is considered. Finally, we anticipate that further development of such computational approaches will eventually support the creation of individualized models to better characterize injury profiles and assist in development of therapeutic strategies.

Methods

Experimental data

In this modeling study, we used previously published experimental data on gait characteristics of overground locomotion within the full range of speeds in intact rats and following their recovery from either a mild-moderate contusion or a lateral hemisection spinal cord injury (Danner et al., 2023). In those studies, all experiments were performed in accordance with the Public Health Service Policy on Humane Care and Use of Laboratory Animals, and with the approval of the Institutional Animal Care and Use Committee (IACUC, protocol number 19644) at the University of Louisville. No new animals were used in the current study. Locomotion was assessed in a 3-m long Plexiglas tank where the walking surface was coated with sylgard to provide increased grip. A group of 17 Adult female Sprague-Dawley rats (225–250 g) was trained, using food treats, to traverse the length of the tank. Recordings were made using four high-speed (200 Hz) video cameras placed to capture the ventral view of the animal and analyzed using custom software. Then rats were divided into two groups and given either lateral hemisection (nine rats) or midline contusion injuries (eight rats; 12.5 g cm) at the T10 level of the spinal cord. The animals were allowed to recover for at least 4 weeks and were re-introduced to the long-tank set-up and trained again. At this time their locomotor recovery plateaued (contusion: 15.7±1.9, hemisection 18.3±0.8 BBB Open Field Locomotor Scale). The post-recovery recordings were analyzed for each type of injury and compared to those for the intact animals. For more details, see Danner et al. (2023).

Modeling methods

Our model represents a network of interacting neural populations. Each population was simulated as an activity-based (firing rate) neuronal population model. In this description, the voltage variable V represents an average membrane potential of population and the output function f(V) transvers V to the integrated population activity or normalized average firing rate at the corresponding voltage V (Ermentrout, 1994). This description was used in our previous modeling studies (Rubin et al., 2009; Danner et al., 2016, 2017; Ausborn et al., 2019; Latash et al., 2020; Rybak et al., 2024) and allows for inclusion of explicitly represented ionic currents, particularly, the leakage current I_L and the persistent (slowly inactivating) sodium current, I_NaP (Rubin et al., 2009). The latter current was proposed to be responsible for generating rhythmic activity in the spinal cord (Rybak et al., 2006; McCrea and Rybak, 2007; Tazerart et al., 2007, 2008; Brocard et al., 2013) and was included in description of flexor and extensor RG centers.

The variable V for flexor and extensor centers obeyed by the differential equation:

All other populations did not include I_NaP, and in those populations, V was described as:

where C is the membrane capacitance, I_NaP is the persistent sodium current, I_L is the leakage current, I_SynE and I_SynI excitatory and inhibitory synaptic currents, respectively, and I_Noise a noisy current.

The leakage current was described as:

where g_L is the leakage conductance and E_L represents the leakage reversal potential.

The persistent sodium current in the flexor and extensor half-centers was described as:

where is the I_NaP maximal conductance, E_Na is the sodium reversal potential, and m(V) and h(V) are voltage-dependent activation and (slow) inactivation variables of I_NaP. m(V) was instantaneous and its steady state was described as

and the slow I_NaP inactivation, h(V), was modeled by the differential equation

where h_∞(V) is the inactivation steady state and τ_h(V) the inactivation time constant.

In Equations (5–8), V_1/2 and k represent half-voltage and slope of the corresponding variables (m, h, and τ); τ₀ and τ_max are the baseline and maximum of inactivation time constant τ_h, respectively.

In Equation (1) and (2), excitatory and inhibitory synaptic currents (I_SynE and I_SynI, respectively) for population i were described as:

where g_SynE and g_SynI are synaptic conductances and E_SynE and E_SynI are the reversal potentials of the excitatory and inhibitory synapses, respectively; w_ji is the synaptic weight from population j to population i (w_ji > 0 for excitatory connections and w_ji < 0 for inhibitory connections); function S is described as follows:

The output function f(V) (0 ≤ f(V) ≤ 1) was defined as

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

where d_{E,I},i is the slope and b_{E,I},i is the intercept; the free parameter α represent the strength of the brainstem drive.

The noisy current I_Noise,i for each population i was described as an Ornstein–Uhlenbeck process

where τ_Noise is the time constant, σ_Noise the standard deviation of the noise, and W_i a Wiener process.

It was numerically simulated using the finite difference formula

where is ξ_i∼𝒩(0,1) a normally-distributed random number with a variance of 1, and Δt the time step (1 ms).

Model parameters

The general neuron parameters were adapted from our previous models (Danner et al., 2016, 2017) to provide the frequency range of locomotor activity characteristic for rats. The following values of parameters were used: C = 10 pF; g_L = 4.5 nS for RG centers and g_L = 2.8 nS for all other neurons; ; g_synE = g_synI = 10 nS E_L = −62.5 mV for RG centers, E_L = −60 mV for all other populations; E_Na = 50.0 mV ; E_EsynE = −10 mV ; E_EsynI = −75 mV ; V_thr = −50 mV ; V_max = 0 mV ; V_1/2,m = −40.0 mV; k_m = −6 mV; V_1/2,h = −45.0 mV; k_h = 4 mV; τ_max = 400 ms; τ₀ = 150 ms; V_1/2,τ = −35 mV.

The following drive parameters were used for the intact model: d_E = 0.0 and b_E = 0.1 for the extensor centers; d_E = 0.1 and b_E = 0.0 for the flexor centers; d_I = 0.75 and b_I = 0.0 for homologous (within a girdle) V0_D CINs in both cervical and lumbar compartments; d_I = 1.5 and b_I = 0.0 for diagonal V0_D LPNs (dV0_D); d_I = 0.25 and b_I = 0.0 for the homologous cervical (fore) V0_V CINs; d_I = 0.15 and b_I = 0.0 for the homologous lumbar (hind) V0_V CINs. Connection weights for the intact model are listed in Table 1.

Table 1

To simulate rat behavior after recovery from right hemisection, drive to the ipsilesional flexor center was reduced by 10%, inhibitory drives to both cervical and lumbar V0_V CINs were reduced by 50%, and connection weights of LPN connections affected by the injury were reduced to 40% of their pre-injury values (see Table 2 for a comprehensive list of differences to the pre-injury case).

Table 2

To simulate rat behavior after recovery from contusion, all weights of LPN connection were reduced by 95%. In addition, the increasing inhibitory drives to lumbar V0_D and V0_V CINs were eliminated and replaced by constant drives (see Table 3 for a comprehensive list of differences to the pre-injury case).

Table 3

Computer simulations

The set of differential equations were solved using the odeint (Ahnert et al., 2011) implementation of the Cash-Karp 5_th-order Runge Kutta method (Cash and Karp, 1990) of the boost C++ library (version 1.88.0). The custom C++ code was compiled (with clang-1700.0.13.3 for MAC OS 15.3.1) as a Python package. Python 3.11.5 (Python Software Foundation, Wilmington, Delaware, USA) was used to run simulations and further analyze the results. The neural network simulation package is available at https://github.com/SimonDanner/CPGNetworkSimulator.

Data analysis

Activities of the flexor and extensor centers were used to assess the model behavior. Each rhythm generator was considered as being in flexion when output f(V) of its flexor center was greater or equal to 0.1 otherwise it was considered as being in extension. The midextension time was defined as the halfway between onset and offset of extension. The locomotor period (step) was defined as the duration between two consecutive midextension time points; the frequency as the reciprocal of the period; extension and flexion phase durations as the duration between on- and offset of extension and flexion, respectively. Normalized phase differences were calculated as the delay between the midextension time points of a pair of rhythm generators divided by the locomotor period. Six normalized phase differences were calculated: hind left–right (right hind–left hind), fore left–right (right fore–left fore), left homolateral (left fore–left hind), right homolateral (right fore–right hind), and two diagonal (right fore–left hind and left fore–right hind).

Gait classification

Gaits for each step cycle were operationally defined based on phase differences, similarly to the method used in Danner et al. (2023). Each step was represented as a point Φ in the phase space consisting of the three (non-normalized) phase differences with the left hind RG as a reference. Then the distance to all idealized one-, two-, three-, and four-beat gaits (Ψ_g; Table 4), was calculated

Table 4

and the step cycle was assigned to belong to the gait Ψ_g to which it was closest. The duty factor was not considered for gait classification. Gait transitions were analyzed on a step-by-step basis (Danner et al., 2023).

Analysis of model performance

To model step-to-step variability and gait transitions comparable with locomotor bouts of rats in the long-tank (Danner et al., 2023), simulations with increased noisy currents (σ_Noise = 1.1 pA) were performed in all three conditions: intact, post-hemisection, and post-contusion. In these simulations, α was gradually ramped up from 0.55 to 1.05 for the intact condition, or to 1.0 for hemisection and contusion conditions, and then ramped back down to 0.55. Each phase of the ramp lasted 8 seconds, and the entire ramp cycle was repeated 200 times. Variability of phase differences was quantified as the deviation of each phase difference from its circular moving average (cutoff frequency 0.125 step cycles) averaged across ramp cycles (Danner et al., 2023).

To identify stable solutions of the model, bifurcation diagrams (Danner et al., 2017) were built for the six normalized phase differences. For this purpose, α was incrementally increased from 0.075 to 1.1, and then decreased back to 0.075, using a step size of 0.001. At each step, the simulation was run for 10 s, then if the normalized circular standard deviation of the five last step cycles of each phase difference was smaller than 0.005, the final state was used as the initial condition for the next step, otherwise, the same procedure was repeated. If the normalized standard deviation was not met by 100 s of simulation, the model was considered not stable at this specific α-value. Otherwise, the mean phase differences were regarded as the stable solutions. Initial conditions were randomized, and multiple runs were evaluated. To ensure that the model does not remain on an unstable trajectory, these simulations were performed with a weak noisy current (σ_Noise = 5 fA).

Data Availability

Model configuration files, Python code for running simulations, analyzing results, and processing experimental data, as well as scripts to reproduce all simulations and generate all figures presented in the paper, are available at https://github.com/dannerlab/rat-sci-locomotion-model. The neural network simulation package used in this study is available at https://github.com/SimonDanner/CPGNetworkSimulator. The experimental data reproduced from Danner et al. (2023) can be found at https://github.com/dannerlab/rat-sci-locomotion.

Acknowledgements

This work was supported by the National Institutes of Health (NIH) grants R01 NS112304, R01 NS115900, R01 NS130799, and T32 NS121768, and the Kentucky Spinal Cord and Head Injury Research Trust.

Additional information

Funding

National Institutes of Health (R01 NS112304)

National Institutes of Health (R01 NS115900)

National Institutes of Health (R01 NS130799)

National Institutes of Health (T32 NS121768)

Kentucky Spinal Cord and Head Injury Research Trust