Introduction

Locomotion is generated and controlled by three main neural components that interact dynamically 1-6. The spinal network, including circuits of the central pattern generator (CPG), generates the basic locomotor pattern, characterized by alternation of flexor and extensor activity in each limb and coordination of activities related to left and right limbs. Supraspinal structures initiate and terminate locomotion and control voluntary aspects of locomotion. Somatosensory feedback from primary afferents originating in muscle, joint, and skin mechanoreceptors provides information on the state of the musculoskeletal system and the external environment. Although we know that the spinal network generates the basic locomotor pattern, its functional architecture, operating regimes, and the way it interacts with supraspinal signals and somatosensory feedback to produce different locomotor behaviors, including at different speeds, remain poorly understood.

In terrestrial locomotion, the step cycle of each limb consists of two main phases, a swing and a stance, that correspond to the activity of limb flexors and extensors, respectively. A main role of corresponding spinal rhythm-generating circuits is to establish the frequency of oscillations, or cycle duration, and the durations of flexor and extensor phases 7. In mammals, including humans, increasing walking speed leads to a decrease in cycle duration mostly because of a decrease in stance/extensor phase duration while swing/flexor phase duration remains relatively unchanged (reviewed in 8-10). This is observed on a treadmill in intact animals and also following complete spinal thoracic transection (spinal animals). In intact animals, supraspinal signals interact with spinal networks and sensory feedback from the limbs, whereas in spinal animals, supraspinal signals are absent.

As an experimental basis of this study, we used data previously obtained in intact and spinal cats stepping on a treadmill with two independently controlled belts in tied-belt (equal speeds of left and right belts) and split-belt (different speeds of left and right belts) conditions 11-13. Changes in cycle and phase durations during tied-belt and split-belt locomotion with increasing speed and left-right speed differences were similar in intact and spinal cats over a range of moderate speeds (0.4 – 1.0 m/s). This is somewhat surprising because intact cats need to engage supraspinal structures to maintain a fixed position on the treadmill (in relation to the external space), whereas spinal cats can only use limb sensory feedback to adjust locomotion to the belt speed. Also, an interesting observation is that intact cats cannot consistently perform treadmill locomotion at very slow speeds (below 0.3-0.4 m/s), whereas spinal cats have no problem walking at these speeds 12,14.

To investigate the organization and operation of spinal circuits and how different mechanisms interact to control locomotion, we developed a tractable computational model of these circuits operating under the control of both supraspinal drives and sensory feedback. Our model followed the commonly accepted concept that the spinal cord contains CPG circuits that can intrinsically generate locomotor-like oscillations without rhythmic external inputs 4-7,15-20. The concept of spinal mechanisms that intrinsically generate the basic locomotor pattern (CPG prototype) was initially proposed by Thomas Graham Brown 15,16 in opposition to the previously prevailing viewpoint of Charles Sherrington 21,22 that locomotion is generated through a chain of reflexes, i.e., critically depends on limb sensory feedback (reviewed in 23).

The present model is based on our previous models 24-30 and contains two rhythm generators (RGs) that control the left and right hindlimbs and interact through a series of commissural pathways. The RGs receive supraspinal drive, and limb sensory feedback during the extension phase. Each RG consists of flexor and extensor half-centers operating as conditional bursters (i.e., capable of intrinsically generating rhythmic bursting in certain conditions) and inhibiting each other. However, each RG is considered not as a simple neuronal oscillator, but as a neural structure with a dual function: as a state machine, defining operation in each state/locomotor phase independent of the phase-transition mechanisms, and as an actual oscillator, defining the mechanisms of state/phase transitions. We propose and show that, depending on conditions, each RG can operate in three different regimes: a non-oscillating state-machine regime, and in a flexor-driven and a classical half-center oscillatory regimes. In the full model, the generated locomotor behavior depends on the excitatory supraspinal drives to the RGs and on limb sensory feedback. The gain of limb sensory feedback in the intact model is suppressed by supraspinal drive via presynaptic inhibition 31-33 and is released from inhibition following spinal transection.

Our model reproduces and proposes explanations for experimental data, including the dependence of main locomotor characteristics on treadmill speed in intact and spinal cats. Particularly, based on our simulations, we suggest that locomotion in intact cats at low speeds and in spinal cats at any speed is mainly controlled by limb sensory feedback, consistent with Sherrington’s view 21, 22, whereas locomotion in intact cats at higher/moderate speeds is primarily controlled by the intrinsic oscillatory activity within the spinal network, supporting Brown’s concept 15,16.

Results

Cycle and phase durations during tied-belt and split-belt locomotion in intact and spinal cats

Tied-belt locomotion

Our previous studies in intact (Fig. 1A) and spinal (Fig. 1B) cats showed a decrease in cycle duration with increasing speed due to a shortening of the stance phase with a relatively constant swing phase duration during tied-belt locomotion 11-13. This agrees with other studies in mammals 8-10,34. When comparing intact and spinal cats (Fig. 1C), we highlight the following observations: (1) Intact cats do not step consistently during quadrupedal tied-belt locomotion if treadmill speed is at or below 0.3 m/s, whereas spinal cats have no problem performing hindlimb locomotion from 0.1 -0.3 m/s; (2) Speed-dependent changes in swing and stance phase durations from 0.4 to 1.0 m/s are qualitatively similar in intact and spinal cats, but stance duration in intact cats is usually a little longer; (3) With increasing treadmill speed, the duty factor (ratio between stance and step cycle durations) in intact and spinal cats decreases towards 0.5 (equal swing and stance proportion). In spinal cats, the duty factor reaches 0.5 at approximately 0.8 - 0.9 m/s.

Locomotion of intact and spinal cats on a tied-belt treadmill.

