Front. Neurosci. Frontiers in Neuroscience Front. Neurosci. 1662-453X Frontiers Media S.A. 10.3389/fnins.2020.00017 Neuroscience Original Research Contribution of Phase Resetting to Adaptive Rhythm Control in Human Walking Based on the Phase Response Curves of a Neuromusculoskeletal Model Tamura Daiki 1 Aoi Shinya 1 * Funato Tetsuro 2 Fujiki Soichiro 3 Senda Kei 1 Tsuchiya Kazuo 1 1Department of Aeronautics and Astronautics, Graduate School of Engineering, Kyoto University, Kyoto, Japan 2Department of Mechanical Engineering and Intelligent Systems, Graduate School of Informatics and Engineering, The University of Electro-Communications, Tokyo, Japan 3Department of Physiology and Biological Information, School of Medicine, Dokkyo Medical University, Tochigi, Japan

Edited by: Tomohiko Takei, Kyoto University, Japan

Reviewed by: Srinivasa Chakravarthy, Indian Institute of Technology Madras, India; Auke Ijspeert, École Polytechnique Fédérale de Lausanne, Switzerland

 *Correspondence: Shinya Aoi shinya_aoi@kuaero.kyoto-u.ac.jp

This article was submitted to Decision Neuroscience, a section of the journal Frontiers in Neuroscience

 05 02 2020 2020 14 17 03 04 2019 09 01 2020 Copyright © 2020 Tamura, Aoi, Funato, Fujiki, Senda and Tsuchiya. 2020 Tamura, Aoi, Funato, Fujiki, Senda and Tsuchiya

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Humans walk adaptively in varying environments by manipulating their complicated and redundant musculoskeletal system. Although the central pattern generators in the spinal cord are largely responsible for adaptive walking through sensory-motor coordination, it remains unclear what neural mechanisms determine walking adaptability. It has been reported that locomotor rhythm and phase are regulated by the production of phase shift and rhythm resetting (phase resetting) for periodic motor commands in response to sensory feedback and perturbation. While the phase resetting has been suggested to make a large contribution to adaptive walking, it has only been investigated based on fictive locomotion in decerebrate cats, and thus it remains unclear if human motor control has such a rhythm regulation mechanism during walking. In our previous work, we incorporated a phase resetting mechanism into a motor control model and demonstrated that it improves the stability and robustness of walking through forward dynamic simulations of a human musculoskeletal model. However, this did not necessarily verify that phase resetting plays a role in human motor control. In our other previous work, we used kinematic measurements of human walking to identify the phase response curve (PRC), which explains phase-dependent responses of a limit cycle oscillator to a perturbation. This revealed how human walking rhythm is regulated by perturbations. In this study, we integrated these two approaches using a physical model and identification of the PRC to examine the hypothesis that phase resetting plays a role in the control of walking rhythm in humans. More specifically, we calculated the PRC using our neuromusculoskeletal model in the same way as our previous human experiment. In particular, we compared the PRCs calculated from two different models with and without phase resetting while referring to the PRC for humans. As a result, although the PRC for the model without phase resetting did not show any characteristic shape, the PRC for the model with phase resetting showed a characteristic phase-dependent shape with trends similar to those of the PRC for humans. These results support our hypothesis and will improve our understanding of adaptive rhythm control in human walking.

 human walking phase resetting phase response curve central pattern generator muscle synergy neuromusculoskeletal model JP15KT0015 JP26120006 Japan Society for the Promotion of Science10.13039/501100001691 1. Introduction

Humans walk adaptively in varying environments by the skillful control of their complicated and redundant musculoskeletal system. Although many studies have investigated the underlying mechanism for adaptive walking, it remains largely unclear what neural mechanisms determine the walking adaptability.

Because human walking is rhythmic, elucidating the rhythm control strategy is crucial. The central pattern generators (CPGs) in the spinal cord are largely responsible for adaptive rhythm control through sensory-motor coordination (Orlovsky et al., 1999). In particular, it has been reported that locomotor rhythm and phase are regulated by producing phase shift and rhythm resetting (phase resetting) for periodic motor commands in response to sensory feedback and perturbation (Duysens, 1977; Conway et al., 1987; Guertin et al., 1995; Schomburg et al., 1998; Lafreniere-Roula and McCrea, 2005; Rybak et al., 2006a; Frigon and Gossard, 2010). However, such phase resetting behavior has been investigated only with electromyographic and electroneurographic data measured during fictive locomotion in decerebrate cats, and thus it is unclear if human motor control has such a rhythm regulation mechanism during walking. From a modeling approach on the basis of the hypothesis that phase resetting works for the control of walking rhythm in humans, the phase resetting mechanism has been introduced in motor control models of human walking. Although the models demonstrated that it improves stability and robustness of walking through forward dynamic simulations of human musculoskeletal models (Yamasaki et al., 2003a,b; Nomura et al., 2009; Aoi et al., 2010; Aoi and Funato, 2016), they did not necessarily verify whether the hypothesis is true.

