Phase-Dependent Response to Afferent Stimulation During Fictive Locomotion: A Computational Modeling Study

Abstract

Central pattern generators (CPGs) in the spinal cord generate rhythmic neural activity and control locomotion in vertebrates. These CPGs operate under the control of sensory feedback that affects the generated locomotor pattern and adapt it to the animal’s biomechanics and environment. Studies of the effects of afferent stimulation on fictive locomotion in immobilized cats have shown that brief stimulation of peripheral nerves can reset the ongoing locomotor rhythm. Depending on the phase of stimulation and the stimulated nerve, the applied stimulation can either shorten or prolong the current locomotor phase and the locomotor cycle. Here, we used a mathematical model of a half-center CPG to investigate the phase-dependent effects of brief stimulation applied to CPG on the CPG-generated locomotor oscillations. The CPG in the model consisted of two half-centers mutually inhibiting each other. The rhythmic activity in each half-center was based on a slowly inactivating, persistent sodium current. Brief stimulation was applied to CPG half-centers in different phases of the locomotor cycle to produce phase-dependent changes in CPG activity. The model reproduced several results from experiments on the effect of afferent stimulation of fictive locomotion in cats. The mechanisms of locomotor rhythm resetting under different conditions were analyzed using dynamic systems theory methods.

Introduction

The mammalian spinal cord contains neuronal circuitry that can generate a basic locomotor rhythm and produce the alternating flexor and extensor motoneuron activities underlying locomotion. Although this locomotor central pattern generator (CPG) can operate in the absence of sensory feedback (reviewed by Grillner, 1981; Rossignol, 1996; Orlovsky et al., 1999; Rossignol et al., 2006), afferent feedback plays a crucial role in adjusting the locomotor pattern to the motor task, the environment, and the biomechanical characteristics of the limbs and body (e.g., Pearson, 2004; Rossignol et al., 2006). Continuous electrical stimulation of the midbrain locomotor region in an immobilized decerebrate cat produces “fictive locomotion” consisting of rhythmic alternating activation of flexor and extensor motoneurons similar to that occurring during normal locomotion in an intact animal (see Rossignol, 1996). To investigate the effects of afferent inputs on the locomotor pattern generated by the CPG and step cycle timing, researchers often use the fictive locomotor preparations while applying stimulation to flexor or extensor sensory afferents (e.g., Guertin et al., 1995; Perreault et al., 1995; McCrea, 2001; Stecina et al., 2005). These studies revealed that in many cases, afferent stimulation can delay or accelerate the phase transition within the ongoing step cycle with or without changing the timing of the subsequent step cycles (Rybak et al., 2006b; McCrea and Rybak, 2007).

Although the anatomical structure of the CPG circuit remains unclear, the use of relatively simple mathematical models of CPGs allows the study of the general effects of afferent stimulation on CPG operation from a dynamic viewpoint. In particular, half-center type CPG models were previously used to reproduce some effects of sensory afferent stimulation on fictive locomotor pattern in cats (Rybak et al., 2006b).

The goal of the present study was to further investigate the mechanism for the phase-dependent response of the locomotor pattern during fictive locomotion using a simplified half-center CPG model. Specifically, we applied stimulation to the CPG model in different phases of the locomotor cycle and examined how the temporal activity of the CPG changed. The use of a relatively simple CPG model allowed us to apply the dynamic system methods and perform mathematical analysis to fully characterize the phase-dependent responses of the CPG to applied stimulation.

Methods

Model

It has been suggested that the rhythmic pattern of the CPG activity is determined in the rhythm generator (RG) network of the CPG (Rybak et al., 2006a,b). In the present study, the model (Figure 1) consisted of two neuron populations representing RG centers (flexor RG-F and extensor RG-E) and two populations of inhibitory interneurons (In-F, In-E), providing mutual inhibition between the flexor and extensor centers. Each population was described as an activity-based (non-spiking) neuron model (Ermentrout, 1994; Markin et al., 2010; Molkov et al., 2015; Danner et al., 2016, 2017). The state of each neuron was characterized by the membrane potential Vi (i = F, E, IF, IE), where the indexes F and E are used for the RG-F and RG-E neurons, respectively, and the indexes IF and IE are used for the In-F and In-E neurons, respectively. The RG-F and RG-E neurons incorporated a persistent (slowly inactivating) sodium current that defined intrinsic rhythmogenic properties of these neurons. The intrinsic oscillation in each RG neuron depended on the variable hi (i = F, E) that defined slow inactivation of the persistent sodium channels. Each RG center could produce rhythmic activities; however, if uncoupled, the extensor center was in the tonic regime due to a supraspinal drive and produced sustained activity. Rhythmic oscillations of the RG were defined by the flexor centers, which provided rhythmic inhibition of the extensor center through In-F. The supraspinal drive to the flexor center determined the oscillation frequency. Synaptic interactions between all neurons in the model are shown in Figure 1. For the state variable of this model, we usedV=[VF,VE,VIF,VIE]Tandh=[hF,hE ]T.

FIGURE 1