A and B. Step cycle, stance and swing phase durations for the right hindlimb during tied-belt treadmill locomotion of intact (A, from Frigon et al., 2015; Latash et al., 2020) and spinal (B, from Frigon et al., 2017; Latash et al., 2020) cats with an increasing treadmill speed. Data were obtained from 6-15 cycles in seven intact and six spinal cats (one cat was studied in both states). Modified from Fig. 3C, D of Latash et al., 2020, under the license CC-BY-4. C. Superimposed curves from A and B to highlight differences.

Split-belt locomotion

We used data from our previous studies in intact (Fig. 2A) and spinal (Fig. 2B) cats during split-belt locomotion, where the left (slow) hindlimb stepped at 0.4 m/s and right (fast) hindlimb stepped from 0.5 to 1.0 m/s 11-13. When comparing intact and spinal cats, we highlight the following observations: (1) Increasing the speed of the fast belt leads to a decrease in cycle duration in intact and spinal cats, but the slow and fast hindlimbs maintain equal cycle duration; (2) In the slow hindlimb, increasing the speed of the fast belt leads to a small decrease in stance and swing durations in intact cats and in swing duration of spinal cats; (3) In the fast hindlimb, increasing the speed of the fast belt produces a decrease in the stance duration and an increase in swing duration in both intact and spinal cats. The duty factor reaches 0.5 at 0.9 m/s and 0.8 m/s in intact and spinal cats, respectively, and goes below this value at 0.9 – 1.0 m/s in spinal cats, where swing duration occupies a greater proportion of the cycle.

Locomotion of intact and spinal cats on a split-belt treadmill.

A. Step cycle, stance and swing phase durations for the left (slow) and right (fast) hindlimbs during split-belt treadmill locomotion of intact cats (from Frigon et al., 2015; Latash et al., 2020). B. Changes in the same characteristics for the left (slow) and right (fast) hindlimbs during split-belt treadmill locomotion in spinal cats (from Frigon et al., 2017; Latash et al., 2020). In both series of experiments the left (slow) hindlimb was stepping at 0.4 m/s while the right (fast) hindlimb stepped with speeds from 0.5 to 1.0 m/s with 0.1 m/s increments. Data were obtained from 6–15 cycles in seven intact and six spinal cats (one cat was studied in both states). Modified from Fig. 6A,B of Latash et al., 2020, under the license CC-BY-4.

Conceptual framework and model description

Basic model architecture

In this study, we fully accept the concept that the mammalian spinal cord contains neural circuits (i.e. CPG) that can (in certain conditions) intrinsically (i.e. without rhythmic or patterned external inputs) generate the complex coordinated pattern of locomotor activity 4,6,7,15,17,18,35,36. Following our previous computational models 24-26,28-30, we assume that the circuitry of the spinal locomotor CPG includes (a) rhythm-generating (RG) circuits that control states and rhythmic activity of each limb, and (b) rhythm-coordinating circuits that mediate neuronal interactions between the RG circuits and define phase relationships between their activities (coupling, synchronization, alternation). Also, following the classical view, we assume that each RG consists of two half-centers that mutually inhibit each other and define the two major states/phases of the RG, the flexor and extensors phases, in which the corresponding sets of limb muscles are activated 6,7,15,35.

Modeling of a single rhythm generator (RG) and its operation regimes

The ability of a RG to generate rhythmic activity can be based on the intrinsic rhythmic bursting properties of one or both half-centers or can critically depend on the inhibitory interactions between half-centers, so that each half-center cannot intrinsically generate rhythmic bursting, as in the classical half-center concept 15,16. However, in the isolated mouse spinal cord using optogenetic stimulation, rhythmic activity can be induced independently in flexor- and extensor-related spinal circuits 37, suggesting that both flexor and extensor half-centers can, in certain conditions, generate independent rhythmic bursting activity (i.e. operate as conditional bursters). Previous mathematical models of neurons defined and analyzed the conditions allowing these neurons to generate rhythmic bursting based on nonlinear voltage-dependent properties of different ionic channels 38-42. Several models that described bursting activity in neurons of the medullary pre-Bötzinger complex for respiration or in the spinal cord for locomotion, suggested that this activity is based on a persistent (slowly-inactivating) sodium current, INaP 24,43-48. We implemented a similar INaP -dependent mechanism for conditional bursting in our RG half-center model (see Methods).

Figure 3 illustrates the behavior of a neuron model, describing a single conditional burster (Fig. 3A) with the output representing changes in spike frequency (see Methods). At low excitatory drive, the neuronal output is equal to zero (“silence”), but when the drive exceeds some threshold, the model switches to a “bursting” regime, during which the frequency of bursts increases with increasing drive (Fig. 3B). At a certain drive, the model switches to sustained “tonic” activity.

