Contribution of Afferent Feedback to Adaptive Hindlimb Walking in Cats: A Neuromusculoskeletal Modeling Study

Mammalian locomotion is generated by central pattern generators (CPGs) in the spinal cord, which produce alternating flexor and extensor activities controlling the locomotor movements of each limb. Afferent feedback signals from the limbs are integrated by the CPGs to provide adaptive control of locomotion. Responses of CPG-generated neural activity to afferent feedback stimulation have been previously studied during fictive locomotion in immobilized cats. Yet, locomotion in awake, behaving animals involves dynamic interactions between central neuronal circuits, afferent feedback, musculoskeletal system, and environment. To study these complex interactions, we developed a model simulating interactions between a half-center CPG and the musculoskeletal system of a cat hindlimb. Then, we analyzed the role of afferent feedback in the locomotor adaptation from a dynamic viewpoint using the methods of dynamical systems theory and nullcline analysis. Our model reproduced limb movements during regular cat walking as well as adaptive changes of these movements when the foot steps into a hole. The model generates important insights into the mechanism for adaptive locomotion resulting from dynamic interactions between the CPG-based neural circuits, the musculoskeletal system, and the environment.

1 Introduction

Mammalian locomotion is generated by the central pattern generators (CPGs) located in the spinal cord, which control the movements of each limb (Grillner, 1981; Rossignol, 1996; Orlovsky et al., 1999; Rossignol et al., 2006). CPGs produce the basic locomotor rhythm and create alternating flexor and extensor motoneuron (Mn) activities. Furthermore, they integrate afferent feedback signals to achieve adaptive locomotion.

CPGs can operate without afferent feedback and continuous electrical stimulation of the midbrain locomotor region in immobilized decerebrate cats produces “fictive locomotion” consisting of rhythmic alternating activation of flexor and extensor Mns similar to that occurring during normal locomotion in intact animals (Rossignol, 1996). Such fictive locomotor preparations have been used to investigate the mechanism for the adaptive regulation of locomotor patterns by somatosensory afferent feedback. These studies have shown that stimulation of flexor and extensor sensory afferents can delay or advance flexion-to-extension or extension-to-flexion transitions depending on the timing of the afferent stimulation (Guertin et al., 1995; Perreault et al., 1995; McCrea, 2001; Stecina et al., 2005). In our previous modeling work (Fujiki et al., 2019), we used a half-center type CPG model and analyzed its responses to afferent stimulation using dynamic systems theory based on nullclines. This previous study was limited to the consideration of only neural responses without interaction between the neural and musculoskeletal systems.

However, the mammalian locomotion is a complex phenomenon involving dynamic interactions between the neural circuits, the musculoskeletal system, and the environment. To investigate these interactions and the roles of afferent feedback in adaptive locomotion, investigators studied locomotion in both intact and decerebrate or spinal animals walking on a treadmill by applying various disturbances or introducing holes and obstacles on the treadmill to disturb normal locomotion (Hiebert et al., 1994; Lam and Pearson, 2001; Drew et al., 2002). Although these studies have revealed adaptive responses to such disturbances, the adaptation mechanisms based on the interactions between the neural system, musculoskeletal system, and environment remain poorly understood.

Here, we extended our previous models (Markin et al., 2010, 2016), which integrated a two-level half-center CPG and a musculoskeletal model of the cat hindlimbs and simulated steady walking, to investigate the mechanism of locomotor adaptation in response to perturbations from a dynamical viewpoint. We determined the necessary model parameters through optimization that allows the model to reproduce regular walking movements on a treadmill. We considered locomotion with the presence of holes in the walking surface, which evoked locomotor disturbances that allowed us to analyze the roles of afferent feedback from flexor and extensor muscles to stabilize locomotion. Particularly, the model realistically reproduced adaptive changes in locomotor characteristics observed when the paw of a cat hindlimb stepped into a hole (Hiebert et al., 1994). The model suggests an afferent feedback-based mechanism for locomotor adaptation, which we analyzed based on dynamical systems theory methods using nullclines (Fujiki et al., 2019). Our simulation and analysis provide important insights into the mechanisms for adaptive locomotion based on dynamic interactions between the neural system, the musculoskeletal system, and the environment.

2 Model

Figure 1 shows our neuromusculoskeletal model, which consists of a musculoskeletal model of a cat hindlimb and a spinal CPG model to drive the musculoskeletal model.

FIGURE 1

FIGURE 1. Neuromusculoskeletal model of a cat hindlimb composed of a neural network and a musculoskeletal model. The neural network model consists of a two-level central pattern generator (CPG) with rhythm generating (RG) and pattern formation (PF) circuits and motoneurons (Mns). This network generates motor commands through Mns and receives afferent feedback signals from the musculoskeletal model (orange arrows). The musculoskeletal model, which consists of three rigid links (black lines) and seven muscles (red and orange lines), walks on a treadmill.

2.1 Musculoskeletal Model

The skeletal model is two-dimensional and consists of three rigid links representing the thigh, crus, and foot. These links are connected by the knee and ankle joints, and the hip joint is fixed above a treadmill. The model walks on the treadmill with a belt speed of 0.4 m/s based on Prilutsky et al. (2016). When the thigh, crus, and foot are in a straight line and parallel to the vertical line, the hip angle is 135° and the knee and ankle angles are both 180°. The joint angles increase as the joints are extending. The contact between the limb tip and treadmill were modeled using viscoelastic elements. We derived the equations of motion for the skeletal model using Lagrangian equations, where we used the same physical parameters for the skeletal model as those in Ekeberg and Pearson (2005), and solved the equations numerically using the fourth-order Runge-Kutta method with a time step of 0.04 ms.

The skeletal model is driven by seven Hill-type muscles including five uni-articular muscles which are: hip flexor (iliopsoas, IP), hip extensor (gluteus maximus, GM), knee extensor (vastus lateralis, VL), ankle flexor (tibialis anterior, TA), ankle extensor (soleus, SO), and two bi-articular muscles including hip extensor/knee flexor (biceps femoris, BF) and knee flexor/ankle extensor (gastrocnemius, GA). We assumed that the moment arms of all muscles are constant. Each muscle generates the muscle tension through contractile and passive elements. The muscle model consisting of contractile and passive elements is based on the same description and parameters in Ekeberg and Pearson (2005). Specifically, the muscle force $F_m\left(m\in \left\{\mathrm{M}\right\}=\left\{\mathrm{I}\mathrm{P},\mathrm{G}\mathrm{M},\mathrm{V}\mathrm{L},\mathrm{T}\mathrm{A},\mathrm{S}\mathrm{O},\mathrm{B}\mathrm{F},\mathrm{G}\mathrm{A}\right\}\right)$ is given by

$F_m=F_m^{\text{max}}\left(a_m{F}_m^{\mathrm{l}}{F}_m^{\mathrm{v}}+F_m^{\mathrm{p}}\right)(1)$

where $F_m^{\text{max}}$ is the maximum isometric force, a_m is the muscle activation $\left(0\le a_m\le 1\right)$, $F_m^{\mathrm{l}}$ is the force-length relationship, $F_m^{\mathrm{v}}$ is the force-velocity relationship, and $F_m^{\mathrm{p}}$ is the passive component. The muscle lengths were normalized by $l_m^{\text{max}}$, which was set so that all uni-articular muscles had a length of 85% of $l_m^{\text{max}}$ and all bi-articular muscles were at 75% at a neutral posture with the hip at 65°, the knee at 90°, and the ankle at 100°. In addition, 2° of joint motion corresponded to 1% of muscle length change, except for the GA muscle, where 1.5° at the ankle or 4.5° at the knee was required. The muscle contractile velocities were normalized by $l_m^{\text{max}}$ as well.