To investigate rhythm regulation mechanisms in biological and natural phenomena, researchers have applied the phase response curve (PRC) in the phase reduction theory, which explains how the phase of a limit cycle oscillator shifts by a perturbation at an arbitrary phase (Kuramoto, 1984; Winfree, 2001). In our previous work (Funato et al., 2016), we assumed human walking as a limit cycle oscillator and identified the PRC from kinematic measurements by changing the belt speed of a treadmill during human walking, which clarified how human walking rhythm is regulated by perturbations. In this study, to examine the hypothesis, we integrated two previous different approaches that used a physical model and identification of the PRC. More specifically, we performed forward dynamic simulations with our previous neuromusculoskeletal model (Aoi et al., 2010) to walk on a treadmill and disturbed the belt speed at arbitrary phases in the same way as our previous experiments with humans (Funato et al., 2016). In particular, we obtained the PRC for two different cases with and without phase resetting in our motor control model and compared the results with the measured PRC in humans. Based on these results, we discuss the contribution of phase resetting to adaptive rhythm control in human walking.

 2. Methods 2.1. Model

In this study, we used the same neuromusculoskeletal model that we developed in our previous work (Aoi et al., 2010). We briefly explain the model below.

2.1.1. Musculoskeletal Model

Our musculoskeletal model is two-dimensional (Figure 1A), and the physical parameters were determined from data obtained from measurement of human walking (Davy and Audu, 1987; Winter, 2004). The skeletal part of our model has seven rigid links: trunk (head, arms, and torso) and thigh, shank, and foot of each leg, and has nine degrees of freedom: hip, knee, and ankle joint angles of each leg and horizontal and vertical translations and rotation of the trunk. Each joint has a linear viscous element, and the knee and ankle joints are subject to large linear elastic and damping torques when these joint angles exceed their limits. We used four contact points on each sole to receive reaction forces from the treadmill belt (toe, heel, and 4.0 cm inside from the toe and from the heel). The reaction force is modeled by a linear spring and damper system for each horizontal and vertical direction. Our model contains nine principal muscles to achieve the necessary motions in each leg. Six muscles produce uniarticular motion: hip flexion [iliopsoas (IL)], hip extension [gluteus maximus (GM)], knee extension [vastus (VA)], knee flexion [biceps femoris short head (BFS)], ankle flexion [tibialis anterior (TA)], and ankle extension [soleus (SO)]. Three muscles produce biarticular motion: hip flexion and knee extension [rectus femoris (RF)], hip extension and knee flexion [biceps femoris long head (BFL)], and knee flexion and ankle extension [gastrocnemius (GC)]. The muscle model consists of contractile and passive elements. The contractile part depends on force-length and force-velocity relationships and the muscle activation, which is determined through a low-pass filtering of the motor command um (m = IL, GM, VA, BFS, TA, SO, RF, BFL, and GC) from the motor control model. The equations of motion in this model were derived using Lagrangian mechanics and solved using the fourth-order Runge-Kutta method with time steps of 2 × 10−7 s for the forward dynamic simulation.

Neuromusculoskeletal model: (A) musculoskeletal model walking on a treadmill and (B) motor command composed of the linear combination of five rectangular pulses based on the muscle synergy hypothesis, and identification of the muscles activated by each pulse. Each sole has four contact points (two for the toe part and the others for the heel part) to receive reaction forces from the treadmill belt through linear spring and damper systems for each point.

 2.1.2. Motor Control Model

Our motor control model consists of a hypothetical two-layered CPG model at the spinal cord level, which incorporates phase resetting, and a movement regulation model at the brainstem and cerebellum levels. The CPGs in the spinal cord have been suggested to consist of hierarchical networks that include rhythm generator (RG) and pattern formation (PF) networks (Burke et al., 2001; Lafreniere-Roula and McCrea, 2005; Rybak et al., 2006a,b). The RG network generates the basic rhythm and alters it by producing phase shifts and rhythm resetting in response to sensory feedback, while the PF network shapes the rhythm into spatiotemporal patterns of motoneuron activities. For the RG model, we used two simple oscillators whose phase is ϕi (0 ≤ ϕi < 2π, i = right, left) to produce the basic rhythm of the corresponding leg and incorporated phase resetting as explained below. For the PF model, we determined motor commands necessary to produce periodic leg movements in accordance with the oscillator phase based on the muscle synergy hypothesis, which suggests that the linear combination of five basic signals produces a large portion of the motor commands for human locomotion (Ivanenko et al., 2006). More specifically, we used five rectangular pulses pi(ϕ) (i = 1, …, 5) for each leg (Figure 1B), which are given by:

pi(ϕ)={1Φi<ϕ≤Φi+ΔΦi0otherwisei=1,…,5

where Φi and ΔΦi (i = 1, …, 5) are the onset phase and duration, respectively, of the rectangular pulses, and we omitted the suffix of ϕ. We determined the muscle synergy-based motor command umSyn by:

umSyn=∑i=15wm,iΛipi(ϕ)

where wm, i (i = 1, …, 5) is the weighting coefficient of five rectangular pulses (wm, i ≥ 0) and Λi (i = 1, …, 5) is the tuning parameter of the amplitude of the rectangular pulses for different belt speeds.

To emulate the phase shift and rhythm resetting behavior, we incorporated the phase resetting mechanism in the RG model. More specifically, we reset the oscillator phase to a nominal value based on foot contact information by using the following phase dynamics:

ϕ∙i=ω-Kϕ sin(Δϕi-π)-(ϕi-ϕFC)δ(t-tiFC-τFC)

where,

Δϕi={ϕright-ϕlefti=rightϕleft-ϕrighti=left

ω is the basic frequency, Kϕ is the gain parameter, tiFC is the foot-contact time, τFC (= 50 ms) is the transmission delay in receiving the foot-contact information, ϕFC is the phase value to be reset at the foot contact, and δ(·) is the Dirac delta function. The second term on the right-hand side maintains interlimb coordination so that the legs move in antiphase. The third term of the right-hand side corresponds to phase resetting, which resets the oscillator phase ϕi to ϕFC to modulate the timing of the muscle synergy-based motor command based on the foot contact information. The second and third terms regulate the gait frequency and contribute to the generation of stable limit cycle for walking.

In addition to the CPG model at the spinal cord level, we used a movement regulation model at the brainstem and cerebellum levels based on somatosensory information, where only two crucial factors were incorporated for simplicity: maintenance of an upright posture and the desired forward speed. For the maintenance of an upright posture, a simple feedback control regulates the balance of the trunk pitch to prevent it from falling over using antagonistic uniarticular muscles in the hip of the standing leg.

pmTrunk={-κm(θ-θ^)-σmθ∙in stance phase0otherwise

where θ is the trunk pitch angle, θ∙ is the angular rate, θ^ is the reference angle, and κm and σm are the gain parameters (κm = σm = 0 when m ≠ IL or GM). For maintenance of the speed, a simple feedback control is used to increase the ankle push-off when the speed is lower than desired and suppress the pushing force in the opposite case by antagonistic uniarticular muscles in the ankle of the standing leg.

pmSpeed={-λm(v-v^)in stance phase0otherwise

where v is the forward speed, v^ is the target forward speed, and λm is the gain parameter (λm = 0 when m ≠ TA or SO). Because these regulations operate at the brainstem and cerebellar levels, the command signals are delayed and the motor command umReg is given by:

umReg(t)=pmTrunk(t-τReg)+pmSpeed(t-τReg)

where τReg (= 80 ms) is the delay in receiving transmissions of somatosensory information at the brainstem and cerebellar levels and sending the motor command to the spinal cord level.

The motor command um is given by the summation of the muscle synergy-based motor command umSyn and the motor command by the movement regulation umReg.

um=umSyn+umReg 2.1.3. Model Parameters

While the model in our previous work (Aoi et al., 2010) walked over the ground, the model in this study walked on a treadmill, as explained below. Therefore, we slightly modified the values of the motor control parameters so that the model achieved steady walking on the treadmill whose belt speed was 1.3 m/s as follows: the onset phase and duration of rectangular pulses were Φ1 = 6.12 rad, Φ2 = 1.48 rad, Φ3 = 2.56 rad, Φ4 = 3.51 rad, Φ5 = 5.38 rad, ΔΦ1 = 0.70 rad, ΔΦ2 = 0.90 rad, ΔΦ3 = 0.90 rad, ΔΦ4 = 1.07 rad, and ΔΦ5 = 0.96 rad, where we set ϕ = 0 rad at foot contact; the amplitudes and weighting coefficients of the rectangular pulses were Λi = 1.0 (i = 1, …, 5), wVA,1 = 0.42, wTA,1 = 0.35, wSO,2 = 1.26, wGC,2 = 0.87, wIL,3 = 1.02, wBFS,3 = 1.09, wRF,3 = 0.10, wVA,4 = 0.17, wTA,4 = 0.21, wGM,5 = 0.61, wBFS,5 = 0.20, wBFL,5 = 0.20, and the other wm,i = 0; the parameters for the oscillator phase dynamics were ω = 2π/1.0 rad/s, Kϕ = 1.7, and ϕFC = 0.36 rad; and the parameters for the movement regulation were κIL = −1.0, κGM = 2.0, σIL = −0.20, σGM = 0.40, λTA = −0.20, λSO = 0.12, θ^=-0.012 rad, and v^=0.1 m/s.

For different belt speeds, we changed Φ2, ω, Λi (i = 1, …, 5), and ϕFC in a similar way to Aoi et al. (2019) as follows: Φ2 = 1.46 rad, ω = 2π/0.9 rad/s, Λ1 = 1.04, Λ2 = 1.14, Λ3 = 1.10, Λ4 = 1.03, Λ5 = 1.18, and ϕFC = 0.48 rad when the belt speed was increased by 0.02 m/s, and Φ2 = 1.50 rad, ω = 2π/1.1 rad/s, Λ1 = 0.96, Λ2 = 0.90, Λ3 = 0.90, Λ4 = 0.98, Λ5 = 0.82, and ϕFC = 0.04 rad when the belt speed was decreased by 0.02 m/s.

 2.2. Phase Response Curve

In the phase reduction theory (Kuramoto, 1984; Winfree, 2001), for a limit cycle oscillator whose period is τ and closed orbit is C on the phase space (Figure 2), we can define ψ on C, which follows the dynamics:

ψ∙=2πτ

To apply the phase dynamics to the neighborhood of the limit cycle, we assume that the point P on C and the point Q close to C have the same phase when they converge to the same point on C for t → ∞. The surface (curve) with the same phase (ψ = ψ0= const.) is called an isochron.

Limit cycle orbit C and isochron. Point P on C and point Q close to C converge to the same point on C for t → ∞ and are included in the same isochron. Poincaré section S, which determines the cycles, generally mismatches with any of the isochrons.

When a perturbation I(t) is added to the limit cycle oscillator, ψ follows the dynamics:

ψ∙=2πτ+Z(ψ)I(t)

where Z(ψ) is the PRC and explains the phase-dependent rhythm change due to the perturbation. We determine cycles using Poincaré section S, as shown in Figure 2. We assume that the trajectory converges to C before I(t) is added. We define t = 0 for the time at the last intersection of C with S before I(t) is added and ψ(0) = 0. We also define t = tn (n = 1, 2, … ) for the time at the nth intersection of the disturbed trajectory with S after I(t) is added. The integration of (9) from 0 to tn gives:

∫0tn(ψ∙-2πτ)dt=∫0tnZ(ψ)I(t)dt

The Poincaré section S generally mismatches with the isochron of ψ = 0, as shown in Figure 2, which induces the difference of ψ between the Poincaré section and isochron and thus ∫0tnψ∙dt≠2nπ (Imai and Aoyagi, 2016). However, because the disturbed trajectory approaches C as t → ∞, ∫0tnψ∙dt=2nπ approximately for sufficiently large n. This gives:

∫0tnZ(ψ)I(t)dt=2πnτ-tnτ

The right-hand side can be obtained from the phase shift by the perturbation, as shown in Figure 3. For an impulsive perturbation at t = s (0 ≤ s < τ), which is given by I(t) = μδ(t − s) when μ is constant, (11) becomes:

μ∫0tnZ(ψ)δ(t-s)dt=2πnτ-tnτ

This gives,

Z(ψ(s))=2πμnτ-tnτ

In this study, we calculated the PRC from (13) using our neuromusculoskeletal model, where we used the foot contact condition of the right leg (any of four contact points of the right foot is below the treadmill belt) for the Poincaré section S. In particular, our previous work (Funato et al., 2016) obtained the PRCs from human walking measurements by accelerating or decelerating the belt speed of a treadmill independently. To compare the simulation results with the human measurements, our model also walked on a treadmill and we obtained the PRCs for the acceleration and deceleration perturbations in the belt speed separately. More specifically, after the model achieved steady walking on a treadmill, we increased or decreased the belt speed by 0.1 m/s for 0.001 s once per trial (μ=2π±0.1·0.001ντ), where ν is the belt speed, and n = 50 so that the model achieved steady walking after being disturbed. We performed 100 trials by changing the perturbation phase to obtain the PRC. Furthermore, we used the models with and without phase resetting in the motor control model and compared the PRCs calculated from these models by referring to the PRCs obtained from measurements of human walking.

Phase shift by perturbation of a limit cycle oscillator at t = s. In this case, a positive peak condition is used for the Poincaré section.

 3. Results

By modifying the motor control parameters from those in our previous work (Aoi et al., 2010), both the models with and without phase resetting achieved steady walking on a treadmill whose belt speed was 1.3 m/s. The locomotor behavior, especially the joint kinematics and muscle activities, were almost identical to those in our previous work (Aoi et al., 2010) except for the difference between walking over ground and on a treadmill. Figures 4A,B show representative responses of the forward speed for the models without and with phase resetting, respectively, after the models were disturbed. For both models, the forward speed fluctuated after the perturbation and then recovered to steady periodic behavior. Although the model without phase resetting had no shift of the locomotion phase after the recovery, the model with phase resetting had a phase shift.

Responses of forward speed for the model (A) without phase resetting and the model (B) with phase resetting.

From 100 trials with different perturbation phases, we obtained the PRCs. Figures 5A–C show the PRCs calculated by acceleration and deceleration perturbations for the model without phase resetting, the model with phase resetting, and kinematic measurements of human walking, respectively. Although the PRCs for the model without phase resetting were zero irrespective of the perturbation phase for both types of perturbation, the PRCs for the model with phase resetting showed characteristic phase-dependent shapes. In particular, they intersected with the horizontal axis around the foot-contact timings and mid-stance phases. They had steep positive slopes around the foot contact and positive peaks in the double-stance phase. They also had gentle negative slopes after the double-stance phase and negative peaks before the next foot contact. Furthermore, the PRCs for the acceleration and deceleration perturbations were almost identical. These trends are similar to those in the PRCs obtained from the measurement of human walking.

PRCs calculated for acceleration and deceleration perturbations for (A) the model without phase resetting, (B) the model with phase resetting, and (C) kinematic measurements of human walking. (C) is modified from Funato et al. (2016). 0 and 100% of the gait cycle represent right foot contact, and gray regions indicate the double-stance phase.

To investigate the robustness of the obtained results, we examined how the PRC changes for different belt speeds and different motor control parameters, such as the gait frequency. Figures 6A,B show the PRCs for the model with phase resetting when the steady belt speed was increased and decreased by 0.02 m/s, respectively. Although there are some differences, the characteristic properties mentioned above remain unchanged.

PRCs for acceleration and deceleration perturbations for the model with phase resetting when the steady belt speed was (A) increased and (B) decreased by 0.02 m/s. 0 and 100% of the gait cycle represent right foot contact, and gray regions indicate the double-stance phase.

 4. Discussion 4.1. Contribution of Phase Resetting Based on Foot Contact

Our motor control model incorporated phase resetting induced by foot-contact information based on physiological evidence. In particular, cutaneous feedback was observed to contribute to phase shift and rhythm resetting behaviors in fictive locomotion of decerebrate cats (Duysens, 1977; Schomburg et al., 1998). Furthermore, spinal cats walking on a treadmill changed their gait, such as walking, trotting, and galloping, in accordance with the belt speed (Forssberg and Grillner, 1973; Orlovsky et al., 1999), which suggests that the tactile sensory information obtained by their feet from the belt influenced the locomotion phase and rhythm generated by the CPG (Duysens et al., 2000). In human walking, the timing of basic muscle activation patterns was also strictly linked to the foot-contact event (Ivanenko et al., 2006). Our model integrated the phase resetting mechanism with the muscle synergy hypothesis that suggested that a large portion of motor commands for walking is generated by the linear combination of five basic signals to solve the redundancy problem in motor control (Ivanenko et al., 2006). More specifically, phase resetting in our model just controlled the timing of the basic signals to determine the motor command using the foot-contact information, which is a very simple strategy. Despite the simple strategy, this timing regulation of the basic signals has been reported to produce various locomotor functions, such as the walk-run transition (Cappellini et al., 2006; Aoi et al., 2019), stepping over an obstacle (Ivanenko et al., 2005; Aoi et al., 2013), and split-belt treadmill walking (MacLellan et al., 2014; Fujiki et al., 2018). Hodgkin-Huxley style neuron model showed that the phase of the neurons' activity rapidly changed by external signals (Rybak et al., 2006a,b), which suggests that the neural system has a mechanism to quickly move the phase of neurons' activity. We would like to incorporate a more biologically detailed neuron model to further investigate the contribution of phase resetting in the future.

Electrical stimulation to the swing legs in cats (Forssberg et al., 1975; Forssberg, 1979) and humans (Belanger and Patla, 1984; Duysens et al., 1990) and mechanical stimulation in humans (Schillings et al., 1996, 2000) showed phase-dependent responses. In particular, stimulation early in the swing phase enhanced flexor muscle activities and extended the swing phase (elevating strategy), while stimulation late in the swing phase enhanced extensor muscle activities and advanced the foot-contact timing (lowering strategy) (Eng et al., 1994). From the intersection of the obtained PRCs with zero (Figure 5C), our previous work (Funato et al., 2016) showed that the mid-single stance phase extended in response to acceleration perturbations and the foot-contact timing advanced in response to deceleration perturbations, which correspond to the elevating and lowering strategies, respectively. In this study, we incorporated the phase resetting mechanism that modulates the timing of the motor command based on the foot-contact information. This mechanism is related to the lowering strategy. Despite not incorporating the elevating strategy, our model had a PRC shape similar to that for humans not only for deceleration perturbations but also for acceleration perturbations (Figures 5B,C). That is, application of only one of these two strategies allowed the model to reproduce the PRCs for acceleration and deceleration perturbations in humans. In the future, we would like to incorporate the elevating strategy in our motor control model to further clarify the adaptive rhythm control mechanism in human walking.

 4.2. Calculation of PRC

To calculate the PRC from kinematic measurements, mainly two methods have been proposed. One is the impulse method that uses single-impulse perturbation, and the other is the weighted spike-triggered average (WSTA) method, which uses sequential pulse perturbation with zero mean and no temporal correlation (Ota et al., 2009). Our previous work (Funato et al., 2016) used both of these methods to calculate the PRCs for human walking. Because the impulse method required many trials that exhausted the subjects, the obtained PRCs had large deviations and low temporal resolution and could not show characteristic properties. The WSTA method improved the PRCs, and clear phase-dependent shapes could be resolved (Figure 5C). However, it still has limitations with regard to obtaining precise PRCs. For example, two positive peak timings of the PRC for the deceleration perturbation differed. In addition, the acceleration and deceleration perturbations showed some differences in the PRC, and it was difficult to determine whether they were actually different or due to limitations of the method. In particular, although the PRC was analytically derived under the assumption that the perturbation is sufficiently small, the perturbation must be large to reduce the influence of measurement noise in human experiments. In this study, we used a mathematical model to obtain the PRCs for human walking. Because of the high reproducibility of the simulation results, we obtained accurate PRCs for the model using arbitrarily small and short perturbations by the impulse method. Our model showed identical PRCs for the acceleration and deceleration perturbations. The modeling approach using the PRC has an advantage for improving our understanding of the underlying rhythm control mechanism.

 4.3. Limitations of Our Model and Future Work

The PRC for the model with phase resetting had a similar shape to that of the PRC for humans, and it supports the hypothesis that phase resetting contributes to adaptive rhythm control in human walking in comparison with the PRC for the model without phase resetting. However, our model has limitations for accurately reproducing the PRC for human walking. For example, the PRC for the model with phase resetting had much steeper positive peaks in the double-stance phases compared to the PRC for humans (Figures 5B,C). This is possibly because four discrete points on each sole were used for the foot-contact model. Due to the discrete points, perturbations in double-stance phases induced sudden changes in locomotor behavior and caused the steep positive peaks in the PRC. In addition, our model showed short double-stance phases compared to actual human walking (Figures 5B,C), which is mainly due to no phalangeal joint in our foot model. However, the PRCs for the model and humans had similar characteristics in the double-stance phase, such as steep positive slopes around the foot contact and positive peaks located in the double-stance phase.

Although we incorporated the phase resetting mechanism in the motor control model, other sensory-motor coordination mechanisms also play a role in human walking. For example, although we focused on the swing-to-stance phase transition using the foot-contact information, the stance-to-swing phase transition has been suggested to include important sensory-motor coordination mechanisms (Ekeberg and Pearson, 2005; Pearson et al., 2006; Dzeladini et al., 2014; Song and Geyer, 2015), such as the unloading rule that uses force-sensitive afferents in the ankle extensor muscles (Duysens and Pearson, 1980; Whelan et al., 1995) and the hip extension rule that uses position-sensitive afferents from the hip (Grillner and Rossignol, 1978; Hiebert et al., 1996). In addition, although this study changed the belt speed of a treadmill to disturb human locomotor behavior, other types of perturbations, such as pulling on the swing leg, have been used (Kobayashi et al., 2000; Nessler et al., 2016). Because the PRC depends on the perturbation, we would like to incorporate other sensory-motor coordination mechanisms and perturbations to further clarify adaptive rhythm control in human walking.

 Data Availability Statement

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

 Author Contributions

SA developed the study design. DT performed computer simulations and analyzed the data in consultation with SA, TF, SF, KS, and KT. DT and SA wrote the manuscript and all the authors reviewed and approved it.

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. The handling Editor declared a shared affiliation, though no other collaboration, with some of the authors DT, SA, KS, and KT.

 References Aoi S. Funato T. (2016). Neuromusculoskeletal models based on the muscle synergy hypothesis for the investigation of adaptive motor control in locomotion via sensory-motor coordination. Neurosci. Res. 104, 88–95. 10.1016/j.neures.2015.11.00526616311 Aoi S. Kondo T. Hayashi N. Yanagihara D. Aoki S. Yamaura H. . (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. 10.1007/s00422-013-0546-623430278 Aoi S. Ogihara N. Funato T. Sugimoto Y. 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. 10.1007/s00422-010-0373-y20217427 Aoi S. Ohashi T. Bamba R. Fujiki S. Tamura D. Funato T. . (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. 10.1038/s41598-018-37460-330674970 Belanger M. Patla A.E. (1984). Corrective responses to perturbation applied during walking in humans. Neurosci. Lett. 49, 291–295. 6493612 Burke R.E. Degtyarenko A.M. Simon E.S. (2001). Patterns of locomotor drive to motoneurons and last-order interneurons: clues to the structure of the CPG. J. Neurophysiol. 86, 447–462. 10.1152/jn.2001.86.1.44711431524 Cappellini G. Ivanenko Y.P. Poppele R.E. Lacquaniti F. (2006). Motor patterns in human walking and running. J. Neurophysiol. 95, 3426–3437. 10.1152/jn.00081.200616554517 Conway B.A. Hultborn H. Kiehn O. (1987). Proprioceptive input resets central locomotor rhythm in the spinal cat. Exp. Brain Res. 68, 643–656. 3691733 Davy D.T. Audu M.L. (1987). A dynamic optimization technique for predicting muscle forces in the swing phase of gait. J. Biomech. 20, 187–201. 10.1016/0021-9290(87)90310-13571299 Duysens J. (1977). Fluctuations in sensitivity to rhythm resetting effects during the cat's step cycle. Brain Res. 133, 190–195. 902087 Duysens J. Clarac F. Cruse H. (2000). Load-regulating mechanisms in gait and posture: Comparative aspects. Physiol. Rev. 80, 83–133. 10.1152/physrev.2000.80.1.8310617766 Duysens J. Pearson K.G. (1980). Inhibition of flexor burst generation by loading ankle extensor muscles in walking cats. Brain Res. 187, 321–332. 7370733 Duysens J. Trippel M. Horstmann G.A. Dietz V. (1990). Gating and reversal of reflexes in ankle muscles during human walking. Exp. Brain Res. 82, 351–358. 2286237 Dzeladini F. van den Kieboom J. Ijspeert A. (2014). The contribution of a central pattern generator in a reflex-based neuromuscular model. Front. Hum. Neurosci. 8:371. 10.3389/fnhum.2014.0037125018712 Ekeberg Ö 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, 4256–4268. 10.1152/jn.00065.200516049149 Eng J.J. Winter D.A. Patla A.E. (1994). Strategies for recovery from a trip in early and late swing during human walking. Exp. Brain Res. 102, 339–349. 7705511 Forssberg H. (1979). Stumbling corrective reaction: a phase-dependent compensatory reaction during locomotion. J. Neurophysiol. 42, 936–953. 479924 Forssberg H. Grillner S. (1973). The locomotion of the acute spinal cat injected with clonidine i.v. Brain Res. 50, 184–186. 4690545 Forssberg H. Grillner S. Rossignol S. (1975). Phase dependent reflex reversal during walking in chronic spinal cats. Brain Res. 85, 103–107. 1109686 Frigon A. Gossard J.-P. (2010). Evidence for specialized rhythm-generating mechanisms in the adult mammalian spinal cord. J. Neurosci. 30, 7061–7071. 10.1523/JNEUROSCI.0450-10.201020484648 Fujiki S. Aoi S. Funato T. Sato Y. Tsuchiya K. Yanagihara D. (2018). Adaptive hindlimb split-belt treadmill walking in rats by controlling basic muscle activation patterns via phase resetting. Sci. Rep. 8:17341. 10.1038/s41598-018-35714-830478405 Funato T. Yamamoto Y. Aoi S. Imai T. Aoyagi T. Tomita N. Tsuchiya K. (2016). Evaluation of the phase-dependent rhythm control of human walking using phase response curves. PLoS Comput. Biol. 12:e1004950. 10.1371/journal.pcbi.100495027203839 Grillner S. Rossignol S. (1978). On the initiation of the swing phase of locomotion in chronic spinal cats. Brain Res. 146, 269–277. 274169 Guertin P. Angel M.J. Perreault M.-C. 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. 7473249 Hiebert G.W. Whelan P.J. Prochazka A. Pearson K.G. (1996). Contribution of hindlimb flexor muscle afferents to the timing of phase transitions in the cat step cycle. J. Neurophysiol. 75, 1126–1137. Imai T. Aoyagi T. (2016). Improvement effect of measuring phase response curves by using multicycle data. NOLTA IEICE 7, 58–65. 10.1587/nolta.7.58 Ivanenko Y.P. Cappellini G. Dominici N. Poppele R.E. Lacquaniti F. (2005). Coordination of locomotion with voluntary movements in humans. J. Neurosci. 25, 7238–7253. 10.1523/JNEUROSCI.1327-05.200516079406 Ivanenko Y.P. Poppele R.E. Lacquaniti F. (2006). Motor control programs and walking. Neuroscientist 12, 339–348. 16840710 Kobayashi M. Nomura T. Sato S. (2000). Phase-dependent response during human locomotion to impulsive perturbation and its interpretation based on neural mechanism. Jpn. J. Med. Electr. Biol. Eng. 38, 20–32. 10.11239/jsmbe1963.38.20 Kuramoto Y. (1984). Chemical Oscillations, Waves, and Turbulence. New York, NY: Springer. Lafreniere-Roula M. McCrea D.A. (2005). Deletions of rhythmic motoneuron activity during fictive locomotion and scratch provide clues to the organization of the mammalian central pattern generator. J. Neurophysiol. 94, 1120–1132. 10.1152/jn.00216.200515872066 MacLellan M.J. Ivanenko Y.P. Massaad F. Bruijn S.M. Duysens J. Lacquaniti F. (2014). Muscle activation patterns are bilaterally linked during split-belt treadmill walking in humans. J. Neurophysiol. 111, 1541–1552. 10.1152/jn.00437.201324478155 Nessler J.A. Spargo T. Craig-Jones A. Milton J.G. (2016). Phase resetting behavior in human gait is influenced by treadmill walking speed. Gait Posture 43, 187–191. 10.1016/j.gaitpost.2015.09.02126475759 Nomura T. Kawa K. Suzuki Y. Nakanishi M. Yamasaki T. (2009). Dynamic stability and phase resetting during biped gait. Chaos 19:026103. 10.1063/1.313872519566263 Orlovsky G.N. Deliagina T. Grillner S. (1999). Neuronal Control of Locomotion: From Mollusc to Man. Oxford: Oxford University Press. Ota K. Nomura M. Aoyagi T. (2009). Weighted spike-triggered average of a fluctuating stimulus yielding the phase response curve. Phys. Rev. Lett. 103:24101. 10.1103/PhysRevLett.103.02410119659207 Pearson K. Ekeberg Ö Büschges A. (2006). Assessing sensory function in locomotor systems using neuro-mechanical simulations. Trends Neurosci. 29, 625–631. 10.1016/j.tins.2006.08.00716956675 Rybak I.A. Shevtsova N.A. Lafreniere-Roula M. McCrea D.A. (2006a). Modelling spinal circuitry involved in locomotor pattern generation: insights from deletions during fictive locomotion. J. Physiol. 577, 617–639. 10.1113/jphysiol.2006.11870317008376 Rybak I.A. Stecina K. Shevtsova N.A. McCrea D.A. (2006b). Modelling spinal circuitry involved in locomotor pattern generation: insights from the effects of afferent stimulation. J. Physiol. 577, 641–658. 10.1113/jphysiol.2006.11871117008375 Schillings A.M. van Wezel B.M. Duysens J. (1996). Mechanically induced stumbling during human treadmill walking. J. Neurosci. Methods 67, 11–17. 8844520 Schillings A.M. van Wezel B.M. Mulder T. Duysens J. (2000). Muscular responses and movement strategies during stumbling over obstacles. J. Neurophysiol. 83, 2093–2102. 10.1152/jn.2000.83.4.209310758119 Schomburg E.D. Petersen N. Barajon I. Hultborn H. (1998). Flexor reflex afferents reset the step cycle during fictive locomotion in the cat. Exp. Brain Res. 122, 339–350. 9808307 Song S. Geyer H. (2015). A neural circuitry that emphasizes spinal feedback generates diverse behaviours of human locomotion. J. Physiol. 593, 3493–3511. 10.1113/JP27022825920414 Whelan P.J. Hiebert G.W. Pearson K.G. (1995). Stimulation of the group I extensor afferents prolongs the stance phase in walking cats. Brain Res. 103, 20–30. 10.1007/BF002419617615034 Winfree A.T. (2001). The Geometry of Biological Time. New York, NY: Springer. Winter D.A. (2004). Biomechanics and Motor Control of Human Movement, 3rd edn. New York, NY: Wiley. Yamasaki T. Nomura T. Sato S. (2003a). Possible functional roles of phase resetting during walking. Biol. Cybern. 88, 468–496. 10.1007/s00422-003-0402-112789495 Yamasaki T. Nomura T. Sato S. (2003b). Phase reset and dynamic stability during human gait. BioSystems 71:, 221–232. 10.1016/S0303-2647(03)00118-714568223

Funding. This study was supported in part by JSPS KAKENHI Grant Numbers JP15KT0015, JP26120006, and JP17H04914.