Let us now consider a simple half-center RG consisting of two conditional bursters, the flexor half-center (F), receiving an increasing Drive-F, and the extensor half-center (E), receiving a constant Drive-E that maintains it in a state of tonic activity if uncoupled. The two half-centers inhibit each other through inhibitory neurons, InF and InE (Fig. 3C). At low Drive-F there are no oscillations, and the F half-center remains silent, although it could start oscillating if it would not been strongly inhibited by the E half-center. The E half-center shows tonic activity because of the high value of Drive-E. We call this regime of RG operation a state-machine regime (Fig. 3D, left part of the graph). In this regime, the RG maintains a state of extension, until an external signal, sufficiently strong, activates the F half-center or inhibits the E half-center (green arrows in Fig. 3C) to release the F half-center from inhibition, allowing it to generate intrinsic bursts, switching the model to a flexion state. Increasing the excitatory Drive-F releases the F half-center from the E half-center’s inhibition (acting via the InE neuron), switching the RG to a bursting regime. Note that in this regime, the E half-center also exhibits bursting activity in alternation with the F half-center due to rhythmic inhibition from the flexor-half center via the InF neuron, although the E half-center itself, if uncoupled, operates in the tonic mode. We call this regime a flexor-driven regime (Fig. 3D, middle part of the graph). Similar to the intrinsic bursting regime in an isolated conditional burster (Fig. 3B), with an increase in Drive-F the bursting frequency of the RG is increasing (and the oscillation period is decreasing) mostly due to shortening of the extensor bursts with much less reduction in flexor burst duration. Further increasing the excitatory Drive-F leads to a transition of the RG to a classical half-center regime (Fig. 3D, right part of the graph), in which the half-centers cannot generate oscillations if uncoupled, and RG oscillations occur due to, and critically depend on, mutual inhibition between half-centers and an adaptive reduction of their responses. In this regime, with increasing Drive-F, the oscillation frequency (and period) remains almost unchanged, and flexor burst duration increases to compensate for decreasing extensor bursts. To summarize, increasing the excitatory drive to the flexor half-center in this simple RG model demonstrates a sequential transition from a state-machine regime to a flexor-driven regime, and then to a classical half-center regime.

Modeling a single conditional burster and a half-center rhythm-generator.

A. The behavior of a single INaP -dependent conditional burster. B. Changes in the burster’s output when the excitatory input (Drive) progressively increases from 0 to 1.2. With increasing Drive, the initial silence state (zero output) at low Drive values changes to an intrinsic bursting regime with burst frequency increasing with the Drive value (seen in two left insets), and then to a tonic activity (seen in right inset). C. Model of a simple half-center network (RG) consisting of two conditional bursters/half-centers inhibiting each other through additional inhibitory neurons, InF and InE. The flexor half-center (F) receives progressively increasing Drive-F, whereas the extensor half-center (E) receives a constant Drive-E keeping it in the regime of tonic activity if uncoupled. D. Model performance. At low Drive-F values, there are no oscillations in the system. This is a state-machine regime in which the RG maintains the state of extension, until an external (strong enough) signal arrives to activate the F half-center or to inhibit the E half-center (see green arrows) to release the F half-center from E inhibition allowing it to generate an intrinsic burst. Further increasing the Drive-F releases the F half-center from E inhibition and switches the RG to the bursting regime (see two insets in the middle). In this regime, the E half-center also exhibits bursting activity (alternating with F bursts) due to rhythmic inhibition from the F half-center. This is a flexor-driven regime. In this regime, with an increase in Drive-F, the bursting frequency of the RG is increasing (and the oscillation period is decreasing) due to shortening of the extensor bursts with much less reduction in the duration of flexor bursts (see bottom curves and two left insets). Further increasing the excitatory Drive-F leads to a transition of RG operation to a classical half-center oscillatory regime, in which none of the half-centers can generate oscillations if uncoupled, and the RG oscillations occur due to mutual inhibition between the half-centers and adaptive properties of their responses. Also in this regime, with an increase of Drive-F, the period of oscillations remains almost unchanged, and the duration of flexor bursts is increasing partly to compensate for the shortening of extensor bursts, which is opposite to the flexor-driven regime (see bottom curves and right inset).

Model of the spinal locomotor circuitry

The schematic of our model is shown in Fig. 4. The model incorporates neuronal circuits involved in the control of, and interactions between, two cat hindlimbs during locomotion on a treadmill. The spinal circuitry in the model includes two RGs (as described above) interacting via a series of commissural interneuronal (CIN) pathways mediated by different sets of genetically identified commissural (V0D, V0V, and V3) and ipsilaterally projecting (V2a) interneurons, as well as some hypothetical inhibitory interneurons (Ini). The organization of these intraspinal interactions was directly drawn from our earlier models 24-30 that were explicitly or implicitly based on the results of molecular-genetic studies of locomotion in mice or were proposed to explain and reproduce multiple aspects of the neural control of locomotion in these studies. Specifically, V0D and V2a-V0V-Ini pathways secure left-right alternation of RG oscillations during walking and trotting at slow and higher locomotor speeds, respectively 24,25,49. The V3 CINs contribute to left-right synchronization of RG oscillations during gallop and bound 24-26,28-30.

Model of spinal circuits controlling treadmill locomotion.

A. Model of the intact system (“intact model”). The model includes two bilaterally located (left and right) RGs (each is similar to that shown in Fig. 3B) coupled by (interacting via) several commissural pathways mediated by genetically identified commissural (V0D, V0V, and V3) and ipsilaterally projecting excitatory (V2a) and inhibitory neurons (see text for details). Left and right excitatory supraspinal drives (αL and αR) provide activation for the flexor half-centers (F) of the RGs (ipsi- and contralaterally) and some interneuron populations in the model, as well as for the extensor half-centers (E) (γL and γR ipsilaterally). Two types of feedback (SF-E1 and SF-E2) operating during ipsilateral extension affect (excite) respectively the ipsilateral F and E half-centers, and through V3-E neurons affect contralateral RGs. The SF-E1 feedback depends on the speed of the ipsilateral “belt” (βL or βR) and contributes to extension-to-flexion transition on the ipsilateral side. The SF-E2 feedback activates ipsilateral E half-center and contributes to “weight support” on the ipsilateral side. The ipsilateral excitatory drives (αL and αR) suppress (reduce) the effects of all ipsilateral feedback inputs by presynaptic inhibition. B. Model of spinal-transected system. All supraspinal drives (and their suppression of sensory feedback) are eliminated from the schematic shown in A.