The muscle activation a_m $(m\in \left\{\mathrm{M}\right\})$ determines the muscle tension generated by the contractile element and the dynamics of a_m is given by a low-pass filter (Yakovenko et al., 2004) as follows:

${\dot{a}}_m+\frac{1}{{\tau }_{\text{act}}}\left(\frac{{\tau }_{\text{act}}}{{\tau }_{\text{dact}}}+\left[1-\frac{{\tau }_{\text{act}}}{{\tau }_{\text{dact}}}\right]u_m\right)a_m=\frac{1}{{\tau }_{\text{act}}}{u}_m(2)$

where τ_act and τ_dact are activation and deactivation time constants (20 and 32 ms, respectively), and u_m is the motor command determined from the activities of the corresponding Mn of the CPG model.

2.2 CPG Model

The locomotor CPG has been suggested to consist of hierarchical networks, which include rhythm generator (RG) and pattern formation (PF) networks (Rybak et al., 2006a,b). The RG network generates the rhythmic activities while the PF network generates the spatiotemporal patterns of motor commands. For the RG model, we used two neuron populations of flexor and extensor centers (RG-F and RG-E), which receive a supraspinal drive, and two populations of inhibitory interneurons (In-F and In-E), which provide mutual inhibition between the RG-F and RG-E centers. The PF was modeled using two neuron populations of flexor and extensor centers (PF-F and PF-E). The PF-F and PF-E neuron populations receive the excitatory input from the RG-F and RG-E neuron populations, and inhibitory input from the In-E and In-F neuron populations, respectively. Seven Mn populations provide activation for each muscle in the musculoskeletal model (Mn-m, $m\in \left\{\mathrm{M}\right\}$). The Mns of flexor muscles Mn-IP, Mn-TA, and Mn-BF receive excitatory input from the PF-F neuron populations, while those of extensor muscles Mn-GM, Mn-VL, Mn-SO, and Mn-GA receive the excitatory input from the PF-E neuron populations. Synaptic interactions between all neuron populations are shown in Figure 1.

Each population is 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 is characterized by the membrane potential V_i for $i\in \left\{\mathrm{R}\mathrm{G}\right\}$, $\left\{\mathrm{I}\mathrm{n}\right\}$, $\left\{\mathrm{P}\mathrm{F}\right\}$, and $\left\{\mathrm{M}\mathrm{n}\right\}$, where {RG} ={RG-F, RG-E}, {In} ={In-F, In-E}, {PF} ={PF-F, PF-E}, and {Mn}$=\{\mathrm{M}\mathrm{n}‐m|m\in \left\{\mathrm{M}\right\}\}$. The RG, PF, and Mn neurons incorporate a persistent (slowly-inactivating) sodium current that defines the intrinsic rhythmogenic properties of these neurons. The intrinsic oscillations in the RG, PF, and Mn neurons depend on the variable h_i ($i\in \left\{\mathrm{R}\mathrm{G}\right\}$, $\left\{\mathrm{P}\mathrm{F}\right\}$, $\left\{\mathrm{M}\mathrm{n}\right\}$) that defines the slow inactivation of the persistent sodium channel. The RG-F and RG-E neurons can produce rhythmic activities. However, if uncoupled, the RG-E neuron is in the tonic regime due to the supraspinal drive and produces sustained activity. Rhythmic oscillations of the RG neurons are defined by the RG-F neuron, which provides rhythmic inhibition of the RG-E neuron through the In-F neuron. The supraspinal drive to the RG-F neuron determines the oscillation frequency. When the PF and Mn neurons are uncoupled, they do not produce rhythmic activities due to the relatively low maximum conductance of the sodium current. Instead, these neurons produce rhythmic activities through the excitatory inputs from the correspondent RG neurons.

For the state variable for this model, we used $\mathit{V}={[V_{\left\{\mathrm{R}\mathrm{G}\right\}},V_{\left\{\mathrm{I}\mathrm{n}\right\}},V_{\left\{\mathrm{P}\mathrm{F}\right\}},V_{\left\{\mathrm{M}\mathrm{n}\right\}}]}^{\mathrm{T}}$ and $\mathit{h}={[h_{\left\{\mathrm{R}\mathrm{G}\right\}},h_{\left\{\mathrm{P}\mathrm{F}\right\}},h_{\left\{\mathrm{M}\mathrm{n}\right\}}]}^{\mathrm{T}}$. The dynamics of the membrane potential V_i is described as

$C{\dot{V}}_i=\left\{\begin{array}{cc}-I_{\mathrm{N}\mathrm{a}\mathrm{P}}\left(V_i,h_i\right)-I_{\mathrm{L}\mathrm{e}\mathrm{a}\mathrm{k}}\left(V_i\right)-I_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{E}}^i\left(\mathit{V}\right)-I_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{I}}^i\left(\mathit{V}\right)\hfill & i\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{P}\mathrm{F}\right\},\left\{\mathrm{M}\mathrm{n}\right\},\hfill \\ -I_{\mathrm{L}\mathrm{e}\mathrm{a}\mathrm{k}}\left(V_i\right)-I_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{E}}^i\left(\mathit{V}\right)-I_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{I}}^i\left(\mathit{V}\right)\hfill & i\in \left\{\mathrm{I}\mathrm{n}\right\},\hfill \end{array}\right.(3)$

where C is the membrane capacitance, I_NaP is the persistent sodium current, I_Leak is the leakage current, and $I_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{E}}^i$ and $I_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{I}}^i$ are the respective currents in excitatory synapses and inhibitory synapses. The ionic current I_NaP and leakage current I_Leak are described as

$\begin{array}{cc}\hfill I_{\mathrm{N}\mathrm{a}\mathrm{P}}\left(V_i,h_i\right)& ={\widehat{g}}_{\mathrm{N}\mathrm{a}\mathrm{P}}^i{m}_{\text{NaP}}\left(V_i\right)h_i\left\{V_i-E_{\text{Na}}\right\}i\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{P}\mathrm{F}\right\},\left\{\mathrm{M}\mathrm{n}\right\},\hfill \\ \hfill I_{\mathrm{L}\mathrm{e}\mathrm{a}\mathrm{k}}\left(V_i\right)& ={\widehat{g}}_{\text{Leak}}^i\left\{V_i-E_{\text{Leak}}^i\right\}i\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{I}\mathrm{n}\right\},\left\{\mathrm{P}\mathrm{F}\right\},\left\{\mathrm{M}\mathrm{n}\right\},\hfill \end{array}(4)$

where ${\widehat{g}}_{\mathrm{N}\mathrm{a}\mathrm{P}}^i$ and ${\widehat{g}}_{\mathrm{L}\mathrm{e}\mathrm{a}\mathrm{k}}^i$ are the maximum conductances of the corresponding currents, and E_Na and $E_{\mathrm{L}\mathrm{e}\mathrm{a}\mathrm{k}}^i$ the reversal potentials. In addition, m_NaP is the activation of the sodium channel of the RG, PF, and Mn neurons and is described as

$m_{\text{NaP}}\left(V_i\right)=\frac{1}{1+\exp \left(-\frac{V_i+40.0}{6.0}\right)}i\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{P}\mathrm{F}\right\},\left\{\mathrm{M}\mathrm{n}\right\}.(5)$

The dynamics of the inactivation of the sodium channel h_i of the RG, PF, and Mn neurons is given by

$\tau \left(V_i\right){\dot{h}}_i=h_{\infty }\left(V_i\right)-h_{i}i\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{P}\mathrm{F}\right\},\left\{\mathrm{M}\mathrm{n}\right\},(6)$

where

$h_{\infty }\left(V_i\right)=\frac{1}{1+\exp \left(\frac{V_i+45.0}{4.0}\right)},(7)$

$\tau \left(V_i\right)=320+\frac{320}{\cosh \left(\frac{V_i+35.0}{15.0}\right)}\mathrm{m}\mathrm{s}i\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{P}\mathrm{F}\right\},\left\{\mathrm{M}\mathrm{n}\right\}$

