The locomotor gait in limbed animals is defined by the left-right leg coordination and locomotor speed. Coordination between left and right neural activities in the spinal cord controlling left and right legs is provided by commissural interneurons (CINs). Several CIN types have been genetically identified, including the excitatory V3 and excitatory and inhibitory V0 types. Recent studies demonstrated that genetic elimination of all V0 CINs caused switching from a normal left-right alternating activity to a left-right synchronized “hopping” pattern. Furthermore, ablation of only the inhibitory V0 CINs (V0D subtype) resulted in a lack of left-right alternation at low locomotor frequencies and retaining this alternation at high frequencies, whereas selective ablation of the excitatory V0 neurons (V0V subtype) maintained the left–right alternation at low frequencies and switched to a hopping pattern at high frequencies. To analyze these findings, we developed a simplified mathematical model of neural circuits consisting of four pacemaker neurons representing left and right, flexor and extensor rhythm-generating centers interacting via commissural pathways representing V3, V0D, and V0V CINs. The locomotor frequency was controlled by a parameter defining the excitation of neurons and commissural pathways mimicking the effects of N-methyl-D-aspartate on locomotor frequency in isolated rodent spinal cord preparations. The model demonstrated a typical left-right alternating pattern under control conditions, switching to a hopping activity at any frequency after removing both V0 connections, a synchronized pattern at low frequencies with alternation at high frequencies after removing only V0D connections, and an alternating pattern at low frequencies with hopping at high frequencies after removing only V0V connections. We used bifurcation theory and fast-slow decomposition methods to analyze network behavior in the above regimes and transitions between them. The model reproduced, and suggested explanation for, a series of experimental phenomena and generated predictions available for experimental testing.
Central Pattern Generators (CPGs) are neural networks that can produce organized rhythmic motor activities in the absence of rhythmic inputs and feedbacks from other parts of the nervous system. The CPGs that generate basic rhythmic locomotor activities and control locomotion in vertebrates are located in the spinal cord [1–11]. It appears that each limb in mammals is controlled by a separate CPG because cats with chronic thoracic spinal lesions were shown to step on a split treadmill with different speed for the left and right hind limbs [12]. The gait of locomotion is defined by the coordinated limb movements and hence by the phase relationships between rhythmic patterns generated by CPGs controlling each limb. In turn, these relationships are defined by neural circuits within the spinal cord providing direct or indirect interactions between the CPGs. Several computational models of neural circuits were proposed to reproduce and explain the left-right and segmental coordination during locomotion in lamprey [13–16] and left-right hindlimb and quadrupedal limb coordination in the salamander [17–19] and mammals [20–22]. Some computational studies focused on gait transitions analyzing them as bifurcation phenomena [23–25]. Spinal circuits in some models were simulated and analyzed as coupled nonlinear oscillators [26–28]. In these models, each oscillator represented an independent rhythm-generating center (simulated as a single neuron or a neural population), and the main goal of these models was to study how inter-oscillator couplings could affect the synchronization and the phase lags between rhythm-generating centers [18]. It is considered that the locomotor gait and pattern in these models mostly depend on the topology of interactions and properties of neurons involved, rather than on the local rhythm-generation mechanisms [18]. However, even in lower vertebrates such as the lamprey, it is extremely difficult to identify and characterize all neurons and neural circuits that are involved in such coupling/interactions during locomotion. It is generally considered that phase relationships between the activities of corresponding neurons on the left and right sides of the spinal cord are provided by so-called commissural interneurons (CINs), whose axons cross the midline and innervate neurons on the contralateral side of the cord [6,29–32]. There are different types of inhibitory and excitatory CINs, which are involved in the left-right alternating (e. g. , during normal walking) or left-right synchronized (e. g. , during hopping or galloping) locomotor activities [6–8,29–32]. Significant progress in understanding the functional and structural organization of spinal circuits has been achieved due to recently developed combinations of genetic, molecular, and developmental techniques. Several populations of neurons were found to be derived from genetically distinct populations of embryonic neurons in the spinal cord of adult mice [33–36]. Some of these genetically characterized neurons were identified as CINs, including the V0 and V3 interneurons [6,30,33,37,38]. The V0 population contains both excitatory and inhibitory CINs that both contribute to alternating left and right neuronal activities [6,30,37]. The V3 neurons are excitatory and mediate interactions promoting left-right synchronized hopping behavior [38]. The V0 population contains genetically distinct subpopulations of CINs: V0V and V0D. The inhibitory V0D neurons constitute about two-thirds of the V0 population and the excitatory V0V neurons constitute about one-third [6,30,39]. The specific functional roles of V0V and V0D subpopulations in coordinating left and right locomotor activities have been recently studied using isolated spinal cord preparations from genetically transformed mice with either V0V, or V0D, or both V0 subpopulations knocked out [37]. The spinal cord was isolated in vitro to ensure that supraspinal structures were not participating in rhythm generation. The specific contribution of each subpopulation to the left-right coordination has been found to depend on the frequency of locomotor oscillations. Locomotor oscillations were induced by application of a mixture of N-methyl-D-aspartate (NMDA) and 5-hydroxytryptamine (5-HT). These oscillations were characterized by alternating activities recorded from left and right flexor-related (LL2 and RL2) and left and right extensor-related (LL5 and RL5) lumbar roots. The speed of locomotion (oscillation frequency) was regulated by the NMDA concentration [37,40]. Measuring the phase differences between rhythmic activities recorded from LL2 and RL2 roots and from LL5 and RL5 roots allowed for identification of an alternating or synchronized pattern. The results of these studies can be summarized as follows [37]: (i) spinal cord preparations from control mice showed left-right alternating behavior at all locomotor frequencies; (ii) preparations from the mice with genetically ablated V0 populations (both V0V and V0D types) showed synchronized (hopping) activity at all locomotor frequencies; (iii) preparations from the mice with only V0V populations ablated maintained left-right alternation at low frequencies but switched to hopping gait at high frequencies; (iv) preparations from the mice with selectively ablated V0D populations showed hopping at low frequencies and switched to alternation at high frequencies. This model represents a simplified version of the large-scale model [41], which was used for computer simulation of spinal circuits involved in rhythm generation and frequency-dependent left-right CIN-mediated interactions, but was too complicated to use methods of qualitative theory of differential equations and dynamical systems. Therefore, in this study, we formulated and analyzed a simplified mathematical model. The model consisted of four pacemaker neurons representing left and right, flexor and extensor rhythm-generating centers interacting via commissural pathways representing V3, V0D, and V0V CINs. We used this model to investigate different regimes of behavior of the intact circuitry, and following a removal of commissural interactions mediated by V0D, or V0V, or both V0 CIN pathways to mimic recent experimental studies of Talpalar et al. [37]. Then, using methods of bifurcation theory and fast-slow decomposition, we analyzed the essential properties of this system under different conditions and propose plausible explanations for the experimentally observed behavioral transformations.
Most of the recent models of mammalian spinal circuits and locomotor CPGs were based on interacting populations of tens or hundreds of spiking neurons modeled in the Hodgkin-Huxley style and described by 5–10 differential equations per neuron [41–46]. Such large networks of spiking neurons are too complex to be analyzed by traditional mathematical approaches, and hence certain simplifications need to be applied. Such a simplification has been previously proposed in the analysis of the model of the mammalian respiratory CPG [47,48]. In these studies, each population of spiking neurons was described by a single, non-spiking neuron model, in which the voltage variable represented an average voltage for the population and the output activity was described as a nonlinear function of voltage, f (V) [47–49]. It was established that such a simplification can provide a reasonably accurate description of population activity, including transitions between quiescent (silent) and active states in population dynamics [47]. This simplification was used in our present model. The other simplification of the current study was that, in contrast to the previous multi-scale models [42–45], we focused only on rhythm-generating circuits and did not consider motoneurons and interneurons (such as Ia and Renshaw cells) that are not critically involved in rhythmogenesis and left-right coordination. Therefore, the reduced neural circuit that we analyzed included two (left and right) rhythm generators (RGs) bilaterally interconnected via several commissural pathways. Each RG consisted of one flexor (left, LF, or right, RF) and one extensor (left, LE, or right, RE) centers reciprocally inhibiting each other via inhibitory interneurons (Ini1 and Ini2, respectively) (Fig 1A). Three major commissural pathways were considered: (1) the excitatory V3 CINs mediate mutual excitation between the left and right flexor centers (LF and RF); (2) the inhibitory V0D CINs mediate mutual inhibition between these centers; and (3) the excitatory V0V CINs that were also shown to contribute to mutual inhibition between these centers and promote left-right alternation. To provide mutual inhibition, the excitatory V0V CINs should either receive inputs from the ipsilateral flexor centers and affect the contralateral flexor centers via inhibitory interneurons (as suggested in [37,41]) or mediate a crisscross excitation from each extensor center to the contralateral flexor center as shown in Fig 1A (see also [41]). Because no recordings from V0V neurons have been made so far, it is difficult to determine which of the two V0V pathways suggested is more realistic or whether both these pathways are present in the spinal cord (see discussion in [41]). Preliminary large-scale simulations have shown that both pathways lead to qualitatively similar behavior [41]. Therefore, in this study, we only focused on the V0V CIN pathways providing the direct crisscross excitation from each extensor center to the contralateral flexor center (Fig 1A). The final simplification concerned all interneurons (i. e. CINs, Ini1, and Ini2), which, for simplicity, were removed and replaced by the corresponding direct synaptic connections as shown in Fig 1B. Therefore, the final simplified model under investigation consisted of four centers (LF, LE, RF, and RE) with the unilateral inhibitory interactions between flexor and extensor centers on each side via Ini1 and Ini2 pathways and the bilateral V3, V0D, and V0V commissural interactions between the left and right centers as shown in Fig 1B. Each center is described using a reduced activity-based non-spiking neuron model with only two dynamical variables. This model was adapted from a reduced model of the population of intrinsically bursting neurons in the pre-Bötzinger complex developed by Rubin et al. [47,48]. Endogenous bursting in these neurons was suggested to involve the persistent (slowly inactivating) sodium current, INaP, first described by Butera et al. [50] (Model 1), with burst termination based on the slow inactivation of this current. This model was shown able to generate busting activity with a wide range of a parameters defining neuronal excitation and burst frequencies, and exhibited the correct change (reduction) of burst amplitude when the frequency increased [47,48]. This model has been successfully implemented in several previous models of the locomotor CPG [41–46,51,52]. The persistent sodium current was indeed found in spinal interneurons [52–56] and its blockade terminated rhythm generation in the rat spinal cord [53]. Similar to the Rubin et al. models [47,48], our formulation for each center includes an explicit representation of INaP: C⋅V˙=−INaP−IL, (1) where C is the membrane capacitance, INaP represents the persistent sodium current, and IL is the leak current. The currents in Eq (1) are described as follows: INaP=g¯NaP⋅mNaP∞ (V) ⋅hNaP⋅ (V−ENa); (2) IL=g¯L⋅ (V−EL), (3) where g¯NaP and g¯L are the maximal conductances of the persistent sodium and leak channels, and ENa and EL are the reversal potentials for sodium and leak currents, respectively. hNaP is the INaP inactivation gating variable, and mNaP∞ (V) represents voltage-dependent steady state of INaP activation: The INaP activation, mNaP∞ (V), is considered instantaneous; the INaP inactivation is slow and hNaP represents the “slow” dynamical variable in this model described as follows: τhNaP (V) ⋅h˙NaP=hNaP∞ (V) −hNaP, (5) where hNaP∞ (V) and τhNaP (V) represent the voltage-dependent steady state and time constant for inactivation, respectively: In Eqs (4), (6) and (7), VxhNaP and kxNaP for x ∈ {m, h, τ) represent the gating variable’s half-activation voltage and slope, respectively. The values of all parameters used are provided in Table 1. The output of each neuron (center) describing its activity (normalized firing rate) is represented by a piecewise linear function f (V), changing between 0 and 1 such that: f (V) ={0, ifV<Vmin (V−Vmin) / (Vmax−Vmin), ifVmin≤V<Vmax1, ifV≥Vmax (8) where Vmin and Vmax define the voltages at which threshold and saturation are reached, respectively. The persistent (slowly inactivating) sodium current, INaP, provides a neuron with intrinsic rhythmicity. The behavior of this neuron can be considered using a (V, hNaP) phase plane (Fig 2A). The large time constant of INaP inactivation makes the time scales of the two dynamical variables, V and hNaP, substantially different which allows us to apply the fast-slow decomposition technique for qualitative analysis of neuronal behavior. In the phase plane shown in Fig 2A, V- and hNaP-nullclines (red and black, respectively) are calculated by setting the right hand side of Eqs (1) and (6) equal to zero. EL was increased to mimic higher levels of excitation and this was reflected by a downward movement of the V-nullclines on the phase plane. The band of nullclines where intrinsic bursting activity occurs is shaded gray. Within this area the hNaP-nullcline intersects with the V-nullcline at a point between the V-nullcline knees. This corresponds to an unstable critical point and allows for the emergence of a stable limit cycle. This limit cycle is shown in the (V, hNaP) phase plane by a neuron’s trajectory with the direction indicated by arrows. The neuron’s image point moves slowly along the left and right segments of the V-nullcline (slow timescale, silent and active phases) and quickly jumps between them (fast timescale, switching between phases) showing the behavior typical for a relaxation oscillator. When a neuron is silent its image point travels up the left branch of the V-nullcline until reaching the left knee causing it to jump toward the right branch of the V-nullcline in a movement representing activation of a neuron. Once the neuron is active, the corresponding trajectory will travel down the right branch of the V-nullcline until it reaches the right knee where it will then jump back to the left branch and the neuron will once again be silent. Based on the above description, each isolated neuron can be in one of three regimes: silent, bursting (oscillatory), and tonic (constant activity) (Fig 2A) which depends on the level of neuronal excitation that, in the case of an isolated neuron, is defined by the leak reversal potential, EL. As shown in Fig 2B, bursting in the isolated neuron occurs when EL ∈ (−62. 7, −54. 2) mV (gray area). When EL is less than -62. 7 mV, the hNaP-nullcline in Fig 2A intersects with the V-nullcline at a position to the left of its left knee (all nullclines that exist in the area marked" silence" ). This creates a stable critical point on the inactive (left) branch of the V-nullcline where the neuron' s position represents an inactive (silent) state. When EL exceeds -54. 2 mV, the V-nullcline in Fig 2A loses its cubic-like shape. This creates a stable critical point on the V-nullcline (all nullclines that exist in the area marked" tonic" in Fig 2A) and the neuron therefore demonstrates a constant level of activity, which is referred to as tonically active. Therefore, with an increase of excitation defined by EL the neuron can potentially go through silent to bursting and then to tonically active regimes. The model considered here consists of two (left, L, and right, R) bilaterally connected rhythm generators, each of which represents one flexor (F) and extensor (E) center (Fig 1). As described above, each center can be in silent, bursting, or tonic regimes, depending on the level of its excitation. However, we assume that, at least in the conditions of rhythmic bursting (fictive locomotion) evoked in the isolated rodent spinal cord by neuroactive drugs, the rhythm generators are asymmetric so that only the flexor centers (left, LF, and right, RF) operate in the bursting mode, whereas the extensor centers (left, LE, and right, RE), if isolated, operate in the regime of tonic activity (this issue is specifically addressed in the Discussion). The excitabilities of the flexor and extensor centers (defined by their leak reversal potentials, ELF and ELE, respectively) in the model are regulated by a parameter α. By introducing this parameter we did not intend to perform a biologically realistic simulation of the cellular mechanisms responsible for NMDA-induced bursting in spinal interneurons, which are not well understood. Instead, this parameter was used to simulate the general effect of NMDA, whose application in the isolated spinal cord preparations leads to an increase in neuronal activity and locomotor frequency [37,40]. To simulate the effect of the drug on the excitabilities of both centers, we suggest that: ELF=ELFO⋅ (1−βE⋅α); (9) ELE=ELEO⋅ (1−βE⋅α), (10) where ELFO and ELEO define the basal values of the leak reversal potentials of the corresponding centers (at α = 0), and βE defines the rate of their increase with α. To introduce the asymmetry suggested above, we set ELFO close to the border between silence and bursting and ELEO slightly above the bursting area (Fig 2B) so that with increase of α (in the uncoupled case), the flexor centers would always operate in the bursting mode, whereas the extensor center would be in the tonically active regime. The bilaterally symmetrical network of four interacting centers with their interconnections (Fig 1B) can be described as follows: C⋅V˙i=−INaPi−ILi−ISynEi−ISynIi; (11) INaP=g¯NaP⋅mNaP∞ (V) ⋅hNaP⋅ (V−ENa); (12) ILi=g¯L⋅ (Vi−ELi), (13) where i (i = 1,2, 3,4) is the index corresponding to the center’s number as shown in Fig 1B. Currents ISynEi and ISynIi define, respectively, the excitatory and inhibitory synaptic interactions between the centers: ISynE1= (a3⋅f (V3) +aV⋅f (V4) ) ⋅g¯SynE⋅ (V1−ESynE); (14) ISynE3= (a3⋅f (V1) +aV⋅f (V2) ) ⋅g¯SynE⋅ (V3−ESynE); (15) ISynI1= (aD⋅f (V3) +b2⋅f (V2) ) ⋅g¯SynI⋅ (V1−ESynI); (16) ISynI2=b1⋅f (V1) ⋅g¯SynI⋅ (V2−ESynI); (17) ISynI3= (aD⋅f (V1) +b2⋅f (V4) ) ⋅g¯SynI⋅ (V3−ESynI); (18) ISynI4=b1⋅f (V3) ⋅g¯SynI⋅ (V4−ESynI), (19) where g¯SynE and g¯SynI are the maximal conductances of excitatory and inhibitory synaptic channels, ESynE and ESynI are the reversal potentials of these channels, and aD, aV, a3, b1, and b2 define the synaptic weights of, respectively, V0D, V0V, V3, In1, and In2 inputs to the centers in accordance with their interactions shown in Fig 2B. An increase in NMDA in the spinal cord preparation affects not only the level of excitation of neurons representing the centers, but also the level of excitation and recruitment of CINs. To take this into account in our simplified model, in which CINs are not explicitly included and the pathways they mediate are replaced by direct synaptic connections (see Fig 1B vs. 1A), we need to make synaptic weights of these connections to increase with the NMDA concentration represented by the parameter α. This especially concerns the excitatory commissural pathways that are amplified by neuron recruitment via excitatory interactions within each population and may also involve other excitatory neurons, such as the V2a, which are also recruited and activated with an increase in the NMDA concentration and locomotor speed and mediate inputs to the excitatory V0v CINs [37,57,58]. In other words, the increase of synaptic weights of CIN pathways was included in our simplified model to compensate for removal of CIN neurons whose excitability should also be affected by the drug. Therefore, we suggested that: a3=a3O⋅ (1+βa⋅α); (20) aV=aVO⋅ (1+βa⋅α), (21) where aVO and a3O define the basal weights of the corresponding excitatory connections and βa defines the rate of their increase with α. All model parameters are specified in Table 1. Synaptic weights were first set by manual adjustment to match experimental recordings. Weights were then optimized by running several iterations of one-dimensional bifurcation diagrams (see next section), where the robustness of a given regime was assessed by the width of the parameter range across which it displayed proper behavior, i. e. matching experimental phase diagrams and flexor-flexor phase relationships. Simulations were performed using custom C++ programs and visualized using gnuplot. To inspect changes in the model behavior with changing neuronal excitation (defined by the parameter α), phase differences between the left (LF) and right (RF) flexor activities or between ipsilateral flexor (LF) and extensor (LE) activities were calculated as the time differences between onsets of the corresponding bursts divided by the oscillation period measured in the activity of the flexor center. In the chosen range, α ∈ [0,1. 2] was progressively increased or decreased in a step-wise manner with a step size comprising 1/1000 of the range. At each point, the simulation was run and the phase differences were computed. Moreover, at every time step the initial conditions were chosen as a final state of the system from the previous time step to minimize a transient period. After convergence of the phase difference to a steady state with preset accuracy its value was plotted versus the parameter value. For each bifurcation diagram two series were generated: one with α increasing from 0 to 1. 2 in the described manner and another with α decreasing from 1. 2 to 0. The qualitative changes in system behavior (bifurcations) can be seen on such diagrams as jumps, corners, branching and other types of discontinuities. In addition, the existence of non-overlapping branches obtained by changing the parameter in different directions serves as an unequivocal indication of bi- or multi-stability in the system in the corresponding parameter ranges. In addition to the phase differences, the values of amplitude and instantaneous frequency (reciprocal to the oscillation period) were also calculated and the corresponding plots were built. Furthermore, the frequency of flexor activity was used to construct diagrams where phase differences between left (LF) and right (RF) flexors were plotted against the frequency. The experimental data included for comparison with modeling results had been obtained in earlier experimental studies published in Nature [37]. No new animal data were collected. All experiments were approved by the local ethical committee and performed in accordance with European guidelines for the care and use of laboratory animals. Briefly, the experiments were performed in the isolated spinal cords of wild-type (control) and transgenic mice with ablated V0D or V0V or both V0 neuron types (the detailed description of the transgenic lines of mice used can be found in [37]. Mice aged E18. 5 (with genetically deleted V0D neurons) or newborn mice aged 1–2 days (for all other studies) were used. The isolated spinal cords in a chamber were perfused with normal Ringer’s solution. Locomotor-like activity was induced by the exogenous application of mixtures of serotonin (or 5-hydroxytryptamine, 5-HT) and NMDA (N-methyl-D-aspartate). The locomotor frequency was controlled by concentrations of locomotor-inducing drugs (mostly NMDA) [37,40]. All recordings were performed at room temperature (22–24°C). Locomotor activity was recorded with suction electrodes attached to the L2 and L5 lumbar roots on both sides of the cord. The raw activity was band-pass filtered at 100 Hz to 1 kHz. Data points for analyzing cycle periods and phases were taken after the locomotor activity had stabilized 10–15 min after the initial burst of activity. All details of recordings and data processing can be found in [37].
The performance of the intact network is shown in Fig 3. Panel A of the figure shows the changes in the output activity, f (V), of all four centers with α changing from 0 to 1. 2. The vertical dashed lines in this panel indicate that activities of left (LF) and right (RF) flexor centers alternate at all values of α. Panel B shows how the frequency of oscillations (top diagram), the amplitude of flexor activity (second diagram), and the phase differences between the activities of left and right flexor centers (LF-RF) and left flexor and left extensor centers (F-E) (two bottom diagrams) changed with α. Note that in contrast to panel A, panel B shows α being changed in both directions, first forward from 0 to 1. 2 (red) and then backward from 1. 2 to 0 (blue). However, in all diagrams in Fig 3B, the graphs for forward and backward changes of α fully overlapped. Fig 3A and two top diagrams in Fig 3B show that with increasing α (simulating an increase of neuronal excitation in the experimental preparations by the administration of NMDA) the amplitude of flexor center activity in the model monotonically decreases and the locomotor frequency monotonically increases, which generally fits experimental data [40]. The bifurcation diagram depicting the phase difference between flexor centers in the intact model (Fig 3B, Phase diff. LF-RF) shows two intervals where the model exhibits qualitatively different behavior. On the right (α > 0. 55) flexor centers demonstrate perfect anti-phase synchronization (Δφ = 0. 5). This represents a symmetric alternation between the activities of left and right flexor centers that is in contrast to the behavior seen at lower values of α (α > 0. 55) where the phase difference between the alternating flexor center activities is not 0. 5 and a pitchfork bifurcation occurs at α ≈ 0. 55. The branches of the pitchfork reflect direct inhibitory interactions between the flexor centers (rebound) and differ depending on which left or right flexor center activates first. Moreover, the supercritical pitchfork bifurcation results in instability of the symmetric anti-phase regime with a decrease in α. The existence of two symmetric branches both emerging from Δφ = 0. 5, one with Δφ > 0. 5 and the other with Δφ < 0. 5, is concerned with the symmetry of the underlying network. It is worth mentioning that a slightly non-symmetric network would favor one of these two branches and weaken the other. In this case, the pitchfork bifurcation would be replaced with a saddle-node bifurcation resulting in the emergence of the “weak” branch, while the “strong” branch would replace the perfect anti-phase regime. Accordingly, the major qualitative change in the system dynamics as α is decreased below the bifurcation value of α ≈ 0. 55 is concerned with the appearance of two distinct stable regimes with phase shift between flexor center activities greater and less than 0. 5. It is important to note that at high values of α the periods of flexor and extensor activity have the same duration, while at low values of α the flexor phases can be significantly shorter than the extensor ones (see left parts of the graphs in Fig 3A). However, although the activity of left and right flexor centers in this case is not symmetric (does not alternate in perfect anti-phase), the left and right flexor centers are never active at the same time. Therefore, at any value of α (and at any locomotor frequency, see the top diagram in Fig 3B), the intact network exhibits left-right alternation of flexor activity. After both V0 commissural connections are eliminated, the interactions between left and right centers are provided exclusively by the excitatory V3 pathways (see Fig 1). As shown in Fig 4 (panel A and two top diagrams in panel B), similar to the intact model the amplitude of flexor activity monotonically decreases with increasing α, whereas the locomotor frequency monotonically increases. In contrast to the intact case, the activities of left and right centers are fully synchronized at all values of α. This is seen in Fig 4, panel A (indicated by vertical dashed lines) and panel B (the phase difference between the activities of left and right flexor centers remains 0 or 1 for all values of α - see “Phase dif. LF-RF” diagram). This left-right synchronous activity corresponds to hopping behavior. Selective deletion of the V0V subtype of V0 CINs corresponds, in our model, to elimination of V0V excitatory connections from extensor centers to the contralateral flexor centers (see Fig 1). In this case, the left and right sub-networks interact only through the excitatory (V3) and inhibitory (V0D) pathways between the flexor centers. Independent of commissural connections and similar to the intact case, the amplitude of flexor activity monotonically decreases with increasing α, whereas the locomotor frequency monotonically increases (see Fig 5, panel A and two top diagrams in panel B). Depending on the balance between the excitatory V3 and inhibitory V0D pathways, the net interactions can be excitatory (if the excitatory connections prevail) or inhibitory (if mutual inhibition is stronger than excitation). In our model, a pharmacologically-induced increase in excitation is implemented as an increase in the leak reversal potentials in the centers as well as an increase in the weights of excitatory pathways (to account for drug effects on neurons involved in these pathways). Based on our suggestions, the model parameters are set in such a way that at low levels of excitation (at low values of α) the inhibitory connections between flexor centers dominate, while at high excitation (at high values of α) the net interaction becomes excitatory. As a result, the model demonstrates left-right alternation of flexor activity at low values of α and left-right synchronization at high values of α (both indicated by vertical dashed lines in Fig 5A). The bifurcation diagram in Fig 5B (see “Phase dif. LF-RF” diagram) shows that at extremely low levels of excitation (α < X1) the activities of flexor centers alternate similarly to the intact case (compare the two branches left of X1 in “Phase dif. LF-RF” diagram in Fig 5B with the corresponding diagram in Fig 3B). At very high excitation levels (α > X5) left-right synchronous activity occurs similarly to the case when all V0 pathways are removed (compare two blue branches in the right part of “Phase dif. LF-RF” diagram in Fig 5B with the corresponding panel in Fig 4B). For intermediate excitation levels (X1 < α < X5) the system demonstrates a complex transition scenario from alternating to synchronous activity that is explained below using the technique of fast-slow decomposition. As seen in Fig 2A, the left knee on the V-nullcline is affected by changes in excitation to a much greater extent than the right knee (compare, for example, the width of the shaded area in terms of hNaP values at V = –55mV and at V = –30mV). Therefore, the activity of each flexor center can be affected by a synaptic input during a period of its inactivity, i. e. during the ipsilateral extensor phase, more efficiently than during its active phase. Due to intrinsic bursting properties each flexor center is capable of endogenous escape and inactivation. Accordingly, two types of strong synchronizing events are possible: a release of one flexor center on an escape of the other (release-on-escape) and a release of one flexor center when the other deactivates (release-on-shutdown). The former mechanism creates (in-phase) synchronization between flexor centers (“hopping”), while the latter underlies their alternation (anti-phase synchronization). Fig 6A–6C illustrates and allows for a more comprehensive understanding of system behavior after V0V pathways are deleted. Removal of V0D corresponds to termination of reciprocal inhibition between the flexor centers (see Fig 1). Since V0V-dependent connections remain intact, the network has two competing mechanisms of left-right coordination. One of them concerns reciprocal excitation between the flexor centers provided by V3 connections, which tends to synchronize their activity. The other is characterized by excitation from contralateral extensor centers to each flexor center, which contributes to alternation of flexor center activities. The results of simulations in this case are shown in Fig 7A and 7B. Again, similar to the intact case, the amplitude of flexor activity monotonically decreases with increasing α, whereas the locomotor frequency monotonically increases (see panel A and two top diagrams in panel B). In contrast to the previous case, when V0V pathways were removed, the model now demonstrates left-right synchronization of flexor center activity at low values of α and left-right alternation at high values of α (indicated by vertical dashed lines in Fig 7A). The transition scenario from synchronization to alternation is best seen in “Phase dif. LF-RF” diagram in Fig 7B. The above analysis of our model is centered on two major regimes characterized by left-right alternation and left-right synchronization (hopping). We have shown that the operating regime depends on the network integrity (presence or lack of particular commissural interactions), level of excitation in the network, and initial conditions. The general conclusion from this analysis is that is in the intact system, the left-right alternation is secured by the V0D CIN pathways at low excitabilities/frequencies and by V0V CIN pathways at high excitation/frequencies. The regimes of operation and patterns exhibited are summarized in Table 2 below: In the following section we provide a more general description of the network to elucidate the critical elements of the model that lead to the synchronization properties described above. In the model proposed the connections between neurons are relatively weak. In this case, the dynamics of the system may be described in terms of phase synchronization. This approximation implies that the trajectory of each oscillator in its own phase subspace is perturbed negligibly by the interactions between the oscillators, and hence the state of each oscillator can be described by a single variable, phase. Accordingly, taking into account the left-right symmetry of the system, the full set of differential equations describing the system can be reduced as follows in the approximation of asymptotically small connections (see, for example, [59]): ϕ˙1=Ω+F (ϕ1, ϕ2); ϕ˙2=Ω+F (ϕ2, ϕ1), (22) where ϕ1 and ϕ2 are the phases of the oscillators, Ω is an oscillation frequency and the function F (. , .) describes interactions between oscillators. For the phase difference ϕ = ϕ2 − ϕ1 we will have Hypothesizing that the right hand side of this equation depends only on the phase difference, i. e. F (ϕ2, ϕ1) − F (ϕ1, ϕ2) = G (ϕ), allows Eq (23) to be rewritten in a simpler form: The synchronized regimes of this system are the fixed (equilibrium) points of this first order differential equation (i. e. solutions of the equation G (ϕ) = 0) and their stability is defined by the sign of the first derivative dG (ϕ) / dϕ at the equilibrium. Specifically, the point is stable if the first derivative dG (ϕ) / dϕ is negative and the point is unstable if the derivative is positive. If we introduce a potential function as P (ϕ) = −∫G (ϕ) dϕ, then stable equilibrium points in Eq (24) will be local minima (valleys) of the function P (ϕ), and unstable points will be local maxima (hills) of P (ϕ). Because G (ϕ) is a 2π-periodic function that is odd, i. e. maintaining symmetry about the origin such that –G (φ) = G (−φ), its Fourier series expansion may be represented as: Restricting the expansion to the first two terms gives the following equation: This equation always has at least two fixed points at ϕ = 0 and ϕ = π representing synchronization and alternation regimes, respectively. The former is stable whenever A – 2B < 0, and the latter is stable if A + 2B > 0. Accordingly, for A < −2B only the synchronization regime ϕ = 0 is stable, for A > 2B only the alternation regime ϕ = π is stable, and between these values for −2B < A < 2B both regimes are stable. This inequality has solutions only if B > 0, so the second term in the expansion (Eq (26) ) with positive coefficient B is responsible for the existence of bistability in the system. Another qualitative conclusion from these speculations is that greater values of the parameter A correspond to “less stable” synchronization and “more stable” alternation in a manner referring to these behavior' s basins of attraction. Conversely, smaller values of A make synchronization “more stable” and alternation “less stable”. It is straightforward to assume that excitatory connections between oscillators contribute to a decrease in A, while inhibitory connections increase A. In the system under consideration, V3 and V0D commissural pathways mediate the direct excitatory and inhibitory connections between the flexor centers, respectively, and V0V pathways represent excitatory connections from the contralateral extensor to each flexor center, and hence can be considered as effectively inhibitory pathways between the flexor centers. Accordingly, the parameter A can be constructed from synaptic weights of these connections as A = −V3 + V0V +V0D, where V3, V0V, and V0D represent the strengths of the corresponding pathways. Assuming that all considered pathways include interneurons whose excitability is affected by NMDA, we can suggest that the strength of these pathways linearly increases with α: Accordingly, we can rewrite an expression for A as: In these terms the bifurcation scenarios described above allow for an elegant qualitative interpretation (see Fig 9A–9D). In the case when both subtypes of V0 pathways are deleted V0V0=V0D0=0, k0V = k0D = 0, and the parameter A=−V30−α⋅k3 is always negative and becomes even more negative with increasing α due to the negative slope. Accordingly, the synchronization regime remains stable for all values of α and the alternation regime is always unstable (see Fig 9D). Let’s now consider the case when only V0V pathways are deleted. In this case V0V0=0 and k0V = 0. Therefore, A=−V30+V0D0+α⋅ (−k3+k0D). If we assume that the connections mediated by inhibitory V0D pathways have a much stronger basal component (V0D0>V30), but weaker dependence on α than the excitatory V3 pathways (k0D < k3), then the baseline activity of the system (with small α) will be the opposite of the case when both V0 pathways are deleted since A can now be positive and large enough. With an increase in α the α-dependent term will eventually prevail and we will get the same situation as when all V0 pathways are deleted. In Fig 9B this scenario can be seen on the bifurcation diagram shown at the bottom on the (α, ϕ) -plane. We start, when moving from smaller to larger values of α, with the synchronization regime unstable and alternation regime stable. Then the fixed points ϕ = 0,2π become stable through the supercritical pitchfork bifurcation and new unstable fixed points appear which now separate the attraction basins of synchronization and alternation regimes (shown by dashed blue line). As α further increases these points move towards, and eventually merge to, ϕ = π making it unstable. On the surface described by the potential function, P (ϕ), this scenario manifests itself by the valley at ϕ = π (Δφ = 0. 5) where the surface becomes progressively more shallow and finally becomes a hill (see α close to 5). With V0D deleted, the expression (28) takes the form: A=−V30+VV0+α⋅ (−k3+k0V). Here we make an opposite assumption, i. e. V0V interactions have a smaller basal component than V3 pathways, i. e. V0V0<V30, but stronger α-dependent component, i. e. k0V > k3. Accordingly, the baseline in this case coincides with the case when both V0 pathways are deleted, i. e. we see stable synchronous and unstable alternation regimes at low values of α (see Fig 9C). As α increases the parameter A decreases due to its negative slope. This ultimately results in switching to the stable synchronization regime. Qualitatively this scenario is similar to the previous one but reversed (see the bifurcation diagrams on the (α, ϕ) -plane in Fig 9C). Finally, in the intact system both the intercept and the slope in Eq (28) are positive. The former is due to the fact that V0D0>V30, and the latter is because k0V > k3. Accordingly, A is always positive and only the alternating (out-of-phase synchronization) regime when ϕ = π (Δφ = 0. 5) is stable (see Fig 9A). In summary, this simplified model offers an explanation for why elimination of functionally similar interactions provided by V0V and V0D commissural pathways has such dramatically different effects. Based on our analysis, this may happen if the excitatory interactions between left and right rhythm generators are more dependent on (and change with) general neuronal excitation in the network than the inhibitory interactions. One of the most important characteristics of locomotion is the relative change in the flexor and extensor phase durations with changes in the locomotor speed or step cycle period [60,61]. Fig 10A shows how the durations of both phases in our intact model change with an increase in the locomotor period (slowing down locomotor oscillations). These data demonstrate clearly asymmetric changes in phase durations with variations of the step cycle, so that changes in the duration of the flexor phase are significantly less than changes in the duration of the extensor phase. This means that an increase in locomotor frequency mainly occurs due to shortening the extensor phase. Fig 10B shows a similar diagram representing the changes in flexor and extensor phase durations vs. step cycle period built using the recordings from wild-type mouse spinal cord preparations. One can see that the asymmetric changes in phase durations in experimental studies are qualitatively similar to those in our simulations.
The main objective of this theoretical study was to formulate and analyze a simplified model of the spinal circuitry in order to reproduce and explain recent experimental data concerning frequency-dependent differential contributions of two distinct commissural pathways to the left-right alternation of neural activity in the isolated spinal cord [37]. We used a simplified connectivity scheme adapted from both the experimental data [37] and the previous large-scale computer models [41]. Bifurcation diagrams were constructed for each experimental case to observe the system’s behavior over a large range of excitability (Figs 3–5 and 7). These diagrams allowed us to determine the values of parameters leading to the observed behaviors and provide a strong mathematical explanation of the critical changes seen in left-right phase relationships leading to a normal alternating pattern and to switching to synchronous (hopping) behavior. Moreover, the bifurcation diagrams revealed regions of bistability, where two stable states co-existed for a given set of parameters. To complete our understanding of the dynamics underpinning the transitions, we used fast-slow decomposition methods and generated phase plane diagrams that corresponded to removal of V0V or V0D commissural pathways (Figs 6 and 8). Our analysis showed that left-right alternation and synchronization occurred in both cases when one of the pathways (V0 or V0D) was removed, because of the release-on-shutdown and release-on-escape mechanisms, respectively. These findings allowed the interpretation of experimental results in terms of qualitative theory of dynamical systems. In addition, at the intermediate locomotor frequencies, our model exhibited bistable behaviors and predicted a coexistence of both alternating and synchronized regimes, leading to an intermittent alternating/synchronous activity observed in experimental studies (Fig 15). The exact intrinsic cellular mechanisms involved in rhythmic bursting in the spinal cord remain unknown. Previous modeling studies suggested that these mechanisms may involve the persistent (or slowly inactivating) sodium current, INaP, [9,41–46,51,52] and this suggestion has been supported by several experimental studies [52–56]. Although we used the INaP-dependent mechanism in this model, the results of this study concerning network behavior can be considered independent of the exact cellular mechanism employed. At the same time, it should be noted that some particular features of the model fit to known experimental data specifically due to the INaP—dependent mechanism incorporated (although they can be provided by other cellular/network mechanisms). These features include (a) a monotonic increase of burst frequency and (b) monotonic decrease of (flexor) burst amplitude with an increase of neuronal excitation (by increasing NMDA or glutamate concentration, see [37,40] and Figs 3–5 and 7), as well as (c) small changes in the flexor phase duration relative to the extensor phase duration at low frequencies (the same figures). These observations provide additional indirect support for the idea that INaP—dependent mechanisms are involved in the generation of bursting activity in the spinal cord. A recent study using optogenetics [62] has demonstrated that locomotor-like rhythmic bursting can be induced independently in flexor and extensor networks. This means that both flexor and extensor populations can potentially, under certain conditions, generate endogenous locomotor-like, oscillatory activity. However these findings do not necessary imply that both these populations are in the bursting mode under the conditions of isolated spinal cord preparation with locomotor oscillations induced by drugs or other methods. Moreover, analysis of non-resetting deletions (missing bursts) in this preparation showed that missing flexor bursts were always accompanied by a sustained (tonic) extensor activity, whereas missing extensor bursts occurred without an effect on flexor bursting [45]. Finally, the duration of extensor activity between flexor bursts at low burst frequency can be pretty long, exceeding 5–6 seconds when flexor burst duration is usually less than 2. 5 seconds [37], which also implies that intrinsic bursting is unlikely to occur in rodent extensors. Based on the observations described above, we suggest that each flexor or extensor center can potentially, in certain conditions, intrinsically generate bursting activity. However, we think that under the conditions considered here (i. e. , the isolated spinal cord with 5-HT and NMDA-induced locomotor-like oscillations) only the flexor centers (left and right) operate in a bursting mode whereas the extensor centers (left and right) operate in a state of tonic activity (if isolated) and exhibit rhythmic bursting because of rhythmic inhibitory inputs from the corresponding ipsilateral flexor center (see the section “Rhythm-generating Centers and Connections”). This may occur because the mixture of 5-HT and NMDA used to evoke locomotor activity in the isolated spinal cord, excites the extensor centers much stronger than the flexor centers (and hence ELFO < ELEO, see Fig 2B). The above suggestion leads to the concept that each (left and right) CPG is asymmetric, so that the flexor center on each side differs from the extensor center not only by a level of excitation described above, but also by the organization of left-right and flexor-extensor interactions (see Fig 1B and weights of interactions in Table 1). As a result of this flexor-extensor asymmetry, the burst frequency and the flexor phase duration are mostly defined by the flexor centers, and hence at low frequencies the duration of the flexor phase may be significantly less than the duration of the extensor phase. Fig 10A shows that when frequency goes down (with decreasing α) the increase of step cycle period mainly results from an increase of the extensor phase duration with a very small increase of the flexor phase duration. Interestingly, this phase durations vs. period plot built using simulation data fits well to the plot built using data from experimental studies [37] (Fig 10B). This provides additional support for the idea of an asymmetric CPG employed in our study. The concept of an asymmetric CPG organization with a dominant role of flexor centers in rhythm generation has been previously suggested by Duysens and Pearson [63–66] and implemented in recent computational models [41,45]. Such asymmetric flexor-extensor phase relationships are consistent with the previous experimental studies in cats [60] and support the idea that this flexor-extensor asymmetry is an inherent property of the CPG [61] rather than a result of afferent influences on a quasi-symmetric CPG [67]. Our simulations support the earlier suggestion that left-right alternation is provided by two commissural pathways involving the inhibitory V0D and excitatory V0V CINs [37]. Similar to the experimental [37] and previous modelling [41] studies, our simulations demonstrate that the contribution to left-right alternation of the V0D pathway is dominant at low frequencies and reduces as locomotor frequency increases, whereas the contribution of the excitatory V0V pathway is weak at low frequencies but becomes stronger as frequency increases. The analysis of our models suggests that the different frequency-dependent roles of these pathways are based on the following: (1) a relatively weak dependence of V0D CIN pathways on the neuronal excitability (neuroactive drug concentration), which leads to a net reduction of V0D activity due to the reduction of amplitude of rhythm-generating activity when the oscillation frequency increases; (2) a strong dependence of V0V CIN pathways on the neuronal excitability, leading to the increase of their activity with oscillation frequency. This may be based on the existing data that inputs to V0V CINs are often mediated by ipsilaterally-projecting excitatory V2a neurons [37,68,69] These were shown to increase their activity and recruitment with an increase in neuroactive drug concentration [58,70] and their selective ablation mimics the V0V ablation [70]. Our study support the concept that the resultant locomotor pattern and locomotor gait depend on the balance between different commissural interactions, which, in turn, depends on the level of neuronal excitation and locomotor speed dictated by the locomotor task. This conclusion is consistent with a general concept that the expression of a particular left-right coordinated pattern, such as the left-right alternation of activity typical for regular walking or the left-right synchronized activity characteristic for hopping or galloping, depends on the balance between the commissural pathways working for or against synchronization of left and right spinal CPGs. Therefore, under normal conditions this pattern and the corresponding gait can be changed by additional speed-dependent or speed-independent descending or afferent signals to different commissural interneuron populations involved in these pathways.
Movements of left and right limbs in mammals during locomotion are controlled by distinct rhythm-generating neuronal circuits in the spinal cord. Complex interactions between these circuits provide flexible coordination of limb movements in different gaits. It was shown that interactions between left and right spinal circuits are mediated by commissural interneurons. Genetic ablation of a particular type of these interneurons, called V0, leads to switching from a regular, left-right alternating “walking” activity to a left-right synchronous “hopping” pattern. Moreover, the V0 commissural interneurons have excitatory and inhibitory subtypes that appear to play different roles in the left-right coordination depending on locomotor speed. In this theoretical study, we build a simplified mathematical model of spinal circuits that describes left and right rhythm generators interacting bilaterally via several types of commissural connections. Using this model, we simulate different experimental manipulations, analyze the resultant alternating and synchronous regimes of activity, and propose explanations for the results of experimental studies. We show that although both excitatory and inhibitory V0 commissural pathways support left-right alternation, the resultant locomotor pattern and gait depend on the balance between different commissural interactions, which in turn may depend on the level of neuronal excitation and locomotor speed.
Mechanisms of Left-Right Coordination in Mammalian Locomotor Pattern Generation Circuits: A Mathematical Modeling View