The proposed spinal network connectome allowed the previous models to reproduce and explain multiple experimental phenomena observed during fictive and real locomotion, including the speed-dependent expression of different gaits and its supraspinal control 24-26,28-30. The molecular-genetic types of neurons in this connectome and their known and suggested interactions were based on studies of intact and mutant mice obtained by optogenetic methods, most of which are not available for studies in cats. We assume that the neural connectome in the spinal cord of mammals is evolutionarily conserved and findings from studies in mice can be used for understanding the neural control of locomotion in other quadrupedal mammals, such as cats.

Control of spinal locomotor circuits by supraspinal drives and sensory feedback and interactions between them

In our intact model, the frequency of locomotor oscillations (and the speed of locomotion) is primarily controlled by excitatory supraspinal drives to both RGs, particularly to the flexor half-centers (Fig. 4A). These drives simulate the major descending brainstem pathways to spinal neural circuits. Our previous models also suggested that some supraspinal drives activate/inhibit CINs to modify interlimb coordination and perform gait transitions from left-right alternating (walk, trot) to left-right synchronized (gallop, bound) 24-26,28-30. In the present model, this feature is implemented by supraspinal excitatory signals to V3 CINs (Fig. 4A).

As described for the single RG model (Fig. 3B), both left and right RGs start to intrinsically generate rhythmic locomotor activity when the value of supraspinal drives to the flexor half-centers exceeds some threshold. Below this threshold, the RGs can only operate in a state-machine regime that requires an additional excitatory input to the flexor half-centers to initiate flexion. Such additional trigger signals can come from supraspinal structures, such as the motor cortex, or from limb sensory feedback. In the spinal-transected model that lacks supraspinal drives (Fig. 4B), and in the intact model with low drives necessary for a very slow locomotion, both RGs (specifically, their F half-centers) do not receive sufficient excitation and can only operate in a state-machine regime (see Fig. 3D), in which the extension to flexion transition requires additional signals that can be provided by limb sensory feedback.

It is well known that the operation of the spinal network during locomotion is also controlled by inputs from primary afferents originating in muscle spindles (groups Ia and II), Golgi tendon organs (group Ib) and different skin mechanoreceptors (reviewed in 2). Limb sensory feedback controls phase durations and transitions and reinforces extensor activity during stance. We incorporated two types of limb sensory feedback (SF) in our model in a simplified version, both operating during the extensor phase of each limb.

The first feedback, SF-E1, represents an increase of spindle afferent activity from hindlimb flexor muscles as they are stretched with limb extension (an increase of hip angle) during stance 50, which increases with increasing treadmill speed. This feedback directly activates the ipsilateral F half-center (Fig. 4) and promotes the transition from extension to flexion (simulating lift-off and the stance-to-swing transition). The same feedback acting in the model through the V3-E CIN and InE neurons inhibits the contralateral F half-center (Fig. 4) and promotes the transition from flexion to extension in the contralateral RG. The critical role of hip flexor stretch-related feedback for triggering the stance-to-swing (or extension-to-flexion) transition has been confirmed in many studies in cats during real and fictive locomotion 2,51-55. The role of V3 CINs in transmitting primary afferent activity to the contralateral side was supported by a recent mouse study 56.

The second feedback we incorporated in our model, SF-E2, simulates in a simplified form the involvement of force-dependent Ib positive feedback from limb extensor muscles to the ipsilateral extensor half-center, reinforcing extensor activity and weight support during stance (Fig. 4). This effect and the role of Ib feedback from extensor afferents has been demonstrated and described in many studies in cats during real and fictive locomotion 2,57-59.

An important feature of our model is how supraspinal drives and limb sensory feedback interact with each other and with spinal circuits. Specifically, limb sensory feedback from each limb receives presynaptic inhibition from the ipsilateral supraspinal drive (Fig. 4A), which reduces feedback gains in a speed-dependent manner, suppressing the influence of both feedback types. This interaction in the model reflects experimental data on presynaptic and/or direct inhibition of primary afferents by supraspinal signals 31-33,60. In our model, removing supraspinal drives to simulate the effect of a complete spinal transection also eliminates presynaptic inhibition of all inputs from sensory feedback, increasing the influence of both types of sensory feedback on spinal CPG circuits (Fig. 4B).

We believe that intact animals walking on a treadmill use visual cues and supraspinal signals to adjust their speed and maintain a fixed position relative to the external space. In a simplified version (implemented in the model), the movement speed relative to the treadmill belt is mainly controlled by supraspinal drives to the left and right RGs (parameters αL and αR, see Methods) that define the locomotor oscillation frequency. These drives are automatically adjusted to the speed of the simulated treadmill (parameters βL and βR, see Methods). The frequency of oscillations generated by the RGs and the duty factor of the generated pattern are also affected by SF-E1 and SF-E2 acting during the extension phases of each RG, but they are progressively suppressed by supraspinal drives through presynaptic inhibition with increasing speed. Thus, the role of feedback in the control of locomotion decreases with an increase of locomotor frequency defined by the increasing supraspinal drives. In the spinal-transected model, the transition from extension to flexion is fully controlled by limb sensory feedback. The extensor phase duration and the period (and frequency) of oscillations are mainly defined by the rate of SF-E1 increase, which in turn is defined by the speed of the simulated treadmill (parameters βL and βR, see Methods).

Simulation of tied-belt locomotion in intact and spinal cats