The currents generated by the synapses $I_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{E}}^i$ and $I_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{I}}^i$ are given by

$\begin{array}{cc}\hfill I_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{E}}^i& ={\widehat{g}}_{\text{SynE}}\left\{V_i-E_{\text{SynE}}\right\}\left\{\sum _{j\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{I}\mathrm{n}\right\},\left\{\mathrm{P}\mathrm{F}\right\}}{\alpha }_{ij}f\left(V_j\right)+{\gamma }_{i}d+s_i\right\},\hfill \\ \hfill I_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{I}}^i& ={\widehat{g}}_{\text{SynI}}\left\{V_i-E_{\text{SynI}}^i\right\}\left\{\sum _{j\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{I}\mathrm{n}\right\},\left\{\mathrm{P}\mathrm{F}\right\}}{\beta }_{ij}f\left(V_j\right)\right\}i\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{I}\mathrm{n}\right\},\left\{\mathrm{P}\mathrm{F}\right\},\left\{\mathrm{M}\mathrm{n}\right\},\hfill \end{array}(8)$

where ${\widehat{g}}_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{E}}$ and ${\widehat{g}}_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{I}}$ are the maximum conductances of the corresponding currents; E_SynE and $E_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{I}}^i$ are the reversal potentials of the corresponding currents; d is the tonic drive from the supraspinal region; s_i is the afferent feedback from the musculoskeletal model (as determined in Section 4.1); and α_ij, β_ij, and γ_i are the weight coefficients, where β_ij = 0 for $i\in \left\{\mathrm{M}\mathrm{n}\right\}$ and γ_i = 0 for $i\in \left\{\mathrm{P}\mathrm{F}\right\}$ and $\left\{\mathrm{M}\mathrm{n}\right\}$. Moreover, the output function f translates V into the integrated population activity and is given by

$f\left(V_i\right)=\left\{\begin{array}{cc}0\hfill & V_i<V_{\mathrm{t}\mathrm{h}},\hfill \\ \left(V_i-V_{\mathrm{t}\mathrm{h}}\right)/\left(V_{\mathrm{m}\mathrm{a}\mathrm{x}}-V_{\mathrm{t}\mathrm{h}}\right)\hfill & V_{\mathrm{t}\mathrm{h}}<V_i<V_{\mathrm{m}\mathrm{a}\mathrm{x}},\hfill \\ 1\hfill & V_{\mathrm{m}\mathrm{a}\mathrm{x}}<V_i,\hfill \end{array}\right.(9)$

where V_th and V_max are the lower and upper threshold potentials, respectively. The motor command $u_m(m\in \left\{\mathrm{M}\right\})$ is given by $u_m=f\left(V_{\text{Mn}‐m}\right)$. Based on Fujiki et al. (2019) and Markin et al. (2010), we determined the parameters for the CPG model (see Appendix A) except for α_ij, β_ij, and γ_i for $i\in \left\{\mathrm{M}\mathrm{n}\right\}$, which we determined through optimization, as described in Section 4.2.

2.3 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. After the CPG model was integrated with the musculoskeletal model to achieve steady walking, we used a nullcline-based method, as in our previous work (Fujiki et al., 2019), to investigate the mechanism of the response of the CPG model to a disturbance during steady walking. Specifically, the state of the CPG model is given by $\left(\mathit{V},\mathit{h}\right)$, and the nullclines for the RG neurons are given by

$\left\{\begin{array}{c}{N}_i^{\mathrm{V}}=\left\{\left(\mathit{V},\mathit{h}\right)|{\dot{V}}_i=0\right\},\hfill \\ N_i^{\mathrm{h}}=\left\{\left(\mathit{V},\mathit{h}\right)|{\dot{h}}_i=0\right\},\hfill \end{array}\right.i\in \left\{\mathrm{R}\mathrm{G}\right\}.(10)$

To clarify the dynamics of each RG neuron, we focused on the V_i-h_i space $(i\in \left\{\mathrm{R}\mathrm{G}\right\})$ for the nullclines by assuming that the other variables $V_j(j\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{I}\mathrm{n}\right\},\left\{\mathrm{P}\mathrm{F}\right\},\left\{\mathrm{M}\mathrm{n}\right\},j\ne i)$ and $h_k(k\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{P}\mathrm{F}\right\},\left\{\mathrm{M}\mathrm{n}\right\},k\ne i)$ are stably oscillating during steady walking. Therefore, we modify $N_i^{\mathrm{V}}$ and $N_i^{\mathrm{h}}$ in Eq. 10 as

$\begin{array}{l}\hfill \begin{array}{c}\hfill \left\{\begin{array}{c}{\widehat{N}}_i^{\mathrm{V}}=\left\{\left(V_i,h_i\right)|{\dot{V}}_i=0,V_j=V_j^{\ast },h_k=h_k^{\ast }\right\},\hfill \\ {\widehat{N}}_i^{\mathrm{h}}=\left\{\left(V_i,h_i\right)|{\dot{h}}_i=0,V_j=V_j^{\ast },h_k=h_k^{\ast }\right\},\hfill \end{array}\right.\\ \hfill i\in \left\{\mathrm{R}\mathrm{G}\right\},j\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{I}\mathrm{n}\right\},\left\{\mathrm{P}\mathrm{F}\right\},\left\{\mathrm{M}\mathrm{n}\right\},j\ne i,k\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{P}\mathrm{F}\right\},\left\{\mathrm{M}\mathrm{n}\right\},k\ne i,\end{array}\end{array}(11)$

where x^∗ indicates x for the stable oscillation during steady walking.

While ${\widehat{N}}_i^{\mathrm{h}}(i\in \left\{\mathrm{R}\mathrm{G}\right\})$ has a sigmoid shape and does not change with time, ${\widehat{N}}_i^{\mathrm{V}}(i\in \left\{\mathrm{R}\mathrm{G}\right\})$ has a cubic curve shape and changes by afferent feedback and input from other neurons. Specifically, ${\widehat{N}}_i^{\mathrm{V}}$ mainly has the two different situations shown in Figures 2A,B, one with two different inflection points at which the sign of the slope changes, and one with a monotonic variation. In the former case (Figure 2A), the intersection of ${\widehat{N}}_i^{\mathrm{V}}$ and ${\widehat{N}}_i^{\mathrm{h}}$ is a saddle equilibrium point. Because the time constant for the dynamics of $V_i(i\in \left\{\mathrm{R}\mathrm{G}\right\})$ is smaller than that for $h_i(i\in \left\{\mathrm{R}\mathrm{G}\right\})$, the following two features are present (Fujiki et al., 2019; Molkov et al., 2015; Spardy et al., 2011): 1. near where ${\widehat{N}}_i^{\mathrm{V}}$ has positive slope, the state (V_i, h_i) is slowly attracted to the inflection point along ${\widehat{N}}_i^{\mathrm{V}}$ (slow dynamics) and 2. near the inflection points, the state (V_i, h_i) jumps to the opposite part of ${\widehat{N}}_i^{\mathrm{V}}$ with a positive slope (fast dynamics). In the latter case (Figure 2B), the intersection of ${\widehat{N}}_i^{\mathrm{V}}$ and ${\widehat{N}}_i^{\mathrm{h}}$ is a stable node equilibrium point having the following two features (Spardy et al., 2011; Fujiki et al., 2019): 1. when the state (V_i, h_i) is away from ${\widehat{N}}_i^{\mathrm{V}}$, it is quickly attracted to ${\widehat{N}}_i^{\mathrm{V}}$ (fast dynamics), and 2. near ${\widehat{N}}_i^{\mathrm{V}}$, the state (V_i, h_i) is slowly attracted to the stable node along ${\widehat{N}}_i^{\mathrm{V}}$ (slow dynamics). As described above, the switch between fast and slow dynamics depends on the relationship between the state and the nullcline, and adaptive responses are achieved through the changes in ${\widehat{N}}_i^{\mathrm{V}}$ by afferent feedback.

