Contribution of Phase Resetting to Adaptive Rhythm Control in Human Walking Based on the Phase Response Curves of a Neuromusculoskeletal Model

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.

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.
Keywords: human walking, phase resetting, phase response curve, central pattern generator, muscle synergy, neuromusculoskeletal model

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 , 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 . 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.

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.

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:

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 . More specifically, we used five rectangular pulses p i (φ) (i = 1, . . . , 5) for each leg (Figure 1B), which are given by: 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 u Syn m by: where w m,i (i = 1, . . . , 5) is the weighting coefficient of five rectangular pulses (w m,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: where, ω is the basic frequency, K φ is the gain parameter, t FC i 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.
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.
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 u Reg m is given by: 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 u m is given by the summation of the muscle synergy-based motor command u

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:ψ 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.
When a perturbation I(t) is added to the limit cycle oscillator, ψ follows the dynamics: 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 = t n (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 t n gives: 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 t n 0ψ dt = 2nπ . However, because the disturbed trajectory approaches C as t → ∞, t n 0ψ dt = 2nπ approximately for sufficiently large n. This gives: 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: This gives, 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  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.

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.
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 doublestance 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.
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.

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 .  In human walking, the timing of basic muscle activation patterns was also strictly linked to the foot-contact event . 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 . 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  Funato et al. (2016). 0 and 100% of the gait cycle represent right foot contact, and gray regions indicate the double-stance phase.
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(Schillings et al., , 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 footcontact timing (lowering strategy) (Eng et al., 1994). From the intersection of the obtained PRCs with zero ( Figure 5C), our previous work  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.

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  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.

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 doublestance 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 sensorymotor 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 forcesensitive 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.

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