We simulated an increase in treadmill speed in the tied-belt condition by the progressive increase of parameters βLR. As stated above, we assumed that intact cats voluntary adjust supraspinal drives to the left and right RGs to the treadmill speed to maintain a fixed position relative to the external space. Therefore, in the intact model, the parameters αL=αR that characterized supraspinal drives were adjusted to correspond to the parameters βLR characterizing treadmill speed. In intact (Fig. 5A) and spinal-transected (Fig. 5B) models, changes in cycle and phase durations were qualitatively similar to experimental data (compare with Figs. 1A,B). Specifically, the oscillation period (cycle duration) decreases with increasing “speed” due to a shortening of the extensor phase with a relatively constant flexor phase duration, so that the duty factor approaches or reaches 0.5. Comparison of changes in cycle and phase durations of the intact and spinal-transected models demonstrates similar dynamics of changes (Fig. 5C), matching corresponding experimental data (Fig. 1C). It is important to note that in the intact model at “slow treadmill speeds” (low βLR) and the corresponding weak supraspinal drives (αL=αR, approximately at or below 0.4), both RGs are unable to intrinsically generate rhythmic activity and operate in a state-machine regime where switching from the extensor to the flexor phase is mainly controlled by SF-1 and SF-2. At higher values of “treadmill speeds” and corresponding supraspinal drives, locomotion begins to be mainly controlled by supraspinal drives to the RGs. In this case, the role of limb sensory feedback in the control of phase durations is reduced because of increased presynaptic inhibition by increasing supraspinal drives (see Fig. 5C). In contrast to the intact model, the transected model has no supraspinal drives to the RGs and can only operate in a state-machine regime at all treadmill speeds. Thus, the transition from extension to flexion and the duration of the extensor phase are entirely controlled by sensory feedback.

Simulation of locomotion on a tied-belt treadmill using intact and transected models.

A and B. Changes in the durations of locomotor period and flexor/stance and extensor/swing phases during simulated tied-belt locomotion using the intact (Fig. 4A) and transected (Fig. 4B) models with an increasing simulated treadmill speed. C. Superimposed curves from panels A and B to highlight differences.

Simulation of split-belt locomotion in intact and spinal cats

For split-belt locomotion, we simulated the speed of the “slow” RG with a fix value of βL=0.4 and the “fast” RG with βR progressively increasing from 0.5 to 1.0. Again, based on our assumption above, the supraspinal drives to the left and right RGs in the intact model were adjusted to the corresponding speeds of the treadmill belts. We set a constant drive αL=0.4 to the left RG and progressively increased the drive αR to the right RG from 0.5 to 1.0. We show that in the intact model (Fig. 6A), despite differences in the “speeds” of left and right RGs, the left and right oscillation periods are equal, which corresponds to experimental data (Fig. 2A).

Simulation of locomotion on a split-belt treadmill using intact and transected models.

A. Changes in the durations of locomotor period and flexor/stance and extensor/swing phases for the left (slow) and right (fast) sides during split-belt treadmill locomotion using the intact model (Fig. 4A). B. Changes in the same characteristics for the left (slow) and right (fast) sides during simulation of split-belt treadmill locomotion using the transected model (Fig . 4B). In both cases, the speed of the simulated left (slow) belt was constant (βL=0.4) while the speed of the simulated right belt (βR) changed from 0.5 to 1.0 with 0.1 increments.

In the spinal-transected model (Fig. 6B), changes in cycle and phase durations are similar to those in the intact model and correspond to experimental data (Fig. 2B). The oscillation periods on the left and right sides are equal, and changes in cycle and phase durations on the left slow side do not change much with a progressive increase of βR from 0.5 to 1.0. The flexor phase duration of the fast RG increases with increasing βR, which compensates for the decrease of the extensor phase duration.

An increase of flexion phase duration on the fast side with an increase of speed on that side in the spinal-transected model is much steeper than in the intact model (compare Fig. 6B with Fig. 6A) which reproduces corresponding experimental data (see Fig. 2B vs. Fig. 2A). The mechanism of this increase in the spinal case (spinal-transected model) is the following. With the increase of speed, the SF-E1 starts accumulating a tonic component acting during flexion. This accumulated tonic activity provides a direct excitation to the ipsilateral F half-center during flexion, which is compensated by inhibition from the contralateral SF-E1 that inhibits the ipsilateral F half-center during the flexor phase (via contralateral V3-E CIN and ipsilateral InE neuron, see Fig. 4). In tied-belt locomotion (Fig. 5), this contralateral inhibition acts in balance with the accumulated ipsilateral excitation from SF-E1 and limits excitation of each F half-center during flexion, keeping both RGs within the flexor-driven operating regime with relatively constant flexor phase durations (as in Fig 3D, middle part of the graph). In split-belt locomotion, this contralateral inhibition from the slow (left) side remains relatively weak (it is fixed at a value corresponding to the constant slow left-side speed and does not change with the increasing speed on the right fast side), whereas the accumulated excitation on the fast side continues increasing with the speed). Therefore, the F half-center on the fast side becomes overexcited and the RG on the fast side starts operating in the classical half-center regime (like in Fig. 3D, right part of the graph) with an increasing flexor phase duration as the speed on the fast side increases (Fig. 6B). A similar mechanism operates in the intact split-belt case (Fig. 6A), but the resultant increase of the flexor phase duration on the fast side is much weaker due to presynaptic inhibition of left and right SF-E1 feedback.

Effects of feedback removal on treadmill locomotion

We used our model to simulate the effects of removing limb sensory feedback. We hope that the results of these simulations allow us to formulate modeling predictions which can be tested in future experiments providing an additional validation to our model. We specifically focused on simulation of separate removal of SF-E1 feedback, activated with limb extension and involved in the stance-to-swing transition, or SF-E2 feedback involved in the generation of extensor activity and weight support during stance, as well as removal of both feedback types.

In the transected model, removal of SF-E1 (with or without SF-E2) prevented locomotion (not shown). Without supraspinal drive and SF-E1, both RGs could not switch to the flexor phase and generate locomotor activity. Following removal of SF-E2 in the transected model, sufficient extensor activity (necessary for “weight support”) could not be developed, leading to an abnormal locomotion with extremely short extensor phase durations (and oscillation periods), with an inability to adjust to treadmill speed (Fig. 7). In this case locomotion could be recovered by adding an additional activation of extensor half-centers during stance phases (not shown).