FIGURE 2

FIGURE 2. Nullclines ${\widehat{N}}_i^{\mathrm{V}}$ and ${\widehat{N}}_i^{\mathrm{h}}$ for $i\in \left\{\mathrm{R}\mathrm{G}\right\}$ of RG neurons in the phase plane and fast and slow dynamics that the state (V_i, h_i) for $i\in \left\{\mathrm{R}\mathrm{G}\right\}$ follows. (A) ${\widehat{N}}_i^{\mathrm{V}}$ has two different inflection points and intersection with ${\widehat{N}}_i^{\mathrm{h}}$ is saddle equilibrium point. (B) ${\widehat{N}}_i^{\mathrm{V}}$ changes monotonically and intersection with ${\widehat{N}}_i^{\mathrm{h}}$ is stable node equilibrium point. Bold and thin arrows represent fast and slow dynamics, respectively.

3 Verifying the CPG Model by Phase-Dependent Response in Fictive Locomotion

Even when our CPG model is separated from the musculoskeletal model (or receives no feedback signals), it produces rhythmic activities and exhibits stable oscillations, as shown in Figure 3. Based on these oscillations, we define the active phase for each neuron as the time interval during which the neuron’s potential is higher than V_th, and the silent phase as the time interval when the potential is lower than V_th. We also define the cycle period T as the time interval between two consecutive onsets of the active phase of the PF-F neuron and the phase of the oscillation as ϕ = 2πt/T ∈ [0, 2π).

FIGURE 3

FIGURE 3. Change in membrane potential of PF-F neuron by stimulation: (A) No stimulation. (B) Stimulation to the flexor side during the active phase, which increases the duration of the current active phase and cycle (T′ > T). (C) Stimulation to the flexor side during the silent phase, which initiates the active phase and decreases the cycle duration (T′ < T). Red bars indicate the application of stimulation. Gray regions indicate active phases.

To verify our CPG model, we investigated the phase-dependent response of the CPG activities based on previous studies (Demir et al., 1997; Fujiki et al., 2019) and compared the results with those obtained by experiments with fictive locomotion in cats (Duysens, 1977; Schomburg et al., 1998). Specifically, we used only the CPG model (RG, In, and PF neurons) separated from the musculoskeletal model. After the oscillation of the CPG model stabilized, we applied a 200-ms stimulus to the flexor (RG-F, In-F, and PF-F) or extensor (RG-E, In-E, and PF-E) neurons, where the intensity of the stimulation was set as follows: s_RG-F = s_In-F = s_PF-F = 0.2 and s_RG-E = s_In-E = s_PF-E = 0 for the stimulation to the flexor side and s_RG-F = s_In-F = s_PF-F = 0.0 and s_RG-E = s_In-E = s_PF-E = 0.2 for the stimulation to the extensor side in Eq. 8. Suppose that the neuron activity is perturbed by stimulation with a phase ϕ_s ∈ [0, 2π) and the period of the PF-F neuron changes from T to T′(ϕ_s), as shown in Figure 3. To explain the phase shift of the neuron activity in response to the stimulation, we define

$\mathrm{\Delta }\left({\varphi }_{\mathrm{s}}\right)=2\pi \frac{T^{\prime }\left({\varphi }_{\mathrm{s}}\right)-T}{T}(12)$

Figure 4A shows the phase shift Δ of the PF-F neuron activity after the stimulation of sensory inputs on the flexor side at ϕ_s. When the stimulation was applied during the silent phase of the PF-F neuron (2.70 ≤ ϕ_s ≤ 2π), it caused an earlier transition to the active phase, and this advanced start decreased with ϕ_s. In contrast, almost no phase shift occurred when the stimulation was applied at the beginning of the active phase of the PF-F neuron (0 ≤ ϕ_s ≤ 0.44). However, the neuron activity was delayed by the stimulation during the middle and end of the active phase (0.44 ≤ ϕ_s ≤ 2.70). These trends were similar to those observed during fictive locomotion in spinal cats (Schomburg et al., 1998; Frigon et al., 2010), as shown in Figure 4A.

FIGURE 4

FIGURE 4. Phase-dependent response of CPG model to a stimulation. (A) Response of PF-F neuron in our model to stimulation of the flexor side compared with the response against flexor muscle stimulation during fictive locomotion in spinal cats [adapted from Schomburg et al. (1998)], where the dotted lines are approximate functions using a first order polynomial for each flexion and extension phase. (B) Response of PF-F neuron of our model to stimulation of the extensor side compared with the response against extensor muscle stimulation during fictive locomotion in decerebrate cats [adapted from Duysens (1977)], where the dotted line is an approximate function using an eighth order polynomial. R is the correlation coefficient between the active/silent phase of the simulation results and flexion/extension phase of the approximate function.

Figure 4B shows Δ after stimulation of the extensor side. The active and silent phase of the PF-F neuron corresponds to the silent and active phase, respectively, of the PF-E neuron. The neuron activity was advanced at the middle of the silent phase of the PF-E neuron (0.63 ≤ ϕ_s ≤ 2.64) and was delayed at the end of the active phase of the PF-E neuron (4.27 ≤ ϕ_s ≤ 6.16). The response of the stimulation of the extensor side was qualitatively similar to that for the flexor side. Moreover, the trends were consistent with those observed during fictive locomotion in decerebrate cats (Duysens, 1977; Schomburg et al., 1998; Frigon et al., 2010), as shown in Figure 4B. These results verify the validity of our CPG model.

4 Determining the Motor Control Model by Optimization to Reproduce Normal Walking

We integrated the CPG and musculoskeletal models to determine the remaining parameters for the CPG model through an optimization to produce normal walking on the treadmill.

4.1 Afferent Feedback From the Musculoskeletal Model

We determined the afferent feedback s_i in Eq. 8 based on Markin et al. (2010) as follows:

$s_i=\sum _{m\in \left\{\mathrm{M}\right\}}\left(k_{im}^{\mathrm{v}}{\left(v_m^{\mathrm{norm}}\right)}^{0.6}+k_{im}^{\mathrm{d}}{d}_m^{\mathrm{norm}}+k_{im}^{\mathrm{f}}{F}_m^{\mathrm{norm}}\right)i\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{I}\mathrm{n}\right\},\left\{\mathrm{P}\mathrm{F}\right\},(13)$

where

$v_m^{\text{norm}}=\left\{\begin{array}{cc}{v}_m/l_m^{\text{max}}\hfill & v_m>0,\hfill \\ 0\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e};\hfill \end{array}\right.(14)$

$d_m^{\text{norm}}=\left\{\begin{array}{cc}\left(l_m-l_{\text{th}}\right)/l_m^{\text{max}}\hfill & l_m>l_{\text{th}},\hfill \\ 0\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e};\hfill \end{array}\right.(15)$

$F_m^{\text{norm}}=\left\{\begin{array}{cc}{F}_m/F_m^{\text{max}}\hfill & F_m>0,\hfill \\ 0\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e};\hfill \end{array}\right.(16)$

$l_{\text{th}}=0.9l_m^{\text{max}};(17)$

v_m, l_m, and F_m are the velocity, length, and force of muscle m, respectively, and $k_{im}^{\mathrm{v}}$, $k_{im}^{\mathrm{d}}$, and $k_{im}^{\mathrm{f}}$ are the corresponding weight coefficients. The three terms on the right-hand side of Eq. 13 represent the velocity, length, and force feedback from muscle m. The coefficients $k_{im}^{\mathrm{v}}$, $k_{im}^{\mathrm{d}}$, and $k_{im}^{\mathrm{f}}$ are given as follows: for the flexor side (i = RG-F, In-F, PF-F),