Figure 1. Model schematic of the rhythm generator (RG) network and afferent inputs. The RG network is composed of flexor (RG-F) and extensor (RG-E) centers inhibiting each other via inhibitory interneurons In-F and In-E, respectively. The supraspinal drive provides excitation to the RG-F and RG-E neurons defining the frequency of oscillations. Sensory afferents can synaptically excite both RG neurons and inhibitory interneurons.

The dynamics of the membrane potential Vi of the RG neurons (i = F, E) and the interneurons (i = IF, IE) is described as

CV˙i={−INaP(Vi,hi)−ILeak(Vi)−ISynEi( V)−ISynIi( V) i=E,F −ILeak(Vi)−ISynEi(V)−ISynIi( V) i=IF ,IE (1)

where C is the membrane capacitance, INaP is the persistent sodium current, ILeak is the leak current, andISynEiandISynIi are the currents by excitatory synapses and inhibitory synapses, respectively. The ionic current INaP and leak current ILeak are described as

INap(Vi,hi)=g^NapmNap(Vi)hi{Vi−ENa} i=F,E ILeak(Vi)={ g^LeakRG{Vi−ELeakRG} i=F,E g^LeakInRG{Vi−ELeakInRG} i=IF,IE (2)

whereg^Nap,g^LeakRG, andg^LeakInRG are the maximum conductances of the corresponding current, and ENa,ELeakRG, andELeakInRG are the reversal potentials of the corresponding current. In addition, mNap is the activation of the sodium channel of the RG neurons and is described as

mNap(Vi)=11+exp(-Vi+40.06.0) i=F,E (3)

The dynamics of the inactivation of the sodium channel hi of the RG neurons (i = F, E) is given by

τ(Vi)h˙i=h∞(Vi)-hi i=F,E (4)

where

h∞(Vi)=11+exp(Vi+45.04.0) τ(Vi)=320+320cosh(Vi+35.015.0) ms i=F,E (5)

The currents generated by the synapsesISynEiandISynIiare given by

ISynEi(V)=g^SynE{Vi-ESynE}{∑j={F,E,IF,IE}aijf(Vj)+cid+wisi}ISynIi(V)=g^SynI{Vi-ESynI}{∑j={F,E,IF,IE}bijf(Vj)} i=F,E, IE,IF (6)

whereg^SynEandg^SynI are the maximum conductances of the corresponding current, ESynE and ESynI are the reversal potentials of the corresponding current, d is the tonic drive from the supraspinal region, si (i = F, E, IF, IE) is the feedback input from sensory fibers, and aij, bij, ci, and wi (i, j = F, E, IF, IE) are the weight coefficients. Moreover, the output function f translates V into the integrated population activity and is given by

f(Vi)={0Vi<VthVi−VthVmax >Vi≥Vth1Vi≥Vmax i=F,E,IF,IE (7)

where Vth and Vmax are the lower and upper threshold potentials, respectively. The differential equations of Equations (1), (2), and (4) were solved numerically using the fourth-order Runge-Kutta method with a step size of 0.01 ms. The parameter values are shown in Appendix A.

Modeling the Effects of Phase-Dependent Afferent Stimulation

The CPG model produced rhythmic activity and exhibited stable oscillations, as shown in Figure 2. The active phase for each neuron was defined as the time interval during which the neuron's potential was higher than Vth and the silent phase as the time interval when the potential was lower than Vth. The cycle period T was defined as the time interval between two consecutive onsets of the active phase. The phase of oscillation was defined as ϕ = 2πt/T ∈ [0, 2π).

FIGURE 2

Figure 2. Changes of membrane potential of RG-F (top panel) and RG-E (bottom panel) neurons (A) without any stimulation, (B) with stimulation applied to the flexor side during flexor phase, and (C) with stimulation applied to the flexor side during extensor phase. The red bars indicate the application of stimulation. Gray regions indicate active phases. The applied stimulation increased the duration of the current flexor phase and cycle period in (B) and initiated the flexor phase and decreased the cycle period in (C). Both stimulations produced phase shifts.

The CPG also received external (“sensory”) signals (Figure 1). Based on a previous study (Demir et al., 1997), which investigated the response of a single neuron model to stimulations, we used depolarizing stimuli applied at different phases of oscillatory activity. Specifically, after oscillation stabilized, we applied a 200 ms stimulus to the flexor (RG-F and In-F) or extensor (RG-E and In-E) neurons. The intensity of stimulation was sF = sIF = 0.2 and sE = sIE = 0.0 for the flexor side and sF = sIF = 0.0 and sE = sIE = 0.2 for the extensor side in Equation (6). Suppose that the neuron activity is perturbed by the stimulation at phase ϕs ∈ [0, 2π) and the period changes from T to α(ϕs), as shown in Figure 2. To show the phase shift of the neuron activity in response to the stimulation, we define

Δ(ϕs)=2πα(ϕs)-TT (8)

Calculation of Nullcline

The nullcline is a set of points at which the derivative of a differential equation is equal to zero. It reflects the structure of the solution of the differential equation. To investigate the mechanism of the phase-dependent response of the CPG model, we used a nullcline-based method. The state variable of the CPG model is given by (V, h). The nullclines for the RG neurons are given by

NiV={(V,h) | V˙i=0}Nih={(V,h) | h˙i=0}i= F, E (9)

To clarify the dynamics of each RG neuron, we focused on the Vi-hi space (i = F, E) for the nullclines by assuming that the other variables Vj (j = F, E, IF, IE, j≠i) and hk (k = F, E k≠i) are on the stable oscillation with phase ϕ. Therefore, we modifyNiVandNihin Equation (9) as

N^iV(ϕ)={(Vi,hi) | V˙i=0, Vj=Vj(ϕ),hk=hk(ϕ)}N^ih(ϕ)={(Vi,hi) | h˙i=0, Vj=Vj(ϕ),hk=hk(ϕ)} i=F,E j=F,E,IF,IE j≠i k=F, E k≠i (10)

ForN^iV(ϕ)andN^ih(ϕ) , we can write hi = hi(Vi; ϕ), as explained in Appendix B.

Results

Phase-Dependent Response

Figure 3A shows the phase shift Δ of the RG-F neuron activity after stimulation of sensory inputs on the flexor side at ϕs. When stimulation was applied during the silent phase of RG-F (2.51 ≤ ϕs < 2π), it caused the transition to the active phase to occur earlier and this advanced start decreased with ϕs. In contrast, almost no phase shift occurred when stimulation was applied at the beginning of the active phase of RG-F (0 ≤ ϕs < 1.00). However, the neuron activity was delayed by the stimulation during the middle and end of the active phase (1.00 ≤ ϕs < 2.51). These trends were similar to those observed during fictive locomotion in cats (Schomburg et al., 1998; Frigon et al., 2010), as shown in Figure 3B. Figure 3C shows Δ of the RG-F neuron activity by the stimulating sensory fibers of the extensor side. The active and silent phase of the RG-F neuron corresponds to the silent and active phase, respectively, of the RG-E neuron. The neuron activity was advanced at the middle of the silent phase of the RG-E neuron and was delayed at the end of the active phase of the RG-E neuron. The response of the stimulation of the extensor side was qualitatively similar to that of the flexor side. Moreover, these trends were similar to those seen in animal experiments (Schomburg et al., 1998; Frigon et al., 2010). The effects of the stimulation duration and intensity are further investigated in Figure S1 in Appendix C.

FIGURE 3

Figure 3. (A) Phase-dependent response of the RG-F neuron by stimulating sensory fibers of the flexor side. (B) Response against flexor muscle stimulation during fictive locomotion in cats (adapted from Schomburg et al., 1998). The flexion (extension) phase corresponds to the active (silent) phase of the RG-F neuron of the CPG model. (C) Phase-dependent response of the RG-F neuron by stimulating sensory inputs on the extensor side.

Analysis on Nullclines

Even though the oscillatory behavior of the RG-E neuron was similar to that of the RG-F neuron as shown in Figure 2, the oscillating mechanism was different due to different nullclines as suggested in previous studies (Spardy et al., 2011a,b; Molkov et al., 2015). To understand this mechanism, we briefly explain the roles of nullclines in our neuron model. Figure 4A shows the nullclinesN^FV,N^EV,N^Fh, andN^Ehwith the vector field for the case without synaptic connections from other neurons. WhileN^FhandN^Ehare identical and have a sigmoid shape,N^FVandN^EVhave different cubic curves. In particular, whileN^FVhas two distinct inflection points and the sign of the slope changes at the inflection points,N^EVchanges monotonically. Because two eigenvalues at the intersection ofN^FVandN^Fhare positive and negative, the intersection is a saddle, which induces a limit cycle (orange orbit) due to the following three characteristics; (1) the trajectory approachesN^FV, especially its branches with positive slope due to the difference of the time constants between the dynamics of V and h, (2) the trajectory close to the positive branches moves along them until reaching the inflection points, and (3) the trajectory jumps to the opposite positive branch at the inflection points. In contrast to the case for the RG-F neuron, the two eigenvalues at the intersection ofN^EVandN^Ehare both negative and the intersection is stable node. The trajectory is attracted to this node and stays there as long as the node exists. Therefore, the RG-E neuron does not show any oscillatory behavior.

FIGURE 4

Figure 4. Roles of nullclines of RG-F and RG-E neurons to produce oscillatory behaviors. The green lines showN^FhandN^Eh. The red and blue lines showN^FVandN^EV, respectively. Circles indicate intersections of nullclines [filled circles for both negative eigenvalues (stable node) and open circles for negative and positive eigenvalues (saddle)]. (A)N^FVandN^EVwith the vector field for the case without synaptic connections from other neurons. The saddle produces a limit cycle (orange orbit) while stable node does not produce any oscillatory behavior. (B) Schematic illustration of changes inN^FVandN^EVinduced by synaptic connections from other neurons. The intersection ofN^FVandN^Fhalmost remains saddle (burst mode), which induces an oscillatory behavior. On the other hand,N^EV transitions between two positions depending on the inhibitory signal from the contralateral side. This transition produces oscillatory behavior between tonic and silence modes for the extensor side. (C) Detailed illustration of our model at ϕ = 0, 0.89, 1.78, 2.68, 3.88, and 5.08 rad. Red and blue diamonds are (VF, hF) and (VE, hE), respectively, and these points move in accordance with eigenvalues, as indicated by arrows.

The synaptic connections from other neurons changeN^FVandN^EV as shown schematically in Figure 4B, so that both RG-F and RG-E neurons show oscillatory behavior. On the one hand, although the intersection ofN^FVandN^Fhtemporarily forms a stable node, it remains close to the saddle point (burst mode), which produces an oscillatory behavior. On the other hand, while the intersection ofN^EVandN^Ehremains stable,N^EV transitions between two positions due to an inhibitory signal from the contralateral side, one of which has a high V at the intersection (tonic mode) and the other of which has a low V (silence mode). These transitions produce an oscillatory behavior. Figure 4C shows the details of our model at ϕ = 0, 0.89, 1.78, 2.68, 3.88, and 5.08 rad to show how the nullclines changed during one cycle.

Shortening of Activity Duration During Silent Phase

Next, we investigated the mechanism for the phase-dependent response during the silent phase. Figure 5 shows the responses on the VF-hF plane by the stimulation of the flexor side at ϕs = 3.77, 5.03, and 5.53 rad. The disturbed trajectories took a shortcut to the limit cycle at different positions depending on ϕs, which decreased the activity duration and advanced the neural activity. As shown in Equations (1), (4), and (6), while stimulation directly influences the membrane potential Vi (i = F,E), it does not influence the inactivation of the sodium channel hi (i = F,E). Therefore, a shortcut was produced in the direction of Vi. Moreover, for the same reason, as ϕs occurs earlier, the shortcut has a larger truncated trajectory and the neural activity is more advanced. Although the intersection ofN^FhandN^FVbefore the stimulation was in silence or burst mode, it suddenly changed to tonic mode after the stimulation, which attracted the trajectory toward the intersection and shortened the neuron activity.

FIGURE 5

Figure 5. Response of a RG-F neuron on the VF-hF plane by stimulating the flexor side at ϕs = 3.77, 5.03, and 5.53 rad. The black line shows the limit cycle without stimulation. Stimulation is applied at filled orange circles. Disturbed trajectories (orange line) take a shortcut to enter the limit cycle at different positions depending on ϕs. Earlier ϕs has a larger truncated trajectory. The green line showsN^Fh. The red dashed and solid lines showN^FVjust before and after the stimulation, respectively. The intersection ofN^FhandN^FVchanged from silence or burst mode to tonic mode by the stimulation.

Figure 6A shows the response on the VE-hE plane by the stimulation of the extensor side at ϕs = 1.13 rad. WhileN^EVmoved to the right and the intersection ofN^EhandN^EVchanged from the silence to tonic mode just after the stimulation, the movement ofN^EVwas smaller than that ofN^FV when the flexor side was stimulated (Figure 5). After the stimulation, although the disturbed trajectory moved to the right, it did not completely enter the limit cycle (① in Figure 6A). However,N^EVgradually moved to the right and the intersection ofN^EhandN^EV also further moved to the right. As a result, the trajectory eventually took a shortcut to the limit cycle (② in Figure 6A). Although the shortcut was induced by the change of the intersection ofN^EhandN^EV from the silence to tonic mode in the same way as that of the stimulation of the flexor side (Figure 5), it was delayed due to an inhibitory signal from the flexor side just after the stimulation. More specifically, Figure 6B shows the time profiles of the neurons after the onset of the stimulation. Just after the stimulation, the membrane potentials (VE and VIE) of the RG-E and In-E neurons increased immediately and crossed over Vth (① in Figure 6B), which changed the effect on the connected neurons described by Equation (7). The immediate change of the In-E neuron changed the activities of the other neurons. Especially, the membrane potentials (VF and VIF) of the RG-F and In-F neurons decreased due to the inhibitory signal from the In-E neuron and crossed over Vth. The decrease of the inhibitory signal from the flexor side increased VE (② in Figure 6B), which induced the shortcut.

FIGURE 6

Figure 6. (A) Response of RG-E neuron on the VE-hE plane by stimulating the extensor side at ϕs = 1.13 rad. The black line shows the limit cycle without stimulation. Stimulation was applied at filled cyan circle. The blue dashed and solid lines showN^EVjust before and after the stimulation, respectively. WhileN^EVmoved to the right just after the stimulation and the intersection withN^Eh(green line) became tonic mode, the disturbed trajectory (cyan line) moved to the right without entering the limit cycle (①).N^EVgradually further moved to the right (blue double line showsN^EV at 80 ms after the stimulation) and the trajectory was finally cut short to the limit cycle (②). (B) Time profiles of four neurons from the onset of the stimulation to the end of the shortcut. The vertical lines show the onset and 80 ms after the stimulation. The horizontal line shows Vth. After the stimulation, the membrane potentials of the RG-E and In-E neurons rapidly changed and crossed over Vth (①). After that, while the membrane potentials of the RG-F and In-F neurons decreased due to the inhibitory signal from the In-E neuron and crossed over Vth, the membrane potential of the RG-E neuron gradually increased. As a result, the decrease of the inhibitory signal from the flexor side increased the activity of the RG-E, which induced the shortcut (②).

Prolongation of Activity Duration During Active Phase

At the end of the active phase, the neural activity was delayed as shown in Figure 3. In the case without stimulation (VF, hF), of the RG-F neuron swooped down to the right inflection point ofN^FV at the end of the active phase, as shown in the panel for ϕs = 1.78 rad of Figure 4. However, the stimulation at the end of the active phase movedN^FVto the right and changed the intersection ofN^FhandN^FV from burst to tonic mode, as shown in Figure 7A. Furthermore,N^FVshowed almost no change for a while. These inhibited the deactivation of the RG-F neuron and prolonged the activity duration. In addition, the intersection ofN^EVandN^Eh changed from the burst to the silence mode and stayed in the silence mode for a while, which also delayed the neural activity. Figure 7B shows the case of the stimulation of the extensor side at the end of the active phase of the RG-E neuron. The RG-E neuron maintained the tonic mode due to the stimulation and this prolonged the activity duration. This response was similar to the case of flexor stimulation (Figure 7A).

FIGURE 7

Figure 7. Change ofN^FVandN^EV by stimulation at the end of the active phase of (A) the flexor side (ϕs = 2.39 rad) and (B) the extensor side (ϕs = 5.91 rad). Solid and dotted lines show the results of the cases with and without the stimulation, respectively. Red and blue diamonds mark the positions of (VF, hF) and (VE, hE), respectively (filled diamonds for stimulation and open diamonds for non-stimulation). In (A), the stimulation movedN^FVto the right and changed the intersection ofN^FhandN^FVfrom burst to tonic mode. Furthermore,N^FVshowed almost no change for a while. These prolonged the activity duration. In addition, they kept the intersection ofN^EVandN^Ehsilence mode for a while, which also delayed the neural activity. In (B), the stimulation movedN^EVto the right and maintained the intersection ofN^FhandN^FVin tonic mode. These prolonged the activity duration of the RG-E neuron.

Discussion

In the present study, we investigated the underlying mechanism of the phase-dependent response of a half-center CPG model by applying a brief stimulation to it. The simulation results showed trends in the phase-dependent responses similar to those observed during fictive locomotion in cats (Schomburg et al., 1998; Frigon et al., 2010; Figures 3A,B).

It has been reported that the locomotor rhythm is reset to start a new flexion phase by an electrical stimulation to the flexor nerve in animals (Schomburg et al., 1998). Our simulation results suggest that, while the locomotor rhythm is reset to start a new flexion phase by stimulation during the silent phase, its start phase depends on the stimulation phase. The phase shifts of the RG-F neuron during the active phase (silent phase of the RG-E neuron) were also induced by stimulation of the extensor side (Figure 3C). However, in contrast to stimulation of the flexor side, the change in the intersection of the nullclines was smaller and formation of trajectory shortcut did not occur just after the stimulation of the extensor side (Figure 6). Instead, the In-E neuron was activated by the stimulation (we can estimate this using Equation S7 in Appendix D), which deactivated the RG-F and In-F neurons due to the inhibitory signal from the In-E neuron. As a result, the RG-E neuron was activated because of the deactivation of the neurons in the flexor side. These processes delayed the shortcut after the stimulation of the extensor side. Although the shortcut was delayed by the stimulation of the extensor side, the RG-E neuron had the potential to produce an immediate shortcut by stimulation, as in Figure 5, due to the nullcline intersection changing to a tonic mode when the stimulation intensity was larger as illustrated in Figure S2 in Appendix C.

At the end of the active phase, the neural activity was delayed by the stimulation. When the flexor side was stimulated, the intersection of the nullclines of the RG-F neuron changed from burst to tonic mode (Figure 7A). Similarly, the stimulation of the extensor side at the end of the active phase of the RG-E neuron prolonged the active phase by maintaining the tonic mode (Figure 7B). Even though the parameters of synaptic connection were different between the flexor and extensor sides, the mechanism of the active phase prolongation was the same (Figures 7A,B). As Figure S2 in Appendix C shows, the stimulation contributed to the nullcline intersection changing to a tonic mode irrespective of ϕS. From our simulation results, the phase-dependency was caused by these acceleration and prolongation mechanisms, which were commonly induced by the change of the nullcline intersection to a tonic mode.

Contribution of Different Afferent Types

Schomburg et al. (1998) demonstrated the resetting of the locomotor cycle in response to various flexor nerve stimulation during fictive locomotion. They employed both shorter stimulation trains (around 60 ms) at stimulation intensities activating joint and cutaneous afferents and longer stimulation trains (over 200 ms) at intensities activating only group I and II afferents. Other studies investigating the effects of sensory afferents on locomotor modulation also used relatively longer stimulation (for example, Ia and II afferents of extensor and flexor were stimulated for over 125 ms in Frigon et al., 2010; Ia or Ib afferents of extensor were stimulated for over 500 ms in Whelan et al., 1995; and II afferents of flexor were stimulated for over 200 ms in Perreault et al., 1995). Based on the conditions of these experiments, we used a stimulation lasting 200 ms. In addition, the effect of the stimulation intensity was also investigated in those experiments. Therefore, we examined the effect of the stimulation duration and intensity (Figures S1, S2 in Appendix C).

Functional roles of muscle spindles (Ia and II), Golgi tendon organs (Ib), and cutaneous afferent inputs during locomotion have been investigated in previous studies. During the stance phase, feedback from muscle spindles and Golgi tendon organs of extensor muscles prolong the duration of extensor activity (Guertin et al., 1995; Whelan et al., 1995) and muscle spindles in hip flexors contributed to initiation of the swing phase (Hiebert et al., 1996). At the beginning of the swing phase, stimulation of cutaneous nerves prolonged this phase (Duysens, 1977). As indicated above, the different responses depended on the locomotor phase. Yet, it remains unclear how the neural circuit of the CPG interacts with different types if sensory fibers and which neural circuits contributed to the generation of a phase-dependent response. In our present model, we did not identify the relative contributions of different afferent types to the CPG (Figure 1). Nevertheless, our model reproduced a phase-dependent response (Figure 3). Further experimental and computational studies are necessary to delineate anatomically and functionally plausible interactions between the CPG and the sensory afferents.

Functional Roles of the Different Layers in CPGs

Although the anatomical structure of the CPG remains unclear, it has been suggested from modeling studies (Rybak et al., 2006a,b) that the CPG consists of a RG layer and a pattern formation (PF) layer. The PF layer is thought to determine the spatial motor pattern depending on the phase generated in the RG neurons; that is, it determines the distribution of the co-activated α-motoneurons over time. The muscle synergy hypothesis is one candidate for the determination of the distribution (Ivanenko et al., 2004, 2006) and modeling studies have shown that a motor control system based on this hypothesis could generate locomotion using musculoskeletal models (Aoi et al., 2010, 2013, 2019; Fujiki et al., 2018). In those models, the amplitudes of the α-motoneuron activities were determined in the PF layer. Based on this, it is suggested that the neurons in the PF layer modulate their amplitudes, which would be related to the phase-dependent response in terms of amplitude of the electromyography of Hoffmann-reflex during locomotion (Capaday and Stein, 1986; Yang and Stein, 1990). However, the neurons in the RG layer control the temporal aspect of the phase-dependent response as shown in the present study. As physiological experiments have shown, the feedback from muscle spindles contributed to the modulation of the muscle activity strength (Mayer et al., 2018) and the timing of the stance-to-swing and swing-to-stance transitions (Grillner and Rossignol, 1978; Hiebert et al., 1996; Akay et al., 2014). Therefore, the different layers of the CPG may explain the two different types of phase-dependency.

Limitations of Model

In our study, we used the activity-based neuron model (Ermentrout, 1994; Markin et al., 2010; Molkov et al., 2015; Danner et al., 2016, 2017). This neuron model does not show spiking because it omits the potassium and fast-type sodium currents. Instead, this used a persistent sodium current, which enables the neuron model to generate bursting. Ausborn et al. (2018) showed that an activity-based neuron model preserved the principal dynamic features of neural activities as a half-center CPG. Even though our model did not include potassium and fast-type sodium currents, it reproduced the phase-dependent response and contributed to analysis of its dynamic structure.

Interaction Between Body and Neural System During Adaptive Walking

In the present study, we focused on the phase-dependent response of the CPG activity during fictive locomotion. When animals walk, motor commands are sent to the leg muscles from the spinal CPG, and the CPG receives sensory signals from the leg nerves. While fictive locomotion is generated in an open-loop system, actual locomotion is generated in a closed-loop system. In addition to the analysis of fictive locomotion, in the future, we would like to investigate the entrainment mechanism through the dynamics of the CPG circuit, the body mechanical system, and the sensory system. Moreover, it has been suggested that the CPG consists of the RG and PF layers. While the RG layer determines the rhythm pattern of motor commands, the PF layer determines the spatial pattern (Rybak et al., 2006a). In the future, we would like to introduce the PF layer to our model to clarify further neural mechanisms of sensorimotor integration for adaptive locomotion.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

SF and SA developed the study design in consultation with KT, SD, and DY. SF performed simulations and analyzed the data in consultation with SA, KT, and SD. SF, SA, SD, and IR wrote the manuscript. All authors reviewed and approved the manuscript.

Funding

This research was supported in part by JSPS KAKENHI Grant-in-Aid for Young Scientists (B) JP16K16482, Grant-in-Aid for Scientific Research (B) JP15KT0015, and Grant-in-Aid for Scientific Research on Innovative Areas JP26120006 and also by the National Institutes of Health grants R01-NS100928 and R01-NS090919.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://[see PDF]/articles/10.3389/fnins.2019.01288/full#supplementary-material

References

Akay, T., Tourtellotte, W. G., Arber, S., and Jessell, T. M. (2014). Degradation of mouse locomotor pattern in the absence of proprioceptive sensory feedback. Proc. Natl. Acad. Sci. U.S.A. 111, 16877–16882. doi: 10.1073/pnas.1419045111

Aoi, S., Kondo, T., Hayashi, N., Yanagihara, D., Aoki, S., Yamaura, H., et al. (2013). Contributions of phase resetting and interlimb coordination to the adaptive control of hindlimb obstacle avoidance during locomotion in rats: a simulation study. Biol. Cybern. 107, 201–216. doi: 10.1007/s00422-013-0546-6

Aoi, S., Ogihara, N., Funato, T., Sugimoto, Y., and Tsuchiya, K. (2010). Evaluating functional roles of phase resetting in generation of adaptive human bipedal walking with a physiologically based model of the spinal pattern generator. Biol. Cybern. 102, 373–387. doi: 10.1007/s00422-010-0373-y

Aoi, S., Ohashi, T., Bamba, R., Fujiki, S., Tamura, D., Funato, T., et al. (2019). Neuromusculoskeletal model that walks and runs across a speed range with a few motor control parameter changes based on the muscle synergy hypothesis. Sci. Rep. 9:369. doi: 10.1038/s41598-018-37460-3

Ausborn, J., Snyder, A. C., Shevtsova, N. A., Rybak, I. A., and Rubin, J. E. (2018). State-dependent rhythmogenesis and frequency control in a half-center locomotor CPG. J. Neurophysiol. 119, 96–117. doi: 10.1152/jn.00550.2017

Capaday, C., and Stein, R. B. (1986). Amplitude modulation of the soleus H-reflex in the human during walking and standing. J. Neurosci. 6, 1308–1313. doi: 10.1523/JNEUROSCI.06-05-01308.1986

Danner, S. M., Shevtsova, N. A., Frigon, A., and Rybak, I. A. (2017). Computational modeling of spinal circuits controlling limb coordination and gaits in quadrupeds. eLife 6:e31050. doi: 10.7554/eLife.31050

Danner, S. M., Wilshin, S. D., Shevtsova, N. A., and Rybak, I. A. (2016). Central control of interlimb coordination and speed-dependent gait expression in quadrupeds. J. Physiol. 594, 6947–6967. doi: 10.1113/JP272787

Demir, S. S., Butera, R. J. Jr., DeFranceschi, A. A., Clark, J. W. Jr., and Byrne, J. H. (1997). Phase sensitivity and entrainment in a modeled bursting neuron. Biophys. J. 72, 579–594. doi: 10.1016/S0006-3495(97)78697-1

Duysens, J. (1977). Reflex control of locomotion as revealed by stimulation of cutaneous afferents in spontaneously walking premammillary cats. J. Neurophysiol. 40, 737–751. doi: 10.1152/jn.1977.40.4.737

Ermentrout, B. (1994). Reduction of conductance-based models with slow synapses to neural nets. Neural Comput. 6,679–695. doi: 10.1162/neco.1994.6.4.679

Frigon, A., Sirois, J., and Gossard, J. P. (2010). Effects of ankle and hip muscle afferent inputs on rhythm generation during fictive locomotion. J. Neurophysiol. 103, 1591–1605. doi: 10.1152/jn.01028.2009

Fujiki, S., Aoi, S., Funato, T., Sato, Y., Tsuchiya, K., and Yanagihara, D. (2018). Adaptive hindlimb split-belt treadmill walking in rats by controlling basic muscle activation patterns via phase resetting. Sci. Rep. 8:17341. doi: 10.1038/s41598-018-35714-8

Grillner, S. (1981). “Control of locomotion in bipeds, tetrapods, and fish,” in Handbook of Physiology, Vol. II, eds J. M. Brookhart and V. B. Mountcastle (Bethesda, MD: American Physiological Society), 1179–1236.

Grillner, S., and Rossignol, S. (1978). On the initiation of the swing phase of locomotion in chronic spinal cats. Brain Res. 146, 269–277. doi: 10.1016/0006-8993(78)90973-3

Guertin, P., Angel, M. J., Perreault, M.-C., and McCrea, D. A. (1995). Ankle extensor group I afferents excite extensors throughout the hindlimb during fictive locomotion in the cat. J. Physiol. 487, 197–209. doi: 10.1113/jphysiol.1995.sp020871

Hiebert, G. W., Whelan, P. J., Prochazka, A., and Pearson, K. G. (1996). Contribution of hind limb flexor muscle afferents to the timing of phase transitions in the cat step cycle. J. Neurophysiol. 75, 1126–1137. doi: 10.1152/jn.1996.75.3.1126

Ivanenko, Y. P., Poppele, R. E., and Lacquaniti, F. (2004). Five basic muscle activation patterns account for muscle activity during human locomotion. J. Physiol. 556, 267–282. doi: 10.1113/jphysiol.2003.057174

Ivanenko, Y. P., Poppele, R. E., and Lacquaniti, F. (2006). Motor control programs and walking. Neuroscientist 12, 339–348. doi: 10.1177/1073858406287987

Markin, S. N., Klishko, A. N., Shevtsova, N. A., Lemay, M. A., Prilutsky, B. I., and Rybak, I. A. (2010). Afferent control of locomotor CPG: insights from a simple neuromechanical model. Ann. N. Y. Sci. 1198, 21–34. doi: 10.1111/j.1749-6632.2010.05435.x

Mayer, W. P., Murray, A. J., Brenner-Morton, S., Jessell, T. M., Tourtellotte, W. G., and Akay, T. (2018). Role of muscle spindle feedback in regulating muscle activity strength during walking at different speed in mice. J. Neurophysiol. 120, 2484–2497. doi: 10.1152/jn.00250.2018

McCrea, D. A. (2001). Spinal circuitry of sensorimotor control of locomotion. J. Physiol. 533, 41–50. doi: 10.1111/j.1469-7793.2001.0041b.x

McCrea, D. A., and Rybak, I. A. (2007). Modeling the mammalian locomotor CPG: insights from mistakes and perturbations. Progr. Brain Res. 165, 235–253. doi: 10.1016/S0079-6123(06)65015-2

Molkov, Y. I., Bacak, B. J., Talpalar, A. E., and Rybak, I. A. (2015). Mechanisms of left-right coordination in mammalian locomotor pattern generation circuits: a mathematical modeling view. PLoS Comput. Biol. 11:e1004270. doi: 10.1371/journal.pcbi.1004270

Orlovsky, G. N., Deliagina, T., and Grillner, S. (1999). Neuronal Control of Locomotion: from Mollusc to Man. New York, NY: Oxford University Press. doi: 10.1093/acprof:oso/9780198524052.001.0001

Pearson, K. G. (2004). Generating the walking gait: role of sensory feedback. Prog. Brain Res. 143, 123–129. doi: 10.1016/S0079-6123(03)43012-4

Perreault, M., Angel, M. J., Guertin, P., and McCrea, D. A. (1995). Effects of stimulation of hindlimb flexor group II afferents during fictive locomotion in the cat. J. Physiol. 487, 211–220. doi: 10.1113/jphysiol.1995.sp020872

Rossignol, S. (1996). “Neural control of stereotypic limb movements,” in Handbook of Physiology, eds L. B. Rowell and J. Shepherd (Bethesda, MD: American Physiological Society, 173–216.

Rossignol, S., Dubuc, R., and Gossard, J. P. (2006). Dynamic sensorimotor interactions in locomotion. Physiol. Rev. 86, 89–154. doi: 10.1152/physrev.00028.2005

Rybak, I. A., Shevtsova, N. A., Lafreniere-Roula, M., and McCrea, D. A. (2006a). Modelling spinal circuitry involved in locomotor pattern generation: insights from deletions during fictive locomotion. J. Physiol. 577, 617–639. doi: 10.1113/jphysiol.2006.118703

Rybak, I. A., Stecina, K., Shevtsova, N. A., and McCrea, D. A. (2006b). Modelling spinal circuitry involved in locomotor pattern generation: insights from effects of afferent stimulation. J. Physiol. 577, 641–658. doi: 10.1113/jphysiol.2006.118711

Schomburg, E. D., Petersen, N., Barajon, I., and Hultborn, H. (1998). Flexor reflex afferents reset the step cycle during fictive locomotion in the cat. Exp. Brain Res. 122, 339–350. doi: 10.1007/s002210050522

Spardy, L. E., Markin, S. N., Shevtsova, N. A., Prilutsky, B. I., Rybak, I. A., and Rubin, J. E. (2011a). A dynamical systems analysis of afferent control in a neuromechanical model of locomotion. I. Rhythm generation. J. Neural. Eng. 8:065003. doi: 10.1088/1741-2560/8/6/065003

Spardy, L. E., Markin, S. N., Shevtsova, N. A., Prilutsky, B. I., Rybak, I. A., and Rubin, J. E. (2011b). A dynamical systems analysis of afferent control in a neuromechanical model of locomotion. II. Phase asymmetry. J. Neural. Eng. 8:065004. doi: 10.1088/1741-2560/8/6/065004

Stecina, K., Quevedo, J., and McCrea, D. A. (2005). Parallel reflex pathways from flexor muscle afferents evoking resetting and flexion enhancement during fictive locomotion and scratch in the cat. J. Physiol. 569, 275–290. doi: 10.1113/jphysiol.2005.095505

Whelan, P. J., Hiebert, G. W., and Pearson, K. G. (1995). Stimulation of the group I extensor afferents prolongs the stance phase in walking cats. Exp. Brain Res. 103, 20–30. doi: 10.1007/BF00241961

Yang, J. F., and Stein, R. B. (1990). Phase-dependent reflex reversal in human leg muscles during walking. J. Neurophysiol. 63, 1109–1117. doi: 10.1152/jn.1990.63.5.1109