Simulation of the effect of removing SF-E2 feedback in the transected model during simulated tied-belt locomotion.

Changes in the durations of locomotor period and flexor/stance and extensor/swing phases during simulated tied-belt locomotion using the transected model (Fig. 4B) after removal of SF-E2 feedback.

The described above effects of removing limb sensory feedback on locomotion in the spinal-transected model were expected. However, what happens after removal limb sensory feedback in the intact model, in which supraspinal drives is present? Removing SF-E1 in the intact model produced interesting results. It considerably increased the duration of the extensor phase (Fig. 8A). However, the model only demonstrated oscillatory activity starting at some moderate speed, when βLR 0.5. This required the corresponding supraspinal drives to the RGs to be sufficient (αL=αR 0.5) to operate in the flexor-driven regime and produce phase transitions without any contributions from SF-E1. In contrast, removing SF-E2 in the intact model reduced the activity of the extensor half-centers, shortening extensor phase durations but preserved oscillations in the full range of considered “speeds” (Fig. 8B). Therefore, the removal of SF-E2 produced an opposite effect on extensor phase durations compared with removing SF-E1. Finally, removing both feedback types in the intact model shifted the start of oscillatory activity to the left compared to removal of only SF-E1, allowing oscillations to start at βLR 0.35 (Fig. 8C). This allows the prediction that below 0.35 m/s, cats with diminished limb sensory feedback can only perform locomotion with patterned step-by-step supraspinal signals to produce phase transitions. This can also explain why intact cats do not usually walk consistently on a treadmill with a speed at or below 0.3 m/s (see Fig. 1A).

Simulation of the effect of removing SF-E1, or SF-E2, or both feedback types in the intact model during simulated tied-belt locomotion.

Changes in the durations of locomotor period and flexor/stance and extensor/swing phases during simulated tied-belt locomotion using intact model (Fig. 4A) with an increasing simulated treadmill speed after removal of only SF-E1 feedback (A), only SF-E2 feedback (B), and both feedback types (C).

Discussion

Operation regimes of the spinal locomotor network

Despite decades of research, a commonly accepted definition of the spinal locomotor CPG has not been formulated, and different authors used this term in relation to different entities depending on the context. For instance, the term CPG has been used to designate spinal circuits controlling and coordinating rhythmic movements of all limbs, or controlling only rhythmic movements of a single limb, or even a single joint 4,6,7,15-19,36. Here, we use the term CPG for the spinal circuitry controlling and coordinating all limbs, and the term RG (rhythm generator) for the relatively independent part of the CPG that controls rhythmic movements of a single limb. We consider the CPG as a group of RGs, each controlling a single limb and interacting with each other through commissural and/or propriospinal pathways or circuits 3,24-30.

Regardless of the exact name (RG or unit burst generator), a common view is that the function of the RG is to generate locomotor-like rhythmic bursting consisting of two major phases, flexion and extension. In each phase, the controlled limb operates in the corresponding functional state, defining the contraction and relaxation of specific sets of muscles. Thus, the RG is not just an oscillator, but rather a system with a dual function: it functions as a state-machine 61-63, defining the state and operation of the controlled biomechanical system in each phase independent of the exact transition mechanisms, and as an oscillator defining the mechanisms and timing of transitions between states. These transitions may be fully defined by internal properties of the RG and/or require a contribution from external inputs, such as supraspinal signals or somatosensory feedback, that control or adjust the timing of transitions.

Here, we described a relatively simple model of a half-center RG. We identified three operation regimes: state-machine, flexor-driven, and classical half-center. In an intact system, at slow speeds ( 0.35 m/s), the spinal network operates in a state regime and requires external inputs for phase transitions, which can come from limb sensory feedback and/or volitional inputs (e.g. from the motor cortex). At higher speeds and greater supraspinal drives, the spinal network switches to a flexor-driven regime before transitioning to a classical half center regime at higher speeds/drives. Following spinal transection, the spinal network can only operate in the state regime and entirely depends on limb sensory feedback for phase transitions.

Our modeling results also resolve a potential contradiction between the classical half-center and flexor-driven concepts of spinal RG operation. According to the classical half-center concept, both flexor and extensor half-centers are necessary for generating rhythmic activity 7,15,16. An alternative, flexor-driven concept was proposed by Pearson and Duysens (“swing generator model”), where the flexor half-center is intrinsically rhythmic, unlike the extensor half-center, which exhibits rhythmic activity due to rhythmic inhibition from the flexor half-center 64,65. Here we follow our previous modeling studies 13,48 and show that these concepts do not contradict each other, but rather relate to different regimes of RG operation. With an increase of excitatory drive to an RG, the rhythmogenic mechanism changes from a flexor-driven to a classical half-center mechanism (Fig. 3). This transition explains the increase of the flexor phase duration of the fast hindlimb or RG during split-belt treadmill locomotion in intact and spinal cats (Figs. 2 and 6).

Intrinsic spinal rhythm-generating versus reflex-based mechanisms of locomotion

Charles Sherrington proposed that locomotion was generated by chains of reflexes, mostly by studying locomotion in decerebrate and spinal cats 21,22, whereas Thomas Graham Brown clearly showed that an intrinsic spinal mechanism could generate rhythmic alternating activity 15, which later became the CPG concept. Our results support both concepts, depending on locomotor speed and the state of the animal. Based on our simulations, we suggest that in spinal cats at any speed and in intact cats at low speeds, the spinal network operates in a state-machine regime and requires sensory feedback to locomote, consistent with Sherrington’s viewpoint. In contrast, at higher (moderate) speeds ( 0.4 m/s), when supraspinal drive is sufficient, the spinal network can intrinsically generate the basic locomotor activity controlling locomotion, which supports Brown’s concept.