$\left(k_{im}^{\mathrm{v}},k_{im}^{\mathrm{d}},k_{im}^{\mathrm{f}}\right)=\left\{\begin{array}{cc}\left(k_{\mathrm{F}}^{\mathrm{v}},k_{\mathrm{F}}^{\mathrm{d}},k_{\mathrm{F}}^{\mathrm{f}}\right)\hfill & m=\mathrm{I}\mathrm{P},\mathrm{T}\mathrm{A},\mathrm{B}\mathrm{F},\hfill \\ \left(\mathrm{0,0,0}\right)\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e};\hfill \end{array}\right.(18)$

for the extensor side (i = RG-E, In-E, PF-E),

$\left(k_{im}^{\mathrm{v}},k_{im}^{\mathrm{d}},k_{im}^{\mathrm{f}}\right)=\left\{\begin{array}{cc}\left(k_{\mathrm{E}}^{\mathrm{v}},k_{\mathrm{E}}^{\mathrm{d}},k_{\mathrm{E}}^{\mathrm{f}}\right)\hfill & m=\mathrm{G}\mathrm{M},\mathrm{V}\mathrm{L},\mathrm{S}\mathrm{O},\mathrm{G}\mathrm{A},\hfill \\ \left(\mathrm{0,0,0}\right)\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e};\hfill \end{array}\right.(19)$

and for the Mns $(i\in \left\{\mathrm{M}\mathrm{n}\right\})$,