Our results also suggest an important difference in the control of slow exploratory and faster escape locomotion 66-69 . Based on our predictions, slow (conditionally exploratory) locomotion is not “automatic”, but requires volitional (e.g. cortical) signals to trigger step-by-step phase transitions because the spinal network operates in a state-machine regime. In contrast, locomotion at moderate to high speeds (conditionally escape locomotion) occurs automatically under the control of spinal rhythm-generating circuits receiving supraspinal drives that define locomotor speed, unless voluntary modifications or precise stepping are required to navigate complex terrain 70,71.

We also used our model to simulate and predict the effects of removing limb sensory feedback on locomotion. As expected, the spinal-transected model failed to generate locomotion without SF-E1 or normal locomotion without SF-E2 because the spinal network operates in a state-machine regime (Fig. 7). The essential role of somatosensory feedback after spinal cord injury is well known 2,72-76. However, in the intact model, limb sensory feedback, particularly from hip flexor afferents (SF-E1), is also required for locomotion and phase transitions at slow speeds to compensate for low supraspinal drives (Fig. 8).

Another important implication of our results relates to the recovery of walking in movement disorders, where the recovered pattern is generally very slow. For example, in people with spinal cord injury, the recovered walking pattern is generally less than 0.1 m/s and completely lacks automaticity 77-79. Based on our predictions, because the spinal locomotor network operates in a state-machine regime at these slow speeds, subjects need volition, additional external drive (e.g., epidural spinal cord stimulation) or to make use of limb sensory feedback by changing their posture to perform phase transitions.

Model limitations and future directions

In this study, we developed a simplified model of the neural control of cat hindlimb locomotion in tied- and split-belt conditions to simulate and compare locomotion of intact and spinal cats. The major limitation of the present model is the lack of biomechanical elements simulating multi-joint limbs and muscles. Our study and analysis were specifically focused on the proposed organization and operation of spinal CPG circuits, including left-right interactions within these circuits, and their control by supraspinal drives and limb sensory feedback. Another important limitation is that we did not consider possible mechanisms of plasticity following spinal transection. We know that the spinal connectome and sensorimotor interactions change after spinal cord injury 76,80. The precise nature of these changes is not well understood. Nevertheless, we believe that changes (increase) in the gain of sensory feedback, which in the model results from eliminating presynaptic inhibition, in a real situation may include a considerable contribution from plastic changes during the recovery period. However, even with the account of possible plastic changes during the recovery period, the functional role of the increased gain of sensory feedback in spinal-transected animals remains the same independent of the exact mechanisms for this increase (release from presynaptic inhibition, plastic changes during recovery, currently unknown mechanisms, or any combination of the above).

Although the model was based on simplifications and some assumptions, it reproduced and provided explanations for experimental results, such as the main locomotor characteristics (cycle and phase durations) at different speeds and left-right speed differences during tied-belt and split-belt locomotion. We are currently applying our model to simulate cat locomotion following incomplete spinal cord injury (lateral hemisection) during tied-belt and split-belt locomotion based on our recent experimental data 81-83. Using our model, this will allow us to remove supraspinal drives on the hemisected side and determine how the spinal network controls locomotion. We also plan to incorporate a full biomechanical model of the limbs to investigate the neuromechanical control of quadrupedal locomotion and its recovery following incomplete spinal cord injury 84,85.

Methods

Experimental data

For this modeling study, we used previously published data obtained from intact and spinal cats during tied-belt (equal left-right speeds) and split-belt (different left-right speeds) locomotion 11-13. No new animals were used here. In those studies, all procedures were approved by the Animal Care Committee of the Université de Sherbrooke (Protocol 442-18) in accordance with policies and directives of the Canadian Council on Animal Care. Intact cats performed quadrupedal locomotion whereas spinal cats performed hindlimb-only locomotion with the forelimbs on a stationary platform. Intact and spinal cats performed tied-belt locomotion from 0.4 to 1.0 m/s and from 0.1 to 1.0 m/s, respectively, with 0.1 m/s increments. As stated earlier, intact cats cannot perform consistent quadrupedal tied-belt locomotion at or below 0.3 m/s. For split-belt locomotion, the slow belt (left) was 0.4 m/s while the fast belt (right) stepped at speeds of 0.5 to 1.0 m/s in 0.1 m/s increments.

Modeling formalism and model parameters

In our model, we considered spinal circuits as a network of interacting neural populations. Each population is described by an activity-based neuron model (and is sometimes called a neuron in the text), in which the dependent variable V represents an average population voltage and the output f(V) (0 ≤ f(V) 1) represents the average or integrated population activity at the corresponding average voltage 87-88. This description allows an explicit representation of ionic currents, specifically of the persistent (slowly inactivating) sodium current INaP (Rubin et al., 2009), which was proposed to be responsible for generating intrinsic bursting activity in the spinal cord 43,46,47,89-93. Assuming that neurons within each population switch between silence and active spiking in a generally synchronized way, the dynamics of the average voltages are represented within a conductance-based framework used for a single neuron description, but without fast membrane currents responsible for spiking activity.

The dynamics of V for flexor and extensor half-centers, considered as INaP-dependent conditional bursters, is described by the following differential equation:

For other (non-bursting) populations, V is described as:

The output function f(V) transvers V to the integrated population output and is defined as follows:

In Eqs. (1) and (2), C is the membrane capacitance, INaP is the persistent sodium current, IL is the leakage current, ISynE and ISynI represent excitatory and inhibitory synaptic currents, respectively. The leakage current was described as:

where gL is the leakage conductance and EL represents the leakage reversal potential. The persistent sodium current in the flexor and extensor half-centers is described as:

Where is the maximal conductance and ENa is the sodium reversal potential. Its voltage-dependent activation, m(V), is instantaneous and its steady state is described as follows:

The slow INaP inactivation was modeled by the following differential equation:

where h(V) represents inactivation steady state and τh(V) is the inactivation time constant with maximal value τmax.

In Eqs. (6), (8), and (9), V1/2 and k represent half-voltage and slope of the corresponding variables (m, h, and τ). Excitatory and inhibitory synaptic currents (ISynE and ISynI) for population i were described by:

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

The weights of connections wji between populations for the main model (Fig. 4) are shown in Table 1. Di in Eq. (10) for the main model (Fig. 4), represents the total excitatory drive to population i:

where αipsi and αcontra are the ipsi- and contralateral excitatory drives, respectively, which depending on the side, represent left αL and right αR drives, respectively; parameter γ represents the constant drive to the ipsilateral (left or right) extensor half-center; , and define the weights of these drives to population i.

Connection weights.

and in Eq. (10) for the main model (Fig. 4), define the effect of ipsilateral sensory feedback E1 (or SF-E1, see Fig. 4) and E2 (or SF-E2), respectively, to population i. The gain of each feedback to population i (and ) is suppressed (reduced) by the ipsilateral drive αipsi (αL or αR depending on the side):

Where and are the weights of presynaptic inhibition by αipsi of feedback inputs E1 and E2, respectively, to population i (see Fig. 4A).

E1 feedback (SF-E1) represents an increase of activity of length-dependent hip flexor afferents during limb extension. In our model, it is described as:

where kE1 is the parameter of E1 feedback ; βipsi is the parameter characterizing the speed of the ipsilateral treadmill belt (left βL or right βR), f(VE) is the output of the ipsilateral E half-center, thr is the threshold, t0 is the starting time of the ipsilateral extensor burst and t is time (both in ms).

E2 feedback (SF-E2) that represents excitatory feedback from load-dependent afferents of limb extensor muscles to the extensor half-center during extension, described as follows:

where E2o and kE2 are the parameter of SF-E2.

Variables αL and αR characterize left and right supraspinal drives to RGs and define the frequency of locomotor oscillations. The variables βL and βR define the speeds of the left and right treadmill belts. Intermediate coefficients were used to make correspondence between α and β values and between their values and the speed of real treadmill belts in the experiments. For modeling of intact locomotion, we changed left and right α values to simulate “voluntary” locomotion in the intact model via a corresponding adjustment of supraspinal drives. For modeling of locomotion in the spinal-transected model, supraspinal drives were set to 0) and we only manipulated left and right β values simulating speeds of treadmill belts.

The following values of model parameters were used for the main model (Fig. 4): C = 20 pF; gL = 4 nS for RG half-centers and gL = 1 nS for all other neurons; ; gsynE = gsynI = 1 nS; EL = −60.0 mV; ENa = 50.0 mV; EEsynE = 10 mV; EEsynI = 75 mV; Vthr = 50 mV; Vmax = 0 mV; V1/2,m = 40.0 mV; km = −6 mV; V1/2,h = 45.0 mV; kh = 4 mV; τmax = 500 ms; V1/2,τ = 45 mV; kτ = 20 mV; E2o = 2.5; kE1 = 3.12 ∙ 10-3 ms-1; kE2 = 0.5 ∙ 10-3 ms-1; τ1 = 500 ms; τ2 = 300 ms; thr =0.1. All connection weights for the main model (Fig. 4) are specified in Table 1.

In the models shown in Fig. 3: EL = −64.0 mV; w1 = w3 = 3.6; w2 = 3.0; w4 = 5.0; Drive-E = 1,2. All other parameters are the same as in the main model, listed above.

Simulations, Data analysis and availability

All simulations were performed using the custom neural simulation package NSM 2.5.7. The simulation package was previously used for models of spinal circuits 7,25-27,30,46,47. Differential equations were solved using the exponential Euler integration method with a step size of 0.1 ms. Simulation results were saved as ASCII files and represented the output functions for half-centers recorded with a precision of 0.1 ms.

The simulation results were processed using custom Matlab scripts (The Mathworks, Inc., Matlab 2023b). To assess model behavior, the activities of flexor and extensor half-centers were used to determine the onsets and offsets of flexor and extensor bursts and to calculate flexor and extensor phase durations and oscillation period. The timing of onsets and offsets of flexor and extensor bursts was determined at a threshold level of 0.05. The oscillation period was defined as the duration between two consecutive burst onsets in the left extensor half-center. The flexor and extensor phase durations and oscillation period were averaged over the duration of simulation for each value of the parameter α. Duration of individual simulations depended on the value of parameter α to robustly estimate average values of burst durations and oscillation period. For each α value, we omitted the first two transitional cycles to allow stabilization of model variables. Flexor and extensor durations and oscillation period were plotted against the parameters α and β.

The simulation package NSM 2.5.7, the model configuration file necessary to create and run simulations, and the custom Matlab scripts are available at https://github.com/RybakLab/nsm.

Grant support

This study was supported by NIH/NINDS grant R01 NS110550.

Autor contributions

I.A.R. suggested the concept. A.F., B.I.P. and N.A.S. elaborated the details of the work. A.F. provided all experimental data. I.A.R. and N.A.S. developed the basic model. S.N.M. provided and supported simulation tools. N.A.S., I.A.R. and S.N.M. tuned model parameters. I.A.R., A.F. and B.I.P. interpreted the results. I.A.R., N.A.S., and A.F. prepared the figures. I.A.R. and A.F. drafted the manuscript. All authors discussed, reviewed, edited, and approved the final version of the manuscript.

Competing interests

The authors declare no competing interests.