$\left(k_{im}^{\mathrm{v}},k_{im}^{\mathrm{d}},k_{im}^{\mathrm{f}}\right)=\left\{\begin{array}{cc}\left(\eta k_{\mathrm{F}}^{\mathrm{v}},\eta k_{\mathrm{F}}^{\mathrm{d}},\eta k_{\mathrm{F}}^{\mathrm{f}}\right)\hfill & i=\mathrm{M}\mathrm{n}‐m\left(m=\mathrm{I}\mathrm{P},\mathrm{T}\mathrm{A},\mathrm{B}\mathrm{F}\right),\hfill \\ \left(\eta k_{\mathrm{E}}^{\mathrm{v}},\eta k_{\mathrm{E}}^{\mathrm{d}},\eta k_{\mathrm{E}}^{\mathrm{f}}\right)\hfill & i=\mathrm{M}\mathrm{n}‐m\left(m=\mathrm{G}\mathrm{M},\mathrm{V}\mathrm{L},\mathrm{S}\mathrm{O},\mathrm{G}\mathrm{A}\right),\hfill \\ \left(\mathrm{0,0,0}\right)\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\hfill \end{array}\right.(20)$

We determined $k_{\mathrm{F}}^{\mathrm{v}},k_{\mathrm{F}}^{\mathrm{d}},k_{\mathrm{F}}^{\mathrm{f}},k_{\mathrm{E}}^{\mathrm{v}},k_{\mathrm{E}}^{\mathrm{d}},k_{\mathrm{E}}^{\mathrm{f}}$, and η through the optimization technique described in the next section.

4.2 Optimization to Determine Motor Control Parameters

We integrated the CPG and musculoskeletal models and determined the motor control parameters through an optimization to produce normal walking on the treadmill based on the measured kinematic data for cats (Prilutsky et al., 2016). Specifically, we determined the following 14 parameters; α_Mn-IP,PF-F, α_Mn-GM,PF-E, α_Mn-VL,PF-E, α_Mn-TA,PF-F, α_Mn-SO,PF-E, α_Mn-BF,PF-F, and α_Mn-GA,PF-E in Eq. 8 (the other a_ij’s for $i\in \left\{\mathrm{M}\mathrm{n}\right\}$ are 0) and $k_{\mathrm{F}}^{\mathrm{v}},k_{\mathrm{F}}^{\mathrm{d}},k_{\mathrm{F}}^{\mathrm{f}},k_{\mathrm{E}}^{\mathrm{v}},k_{\mathrm{E}}^{\mathrm{d}},k_{\mathrm{E}}^{\mathrm{f}}$, and η in Eqs. 18–20. We used the covariance matrix adaptation evolution strategy (CMA-ES) (Hansen et al., 2003) as the optimization method to minimize the following loss function:

$\epsilon ={\int }_{t_1}^{t_2}\sum _{i\in \left\{\mathrm{J}\right\}}{\left({\theta }_i-{\widehat{\theta }}_i\right)}^{2}dt(21)$

where $\left\{\mathrm{J}\right\}=\{\mathrm{H}\mathrm{i}\mathrm{p},\mathrm{K}\mathrm{n}\mathrm{e}\mathrm{e},\mathrm{A}\mathrm{n}\mathrm{k}\mathrm{l}\mathrm{e}\}$, θ_i and ${\widehat{\theta }}_i$ are the simulated and measured joint angles, and t₁ and t₂ were determined to evaluate the error for two gait cycles after steady walking was achieved.

Figure 5 shows the results obtained by the optimization. The model parameters determined by the optimization are shown in Appendix B. Figure 5A shows the membrane potentials of the RG, PF, and Mn neurons and the velocity, length, and force feedback of the flexor and extensor muscles. The rhythmic activities of the RG-F and RG-E neurons were transmitted to the Mn neurons through the PF neurons and the afferent feedback was used in normal walking. Figures 5B,C show the joint angles and muscle activations, respectively, compared with the data measured for cats (Prilutsky et al., 2016). The simulated joint angles and muscle activities have patterns similar to those for the measured data. These results verify that our model reproduced regular cat walking on a treadmill through the interactions between the neural system, the musculoskeletal system, and the environment by the optimization.

FIGURE 5

FIGURE 5. Simulation results for normal walking obtained by the optimization: (A) Membrane potentials of CPG neurons and afferent feedback from flexor and extensor muscles. (B) Joint angles and (C) muscle activations compared with measured data in cat [adapted from Prilutsky et al. (2016)]. R is the correlation coefficient and S is the cosine similarity. Liftoffs are represented by 0 and 100% in the gait cycle. Vertical lines indicate touchdowns.

5 Adaptive Response to Stepping Into a Hole During Treadmill Walking

To investigate the role of afferent feedback in adaptive locomotion, we applied a disturbance to our model and compared the responses with those measured for cats. Moreover, we clarified the adaptation mechanism based on dynamic systems theory methods.

5.1 Comparison of Simulation Results and Measured Data

Afferent feedback plays important roles in adaptive walking. To clarify these roles, a treadmill with a hole has been used in experiments with cats (Gorassini et al., 1994; Hiebert et al., 1994), where the cat body and forelimbs are supported and the cat walks on the treadmill only with the hindlimbs. When the foot of a hindlimb steps into the hole, it loses ground contact, and an appropriate response based on afferent feedback is required to continue walking. The investigation of the response to stepping into a hole highlights the roles of afferent feedback from the flexor and extensor muscles in generating adaptive walking. Hiebert et al. (1994) showed that when the foot of a hindlimb of a spinal cat entered a hole, the activities of the extensor muscles were shortened and the next onsets of the activities of the flexor and extensor muscles were advanced, which quickly lifted the foot out of the hole. We conducted a simulation of our model in the same situation as Hiebert et al. (1994) using the model parameters determined in the previous section to generate normal walking and compared the results with the measured data from animal experiments. In the simulation, after our model achieved steady walking on a treadmill, we emulated a hole by setting the ground reaction force to zero when the foot touched the treadmill belt until the foot went below the belt and returned back above it.

Figure 6 compares the simulation results when the foot entered a hole with those during normal walking. Specifically, Figure 6A shows the joint angles and stick diagrams. When the foot stepped into the hole, the hip, knee, and ankle joints were first extended. After that, these joints were flexed to lift the foot above the treadmill belt and then the foot landed on the belt. This behavior was consistent with that observed in the animal experiment. Figure 6B shows the activation of the IP and GA muscles representing the flexor and extensor muscles, respectively. The activation of the GA muscle was shortened and the onset of the next activation was advanced compared with those in normal walking. The onset of the next activation of the IP muscle was also advanced. These responses were consistent with those in the animal experiment. These results verify that our model reproduced adaptive responses in cats when a foot steps into a hole during walking, even without changing the model parameters, which were determined to reproduce normal walking.

FIGURE 6

FIGURE 6. Simulation results when on paw enters a hole compared with those during normal walking: (A) Joint angles and stick diagrams compared with measured data [adapted from Hiebert et al. (1994) for entering a hole and from Prilutsky et al. (2016) for normal walking]. These data are shown from the liftoff before entering a hole to the touchdown after entering the hole. R is the correlation coefficient of the results after entering the hole. The stick diagrams are shifted to the left by according to the distance traveled by the treadmill belt. (B) Activities of IP and GA muscles compared with measured data in spinal cat [adapted from Hiebert et al. (1994)]. Upper arrows indicate the activity duration, and lower ones indicate the intervals between the onsets of current and subsequent activities.

5.2 Investigating the Adaptation Mechanism

We investigated the adaptation mechanism via afferent feedback when the foot stepped into a hole based on the nullclines of the RG neurons in the phase plane. Figure 7 compares the stimulation results for the RG neurons between entering a hole and normal walking. Specifically, Figure 7A shows the time profiles of the membrane potentials of the RG-F and RG-E neurons and the velocity, length, and force feedback from the flexor and extensor muscles. Figures 7B,C show the trajectories of the state $\left(V_{\mathrm{R}\mathrm{G}‐\mathrm{F}},h_{\mathrm{R}\mathrm{G}‐\mathrm{F}}\right)$ of the RG-F neuron and the state $\left(V_{\mathrm{R}\mathrm{G}‐\mathrm{E}},h_{\mathrm{R}\mathrm{G}‐\mathrm{E}}\right)$ of the RG-E neuron, and the changes in the nullclines ${\widehat{N}}_{\mathrm{R}\mathrm{G}‐\mathrm{F}}^{\mathrm{V}}$ of the RG-F neuron and ${\widehat{N}}_{\mathrm{R}\mathrm{G}‐\mathrm{E}}^{\mathrm{V}}$ of the RG-E neuron (because ${\widehat{N}}_{\mathrm{R}\mathrm{G}‐\mathrm{F}}^{\mathrm{h}}$ and ${\widehat{N}}_{\mathrm{R}\mathrm{G}‐\mathrm{E}}^{\mathrm{h}}$ do not change, they are not shown), respectively, in their phase planes. Although a stable limit cycle appeared during normal walking, the trajectories diverged after the foot entered the hole. Therefore, just before the foot entered the hole, which corresponds to the touchdown in normal walking (“a” in Figure 7), the trajectories and nullclines were identical between entering a hole and normal walking.

FIGURE 7

FIGURE 7. Comparison of the response of RG neurons to a foot entering a hole with steady oscillation during normal walking. (A) Time profiles of membrane potentials and afferent feedback from flexor and extensor muscles. “a” indicates the time at which the foot enters the hole (touchdown in normal walking). “b”–“e” indicate the lapsed times of 0.1, 0.29, 0.43, and 0.56 s, respectively, after “a”. (B) Trajectories of $\left(V_{\mathrm{R}\mathrm{G}‐\mathrm{F}},h_{\mathrm{R}\mathrm{G}‐\mathrm{F}}\right)$ and $\left(V_{\mathrm{R}\mathrm{G}‐\mathrm{E}},h_{\mathrm{R}\mathrm{G}‐\mathrm{E}}\right)$ in phase planes. (C) Changes of nullclines ${\widehat{N}}_{\mathrm{R}\mathrm{G}‐\mathrm{F}}^{\mathrm{V}}$ and ${\stackrel{‐}{N}}_{\mathrm{R}\mathrm{G}‐\mathrm{E}}^{\mathrm{V}}$.

In the RG-F neuron, the local maximum of the nullcline decreased 0.1 s after entering the hole (“b” in Figure 7) due to the increased velocity and length feedback from the flexor muscles. The trajectory was attracted to the nullcline while following the slow dynamics. As a result, V_RG-F increased slightly. In the RG-E neuron, the nullcline remained almost unchanged from that in “a” due to the loss of the force feedback from the extensor muscles. In other words, the nullcline moved to the left compared to that in normal walking. The trajectory was attracted to the nullcline while following the slow dynamics. As a result, V_RG-E decreased slightly.

In the RG-F neuron, the local maximum of the nullcline further decreased 0.29 s after the foot entered the hole (“c” in Figure 7) due to the increase in the length and force feedback from the flexor muscles. The state entered the fast dynamics by passing over the local maximum. As a result, V_RG-F rapidly underwent a large increase to initiate the next activity. In the RG-E neuron, the nullcline moved significantly to the left due to the increased inhibition from the RG-F neuron, which made the state enter the fast dynamics. As a result, V_RG-E decreased strongly and the activity was rapidly terminated.

As a result of these changes, the onset of the next activity of the RG-F neuron was advanced and the current activity of the RG-E neuron was shortened. The state $\left(V_{\mathrm{R}\mathrm{G}‐\mathrm{F}},h_{\mathrm{R}\mathrm{G}‐\mathrm{F}}\right)$ of the RG-F neuron and the state $\left(V_{\mathrm{R}\mathrm{G}‐\mathrm{E}},h_{\mathrm{R}\mathrm{G}‐\mathrm{E}}\right)$ of the RG-E neuron followed the slow dynamics in “d” in Figure 7 and the trajectories and nullclines in “e” in Figure 7 almost returned to those of “a” in Figure 7 through the fast dynamics to continue normal walking. These results clarify the mechanims how our model produced adaptive responses during stepping into a hole through the interactions between the neural system, the musculoskeletal system, and the environment using dynamic systems theory based on nullclines.

6 Discussion

6.1 Adaptive Responses to Entering a Hole

To continue walking after one of its feet steps into a hole, the cat must quickly pull the foot out of the hole and step onto the treadmill belt again. To provide this behavior the CPG should modify the motor commands based on afferent feedback. Similar to normal walking, the extensor muscles start their activity before touchdown to support the body weight. However, when a foot enters a hole without touchdown, the muscle activity excessively extends the hip, knee, and ankle joints. Therefore, for stable walking it is necessary to quickly terminate extension and start flexion. It has been suggested that unloading of ankle extensor muscles and activation of the Ia and II sensory afferents from flexor muscles play an important role in the transition from stance to swing during normal walking in cats (Duysens and Pearson, 1980; Pearson et al., 1992; Hiebert et al., 1996). Studies of locomotion in spinal cats when holes were present showed that the unloading of the extensor muscle due to the lack of touchdown caused early termination of the extensor activity (Hiebert et al., 1994). In our model, this phase transition was provided by the velocity and length dependent (Ia, II) as well as the force dependent (Ib) feedback. Our simulations have shown that the early termination of extensor activity and early initiation of flexor activity can result from unloading of the extensor muscles (“b” in Figure 7C) and increasing the activity of Ia and II afferent feedback from the flexor muscles of knee and ankle (“c” in Figure 7C) when the foot entering into a hole, which is consistent with the animal experiments. It is important to mention that, although the parameters for afferent feedback were determined during optimization performed to reproduce the normal walking, the afferent feedback incorporated in the model allowed the model to reproduce adaptive changes of locomotor characteristics during stepping into a hole. Furthermore, the mechanism responsible for this adaptive response was explained by the transition from the slow to fast dynamics in the RG neurons induced by the changes in the nullclines due to afferent feedback (Figure 7).

The adaptive locomotor behavior during stepping into a hole can be also interpreted based on the phase-dependent response of the CPG activity. Stimulation of extensor afferents during the active phase in fictive locomotion of decerebrate cats was shown to increase and prolong the extensor activity (Duysens, 1977), and our isolated CPG model was able to reproduce the phase-dependent response (Figure 4B). In the model, during in normal walking, Ib feedback from extensors increased during the stance phase due to the ground reaction force, which increased and prolonged the activities of the extensor half center of the CPG (Figure 5B). However, when the foot entered the hole, the force feedback from the extensor muscles was reduced due to the lack of the ground reaction force and the activity of the extensor half center was reduced (Figures 6B, 7A). Similarly, stimulation of flexors afferents during the silent phase in fictive locomotion in decerebrate cats inhibited the extensor activity and initiated transition to flexion (Schomburg et al., 1998), and our isolated CPG model also reproduced these phase-dependent responses (Figure 4A). The flexor muscles in our model were stretched when the foot entered into a hole and afferent feedback from the flexor side increased (Figure 6B), which inhibited the activities of extensor muscles and shortened the stance phase leading to an advanced onset of the swing phase. As a result, the extension of the joints stopped and they were strongly flexed to pull the foot out of the hole.

6.2 Limitations and Future Works

In this study, we focused only on simulation of walking with a single hole as a source of locomotor disturbance. This allowed us to investigate the roles of afferent feedback from both flexor and extensor muscles because stepping into a hole caused unloading of extensor muscles and stretching of flexor muscles. We have not considered other walking disturbances tested experimentally, such as the presence of obstacles, cutaneous stimulation, and mechanical blockage of limb flexion (Forssberg et al., 1977; Lam and Pearson, 2001; Drew et al., 2002). In particular, while we focused on a hole without a bottom, a foot can have touchdown in a hole with a bottom depending on the depth, which changes the afferent feedback and phase transition. In the future, we would like to simulate walking with these other disturbances to provide mechanistic understanding of adaptation observed during locomotion. In addition, our model would be useful to investigate the responses when afferent feedback does not work well due to the paralysis and provides wrong information.

Because our present model used only one pair of the PF neurons (a PF-F neuron for the flexor side and a PF-E neuron for the extensor side), the motor commands activated the flexor and extensor muscles in strict alternation. However, in the biological system, some muscles show activity in both swing and stance phases and some muscles are active during only part of one phase (Prilutsky et al., 2016). In addition, muscle synergy analysis suggested that muscle activities during animal walking can be explained by the combination of a few basic patterns (Dominici et al., 2011). Several previous CPG models have been proposed to produce such multiple patterns by increasing the complexity of the PF network (Markin et al., 2016). It has been also suggested that some afferent feedback is projected to the contralateral side of the CPG (McCrea, 2001). We will include these components in a more elaborated model in the future.

Also in this work, we only modeled one hindlimb, which did not allow us to consider interlimb coordination. The latter is also known to play an important role in adaptive walking (Forssberg et al., 1980). Particularly, when a foot of one limb enters a hole, the activity of the extensor muscle of the contralateral limb has been shown to increase to compensate for the weight support (Hiebert et al., 1994). In addition, animals change their gait patterns, such as walking, trotting, and galloping, depending on locomotor speed by controlling the interlimb coordination. A neural model simulating coordination between RG circuits controlling each limb has been proposed (Danner et al., 2016, 2017). These interactions will be included in the future extended model simulating locomotion with two forelimbs and two hindlimbs.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author Contributions

SA developed the study design in consultation with the other co-authors. YK performed simulation experiments and analyzed the data in consultation with SA, SF, SD, SM, JA, and IR. YK and SA wrote the manuscript and SD, SM, JA, and IR revised it. All the authors reviewed and approved it.

Funding

This study was supported in part by JSPS KAKENHI Grant Numbers JP19KK0377, JP20H00229, and JP20KK0226; JST FOREST Program Grant Number JPMJFR2021; NIH grants R01NS100928, R01NS110550, R01NS112304, and R01NS115900; and NSF grant 2113069.

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

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 (23), 6947–6967. doi:10.1113/jp272787

PubMed Abstract | CrossRef Full Text | Google Scholar

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, 1–25. doi:10.7554/eLife.31050

CrossRef Full Text | Google Scholar

Demir, S. S., Butera, R. J., DeFranceschi, A. A., Clark, J. W., and Byrne, J. H. (1997). Phase Sensitivity and Entrainment in a Modeled Bursting Neuron. Biophysical J. 72 (2), 579–594. doi:10.1016/s0006-3495(97)78697-1

CrossRef Full Text | Google Scholar

Dominici, N., Ivanenko, Y. P., Cappellini, G., d’Avella, A., Mondì, V., Cicchese, M., et al. (2011). Locomotor Primitives in Newborn Babies and Their Development. Science 334 (6058), 997–999. doi:10.1126/science.1210617

PubMed Abstract | CrossRef Full Text | Google Scholar

Drew, T., Jiang, W., and Widajewicz, W. (2002). Contributions of the Motor Cortex to the Control of the Hindlimbs during Locomotion in the Cat. Brain Res. Brain Res. Rev. 40 (1-3), 178–191. doi:10.1016/s0165-0173(02)00200-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Duysens, J., and Pearson, K. G. (1980). Inhibition of Flexor Burst Generation by Loading Ankle Extensor Muscles in Walking Cats. Brain Res. 187 (2), 321–332. doi:10.1016/0006-8993(80)90206-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Duysens, J. (1977). Reflex Control of Locomotion as Revealed by Stimulation of Cutaneous Afferents in Spontaneously Walking Premammillary Cats. J. Neurophysiol. 40 (4), 737–751. doi:10.1152/jn.1977.40.4.737

PubMed Abstract | CrossRef Full Text | Google Scholar

Ekeberg, Ö., and Pearson, K. (2005). Computer Simulation of Stepping in the Hind Legs of the Cat: An Examination of Mechanisms Regulating the Stance-To-Swing Transition. J. Neurophysiol. 94 (6), 4256–4268. doi:10.1152/jn.00065.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, B. (1994). Reduction of Conductance-Based Models with Slow Synapses to Neural Nets. Neural Comput. 6 (4), 679–695. doi:10.1162/neco.1994.6.4.679

CrossRef Full Text | Google Scholar

Forssberg, H., Grillner, S., and Rossignol, S. (1977). Phasic Gain Control of Reflexes from the Dorsum of the Paw during Spinal Locomotion. Brain Res. 132 (1), 121–139. doi:10.1016/0006-8993(77)90710-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Forssberg, H., Grillner, S., Halbertsma, J., and Rossignol, S. (1980). The Locomotion of the Low Spinal Cat. Ii. Interlimb Coordination. Acta Physiol. Scand. 108 (3), 283–295. doi:10.1111/j.1748-1716.1980.tb06534.x

PubMed Abstract | CrossRef Full Text | Google Scholar

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 (3), 1591–1605. doi:10.1152/jn.01028.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujiki, S., Aoi, S., Tsuchiya, K., Danner, S. M., Rybak, I. A., and Yanagihara, D. (2019). Phase-dependent Response to Afferent Stimulation during Fictive Locomotion: A Computational Modeling Study. Front. Neurosci. 13, 1288. doi:10.3389/fnins.2019.01288

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorassini, M. A., Prochazka, A., Hiebert, G. W., and Gauthier, M. J. (1994). Corrective Responses to Loss of Ground Support during Walking. I. Intact Cats. J. Neurophysiol. 71 (2), 603–610. doi:10.1152/jn.1994.71.2.603

PubMed Abstract | CrossRef Full Text | Google Scholar

Grillner, S. (1981). “Control of Locomotion in Bipeds, Tetrapods, and Fish,” in Handbook of Physiology, Sect. 1. The Nervous System: Motor Control. Editor V. B. Brooks (Bethesda, MD: American Physiological Society), Vol. II, Pt. 2, 1179–1236. doi:10.1002/cphy.cp010226

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, N., Müller, S. D., and Koumoutsakos, P. (2003). Reducing the Time Complexity of the Derandomized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES). Evol. Comput. 11 (1), 1–18. doi:10.1162/106365603321828970

PubMed Abstract | CrossRef Full Text | Google Scholar

Hiebert, G. W., Gorassini, M. A., Jiang, W., Prochazka, A., and Pearson, K. G. (1994). Corrective Responses to Loss of Ground Support during Walking. ii. Comparison of Intact and Chronic Spinal Cats. J. Neurophysiol. 71 (2), 611–622. doi:10.1152/jn.1994.71.2.611

PubMed Abstract | CrossRef Full Text | Google Scholar

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 (3), 1126–1137. doi:10.1152/jn.1996.75.3.1126

PubMed Abstract | CrossRef Full Text | Google Scholar

Lam, T., and Pearson, K. G. (2001). Proprioceptive Modulation of Hip Flexor Activity during the Swing Phase of Locomotion in Decerebrate Cats. J. Neurophysiol. 86 (3), 1321–1332. doi:10.1152/jn.2001.86.3.1321

PubMed Abstract | CrossRef Full Text | Google Scholar

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. Acad. Sci. 1198, 21–34. doi:10.1111/j.1749-6632.2010.05435.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Markin, S. N., Klishko, A. N., Shevtsova, N. A., Lemay, M. A., Prilutsky, B. I., and Rybak, I. A. (2016). “A Neuromechanical Model of Spinal Control of Locomotion,” in Neuromechanical Modeling of Posture and Locomotion. Editors B. I. Prilutsky, and D. H. Edwards (New York: Springer Science), 21–65. doi:10.1007/978-1-4939-3267-2_2

CrossRef Full Text | Google Scholar

McCrea, D. A. (2001). Spinal Circuitry of Sensorimotor Control of Locomotion. J. Physiol. 533 (1), 41–50. doi:10.1111/j.1469-7793.2001.0041b.x

PubMed Abstract | CrossRef Full Text | Google Scholar

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 (5), e1004270–40. doi:10.1371/journal.pcbi.1004270

PubMed Abstract | CrossRef Full Text | Google Scholar

Orlovsky, G., Deliagina, T. G., and Grillner, S. (1999). Neuronal Control of Locomotion: From Mollusc to Man. Oxford, United Kingdom: Oxford University Press.

Google Scholar

Pearson, K. G., Ramirez, J. M., and Jiang, W. (1992). Entrainment of the Locomotor Rhythm by Group Ib Afferents from Ankle Extensor Muscles in Spinal Cats. Exp. Brain Res. 90 (3), 557–566. doi:10.1007/BF00230939

PubMed Abstract | CrossRef Full Text | Google Scholar

Perreault, M. C., 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

PubMed Abstract | CrossRef Full Text | Google Scholar

Prilutsky, B. I., Klishko, A. N., Weber, D. J., and Lemay, M. A. (2016). “Computing Motion Dependent Afferent Activity during Cat Locomotion Using a Forward Dynamics Musculoskeletal Model,” in Neuromechanical Modeling of Posture and Locomotion. Editors B. I. Prilutsky, and D. H. Edwards (New York: Springer Science), 273–307. doi:10.1007/978-1-4939-3267-2_10

CrossRef Full Text | Google Scholar

Rossignol, S., Dubuc, R., and Gossard, J.-P. (2006). Dynamic Sensorimotor Interactions in Locomotion. Physiol. Rev. 86 (1), 89–154. doi:10.1152/physrev.00028.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Rossignol, S. (1996). “Neural Control of Stereotypic Limb Movements,” in Handbook of Physiology, Sect. 12. Exercise: Regulation and Integration of Multiple Systems. Editors L. B. Rowell, and J. T. Sheperd (Oxford: American Physiological Society), 173–216. doi:10.1002/cphy.cp120105

CrossRef Full Text | Google Scholar

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 (2), 617–639. doi:10.1113/jphysiol.2006.118703

PubMed Abstract | CrossRef Full Text | Google Scholar

Rybak, I. A., Stecina, K., Shevtsova, N. A., and McCrea, D. A. (2006b). Modelling Spinal Circuitry Involved in Locomotor Pattern Generation: Insights from the Effects of Afferent Stimulation. J. Physiol. 577 (2), 641–658. doi:10.1113/jphysiol.2006.118711

PubMed Abstract | CrossRef Full Text | Google Scholar

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 (3), 339–350. doi:10.1007/s002210050522

PubMed Abstract | CrossRef Full Text | Google Scholar

Spardy, L. E., Markin, S. N., Shevtsova, N. A., Prilutsky, B. I., Rybak, I. A., and Rubin, J. E. (2011). A Dynamical Systems Analysis of Afferent Control in a Neuromechanical Model of Locomotion: I. Rhythm Generation. J. Neural Eng. 8 (6), 065003–065031. doi:10.1088/1741-2560/8/6/065003

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Yakovenko, S., Gritsenko, V., and Prochazka, A. (2004). Contribution of Stretch Reflexes to Locomotor Control: A Modeling Study. Biol. Cybernetics 90 (2), 146–155. doi:10.1007/s00422-003-0449-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix

A Parameters for CPG Model

Based on Molkov et al. (2015), Danner et al. (2016), Danner et al. (2017), Markin et al. (2010), and Fujiki et al. (2019), we determined the parameters for the synaptic connections α_ij, β_ij, and γ_i $(i\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{I}\mathrm{n}\right\},\left\{\mathrm{P}\mathrm{F}\right\},j\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{I}\mathrm{n}\right\})$ and those for the currents ${\widehat{g}}_{\mathrm{L}\mathrm{e}\mathrm{a}\mathrm{k}}^i$, ${\widehat{g}}_{\mathrm{N}\mathrm{a}\mathrm{P}}^i$, $E_{\text{SynI}}^i$, and $E_{\text{Leak}}^i$ $(i\in \left\{\mathrm{R}\mathrm{G}\right\},\left\{\mathrm{I}\mathrm{n}\right\},\left\{\mathrm{P}\mathrm{F}\right\},\left\{\mathrm{M}\mathrm{n}\right\})$ as shown in Tables 1, 2, respectively. We used the following parameter values: ${\widehat{g}}_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{E}}=10\mathrm{n}\mathrm{S}$, ${\widehat{g}}_{\mathrm{S}\mathrm{y}\mathrm{n}\mathrm{I}}=10\mathrm{n}\mathrm{S}$, E_SynE = − 10 mV, E_Na = 55 mV, C = 20 pF, V_th = − 50 mV, V_max = 0 mV, and d = 1.0.

TABLE 1

TABLE 1. Parameters for synaptic connections.

TABLE 2

TABLE 2. Parameters for currents.

B Parameters Determined by Optimization to Generate Normal Walking

Through the optimization needed to generate normal walking, the parameters for the synaptic connections α_ij $(i\in \left\{\mathrm{M}\mathrm{n}\right\},j\in \left\{\mathrm{P}\mathrm{F}\right\})$ were determined as shown in Table 3 and the weight coefficients for afferent feedback were determined as follows: $k_{\mathrm{F}}^{\mathrm{v}}=0.001$, $k_{\mathrm{F}}^{\mathrm{d}}=0.15$, $k_{\mathrm{F}}^{\mathrm{f}}=0.077$, $k_{\mathrm{E}}^{\mathrm{v}}=0.00039$, $k_{\mathrm{E}}^{\mathrm{d}}=0.058$, $k_{\mathrm{E}}^{\mathrm{f}}=0.2$, and η = 0.27.

TABLE 3

TABLE 3. Optimized parameters for synaptic connections.