
==== Front
Commun Biol
Commun Biol
Communications Biology
2399-3642
Nature Publishing Group UK London

6843
10.1038/s42003-024-06843-w
Article
Interlimb coordination is not strictly controlled during walking
http://orcid.org/0000-0003-3153-5377
Arai Takahiro 1
Ota Kaiichiro 2
Funato Tetsuro 3
Tsuchiya Kazuo 4
Aoyagi Toshio 5
http://orcid.org/0000-0001-9243-2641
Aoi Shinya aoi.shinya.es@osaka-u.ac.jp

6
1 https://ror.org/059qg2m13 grid.410588.0 0000 0001 2191 0132 Center for Mathematical Science and Advanced Technology, Japan Agency for Marine-Earth Science and Technology, Yokohama, 236-0001 Japan
2 Cybozu, Inc., 2-7-1 Nihombashi, Chuo-ku, Tokyo 103-6027 Japan
3 https://ror.org/02x73b849 grid.266298.1 0000 0000 9271 9936 Department of Mechanical Engineering and Intelligent Systems, Graduate School of Informatics and Engineering, The University of Electro-Communications, 1-5-1 Choufugaoka, Choufu, Tokyo 182-8585 Japan
4 https://ror.org/02kpeqv85 grid.258799.8 0000 0004 0372 2033 Department of Aeronautics and Astronautics, Graduate School of Engineering, Kyoto University, Kyoto daigaku-Katsura, Nishikyo-ku, Kyoto, 615-8540 Japan
5 https://ror.org/02kpeqv85 grid.258799.8 0000 0004 0372 2033 Graduate School of Informatics, Kyoto University, Yoshida-Honmachi, Sakyo-ku, Kyoto, 606-8501 Japan
6 https://ror.org/035t8zc32 grid.136593.b 0000 0004 0373 3971 Department of Mechanical Science and Bioengineering, Graduate School of Engineering Science, The University of Osaka, 1-3 Machikaneyama, Toyonaka, Osaka 560-8531 Japan
20 9 2024
20 9 2024
2024
7 115219 2 2023
4 9 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
In human walking, the left and right legs move alternately, half a stride out of phase with each other. Although various parameters, such as stride frequency and length, vary with walking speed, the antiphase relationship remains unchanged. In contrast, during walking in left-right asymmetric situations, the relative phase shifts from the antiphase condition to compensate for the asymmetry. Interlimb coordination is important for adaptive walking and we expect that interlimb coordination is strictly controlled during walking. However, the control mechanism remains unclear. In the present study, we derived a quantity that models the control of interlimb coordination during walking using two coupled oscillators based on the phase reduction theory and Bayesian inference method. The results were not what we expected. Specifically, we found that the relative phase is not actively controlled until the deviation from the antiphase condition exceeds a certain threshold. In other words, the control of interlimb coordination has a dead zone like that in the case of the steering wheel of an automobile. It is conjectured that such forgoing of control enhances energy efficiency and maneuverability. Our discovery of the dead zone in the control of interlimb coordination provides useful insight for understanding gait control in humans.

Analyses of the synchronization mechanism between the legs during walking using two coupled oscillators and Bayesian inference method reveal the presence of a “dead zone,” where the control of interlimb coordination is ineffective.

Subject terms

Nonlinear dynamics
Oscillators
501100003382 MEXT | JST | Core Research for Evolutional Science and Technology (CREST) JPMJCR09U2 MEXT | JST | FOREST Program Grant Number JPMJFR2021501100003382 MEXT | JST | Core Research for Evolutional Science and Technology (CREST) JPMJCR09U2 501100001691 MEXT | Japan Society for the Promotion of Science (JSPS) JP20K21810 501100001691 MEXT | Japan Society for the Promotion of Science (JSPS) JP20H04144 501100001691 MEXT | Japan Society for the Promotion of Science (JSPS) JP20K20520 501100001691 MEXT | Japan Society for the Promotion of Science (JSPS) JP23H04467 issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

During human walking, the left and right legs move alternately, half a stride out of phase with each other1. In general, various locomotion parameters, such as gait frequency, stride length, and duty factor, change if the gait speed varies. However, the antiphase relationship of the leg motion remains unchanged. This is true even if the gait pattern changes to running. By contrast, during walking in which there is left-right asymmetry with regard to the body or environment, such as walking with unilateral leg loading2, walking along a curved path3, or walking on a split-belt treadmill4, the relative phase between the legs shifts from the antiphase condition to compensate for the asymmetry. In addition, the phase relationship fluctuates significantly during walking of elderly people and patients with neurological disabilities, such as those caused by stroke or Parkinson’s disease5–9. These findings indicate that appropriate relative phase (that is, appropriate interlimb coordination) is important for adaptive walking. This seems to suggest that the relative phase is strictly controlled to remain the appropriate relative phase during walking of healthy young people, with the relative phase quickly returning to the appropriate relative phase after being perturbed away from it.

To understand interlimb coordination during human walking, it has been investigated how the relative phase between the motion of the left and right legs depends on the situations described above. However, it remains largely unclear to what extent it is strictly controlled in each situation. This is mainly because there are limitations on the degree to which the control of interlimb coordination can be quantitatively elucidated due to the complexity of neural dynamics, musculo-skeletal dynamics, and interactions with the environment during walking. Elucidating the control of interlimb coordination would not only clarify the adaptive strategy in human walking, but also contribute to the development of rehabilitation techniques for persons with gait disorders and gait assistive devices, such as robotic exoskeletons10–12.

The present study aims to quantitatively clarify the control of interlimb coordination during human walking by taking advantage of both data science and dynamical systems theory, employing the Bayesian inference method13,14 and phase reduction theory15–20. Specifically, we model the rhythmic motion of the left and right legs by two coupled limit-cycle oscillators and describe these dynamics using two coupled phase oscillators. The phase coupling function between the two oscillators represents the nature of the control of the relative phase between the motion of the legs, i.e., interlimb coordination. To identify this phase coupling function in human walking, we first measured the leg motion during walking of subjects on a treadmill with intermittent perturbations in the belt speed. Next, from the measured time-series data, we derived the phase coupling function using the Bayesian inference method. We analyze the characteristics of the control of interlimb coordination during human walking on the basis of this function.

Results

Walking experiment

We measured kinematic data through observation of eight healthy men (Subjects A–H) during walking on a treadmill whose belt speed was intermittently disturbed from 1.0 m/s for a short time (Fig. 1a). More precisely, the belt speed was temporarily increased or decreased by 0.6 m/s over a period of 0.1 s, and then returned to 1.0 m/s over a period of 0.1 s. We regard the disturbance of the belt speed from 1.0 m/s as an external perturbation, written by p(t). This perturbation resulted in the sudden collapse of the antiphase relationship of the leg motion and allowed us to investigate how the relative phase of the leg motion recovered after being disturbed. Each trial lasted approximately 60 s and contained ten disturbances during which the speed was varied. These disturbances were separated by intervals of approximately 5 s, which was sufficiently long to allow the return to normal steady-state walking. We used three types of perturbation conditions for each trial: acceleration, deceleration, and mixed conditions. Only temporary acceleration and deceleration perturbations of the belt speed were applied under the acceleration and deceleration conditions, respectively, while both perturbations were applied randomly under the mixed condition. From the measured kinematic data, we calculated the elevation angle of the limb axis, which connects the hip and toe of each leg. The angles are denoted by φL(t) and φR(t), where the subscripts L and R indicate the left leg and right leg, respectively.Fig. 1 Schematic diagram of the study.

a Measurement of the leg motion, represented by φi(t) (i ∈ (L, R)), during walking on a treadmill with a perturbation in the belt speed, represented by p(t). b Two coupled limit-cycle oscillators with phases ϕi(t) (i ∈ (L, R)) whose dynamics are described by phase equations with phase coupling functions Γij(ϕj − ϕi), phase sensitivity functions Zi(ϕi), external perturbation I(t), and Gaussian noise ξi ((i, j) = (L, R), (R, L)). c Integration of the measured data and phase equations to derive the dynamics of the relative phase ΔLR, which reflects the control of interlimb coordination. When ΔLR = π, corresponding to B, the legs move in an exactly alternating manner. As ΔLR moves away from π, the deviation from this antiphase relationship increases (see A and C, with ΔLR = ΔA and ΔC). If the interlimb coordination is strictly controlled to maintain the antiphase relationship, the function f(ΔLR), which describes the control of interlimb coordination, should intersect the horizontal axis at π with a steep negative slope.

Modeling the control of interlimb coordination using coupled phase oscillators

We model the motion of the left and right legs by two coupled limit-cycle oscillators, whose phases are ϕL and ϕR ∈ (0, 2π], and describe these dynamics using phase equations on the basis of the phase reduction theory15–18 as follows (Fig. 1b):1 ϕ°L(t)=ωL+ΓLR(ϕR−ϕL)+ZL(ϕL)I(t)+ξL(t),

2 ϕ°R(t)=ωR+ΓRL(ϕL−ϕR)+ZR(ϕR)I(t)+ξR(t).

Here, ωL and ωR are the natural frequencies, ΓLR and ΓRL are phase coupling functions, ZL and ZR are phase sensitivity functions, I is an external perturbation (rescaled from p(t)), and ξL and ξR are independent Gaussian white noise terms that represent experimental uncertainty, which satisfy 〈ξi(t)〉 = 0 and ⟨ξi(t)ξj(s)⟩=σi2δijδ(t−s) for (i, j) = (L, R) and (R, L), where σi is the intensity of the Gaussian noise ξi, and δij and δ(t) are the Kronecker and Dirac delta functions, respectively. We determined ϕi(t) (i ∈ (L, R)) from the measured time series of the kinematic angle φi(t), and also determined I(t) from the measured time series of the external perturbation p(t). We determined the best values for ωi, Γij, and Zi ((i, j) = (L, R), (R, L)) from the time series of ϕi(t) and I(t) using the Bayesian inference method14,17,21,22. We did the same for σi in order to determine the extent to which Eqs. (1) and (2) are capable of modeling the dynamics of the observed human walking. We deem that a small intensity indicates that our deterministic model properly describes the dynamics and that the influence of the stochastic process is small. In previous studies14,21,22, ZiI was not used in the identification of the phase equation, and the parameter values were obtained using data measured after external perturbations were applied. However, including ZiI in this identification allowed us to determine the parameter values more accurately through use of data measured while external perturbations were being applied.

Using phase equations derived through the procedure described above, we investigate how the relative phase between the motion of the legs, i.e., interlimb coordination, is controlled during human walking. Specifically, if we ignore the phase response to the external perturbation in Eqs. (1) and (2) and focus on the transient dynamics after the perturbation is applied, we obtain the following equation for the relative phase ΔLR: = ϕR − ϕL:3 Δ°LR=ωR−ωL+ΓRL(−ΔLR)−ΓLR(ΔLR)+ξR−ξL:=f(ΔLR)+ξLR.

Here, we have ξLR: = ξR − ξL, which satisfies 〈ξLR(t)〉 = 0 and ⟨ξLR(t)ξLR(s)⟩=σLR2δ(t−s), where σLR:=σR2+σL2, and f(ΔLR) describes the deterministic dynamics of ΔLR, i.e., the control of interlimb coordination. If the antiphase relationship, ΔLR = π, is tightly maintained, it is required that f(ΔLR) be 0 and have a steep negative slope at ΔLR = π (Fig. 1c).

Figure 2 displays the time series of φi(t), ϕi(t), ΔLR(t), and p(t) (i ∈ (L, R)) for one representative subject (Subject G). There, an increase of ϕi(t) by 2π represents one cycle of φi(t). It is seen that ΔLR(t) is greatly disturbed in response to the external perturbation (highlighted in green) and then returns to the neighborhood of π, the antiphase value. However, it does not converge to π but rather fluctuates about this value. This suggests that the dynamics of the relative phase possesses Lyapunov stability about π, not asymptotic stability. We now proceed to elucidate how f(ΔLR) reflects this stability.Fig. 2 Measured and reconstructed time series for one representative subject (Subject G) under mixed conditions.

The leg motion, represented by φi(t) (i ∈ (L, R)), measured during walking was transformed to the oscillator phases ϕi(t) (i ∈ (L, R)), from which we obtain the relative phase ΔLR. The disturbance of the treadmill belt speed from 1.0 m/s is represented as p(t). The green regions indicate the times at which temporary acceleration or deceleration is applied to the treadmill belt speed.

Figure 3a, b display the functions Zi (i ∈ (L, R)) and f(ΔLR), respectively, evaluated from the data for one representative subject (Subject G). It is seen that both ZL and ZR possess unimodal shapes with peaks near π and are close to zero in the range from 0 to π/2 for all perturbation conditions. Because ΔLR fluctuates within a narrow region around π, as shown in Fig. 2, there was not a sufficient amount of data outside this region for the evaluation of f(ΔLR). Therefore, we limited our evaluation of f(ΔLR) to a range of 3 standard deviations around the mean of the observed ΔLR. Although we expected that f(ΔLR) would intersect the horizontal axis at one point near π with a steep negative slope, in accordance with the strict control of interlimb coordination, as shown in Fig. 1c, the results were not what we expected for any of the perturbation conditions. Specifically, we found that although f(ΔLR) has a steep negative slope sufficiently far from π, it is close to zero near π. Thus, we found that f(ΔLR) has a flat region with a value close to 0 in the neighborhood of the antiphase condition. This result indicates that the relative phase between the legs is not actively controlled until the deviation from the antiphase condition exceeds a certain threshold, where the antiphase relationship is lost as shown in Fig. 3c. We obtained similar results from the other subjects (see Supplementary Note 1 and Supplementary Fig. 1). These characteristic properties of f(ΔLR) remained even when we excluded ZiI from the phase equation in the determination of f(ΔLR) (see Supplementary Note 2 and Supplementary Fig. 2). However, the form of f(ΔLR) obtained with ZiI exhibits these characteristics more clearly than that without ZiI.Fig. 3 Results for the control of interlimb coordination for one representative subject (Subject G).

a Phase sensitivity functions Zi(ϕi) (i ∈ (L, R)) for acceleration (Acc), deceleration (Dec), and mixed (Mix) conditions. b Control of interlimb coordination, represented by f(ΔLR). The vertical dotted lines indicate 3 standard deviations from the mean of the observed ΔLR. In the region outside the vertical dotted line, the estimated function f(ΔLR) is displayed in pale color. The function f(ΔLR) is shifted in order to place the mean of the observed values of ΔLR at π to improve visualization. The histogram displays the noise intensity, σLR. (An expanded histogram, with values 30 times larger than the actual values, is also given.) The mean and standard deviation of the derived noise intensity for all results are 0.0106π and 0.0021π, respectively. These values are too small to cause deviation of ΔLR from the flat region of f(ΔLR). c A schematic diagram for the interpretation of f(ΔLR). The flat region near ΔLR = π at B (neutral stability) and the steep negative slope in the regions away from π reflect the fact that the relative phase between the legs is not actively controlled until the deviation from π exceeds a threshold (ΔA − π at A and ΔC − π at C), where motion of the legs deviates significantly from the antiphase relationship.

Evaluating the flat region in the control of interlimb coordination

To clarify the universality of our findings in the control of interlimb coordination represented by f(ΔLR), specifically, the existence of a flat region in which f(ΔLR) is close to 0 in the neighborhood of ΔLR = π, and steep negative slopes outside this flat region (Fig. 3c), we compared the results for f(ΔLR) among the subjects. For this purpose, we approximated f(ΔLR) by a piecewise linear function to characterize these properties within 3 standard deviations around the mean of the observed ΔLR, as shown in Fig. 4a. There, lC is the length of the flat region, which reflects the extent to which the deviation from the antiphase condition is ignored, and gL and gR are the slopes in the regions to the left and right of the flat region, respectively, which reflect how rapidly the deviation outside the flat region attenuates. These three parameters, lC, gL, and gR, were determined to minimize the discrepancy between the original and approximated functions using the grid search method under the condition that the two functions coincide at the left and right endpoints of the approximated function. When there is no flat region in the original function, lC = 0 is satisfied. Figure 4b displays the approximate forms of f(ΔLR) for the original forms appearing in Fig. 3b. (Similar results for all subjects appear in Supplementary Note 3 and Supplementary Fig. 3.) Figure 4c presents the mean and standard deviation of lC among the subjects. It is seen that lC is approximately 0.4 and the deviation among the subjects is small for all perturbation conditions. We confirmed that lC is significantly larger than 0 using one-tailed t-test (p ≪ 0.01, see Supplementary Note 4). Indeed, the mean value of lC is much larger than its standard deviation, and we thus conclude that lC has a positive value. Figure 4d presents the mean and standard deviation of gL and gR among the subjects. It is seen that both ∣gL∣ and ∣gR∣ are approximately 2.0, sufficiently large that deviations rapidly attenuate for all perturbation conditions. We also confirmed that ∣gL∣ and ∣gR∣ are significantly larger than 0 using one-tailed t-test (p ≪ 0.01, see Supplementary Note 4). These results confirm that f(ΔLR) has a flat region with a value close to 0 in the vicinity of ΔLR = π and steep negative slopes outside this region for all perturbation conditions.Fig. 4 Characteristics of the control of interlimb coordination, represented by f(ΔLR).

a Approximation of f(ΔLR) by a piecewise linear function within 3 standard deviations (3SD) (vertical dotted lines) from the mean of the observed ΔLR. The quantity lC is the length of the flat region in which we have f(ΔLR) = 0, and gL and gR are the slopes in the left and right outside regions, respectively. The three parameters lC, gL, and gR were determined by the left and right endpoints of the flat region (white circles) such that the discrepancy between the original and approximated functions is minimized within the region between the two vertical dotted lines under the condition that the approximated function and the original function coincide at the left and right endpoints of the approximated function (black circles). b Result for the approximation of f(ΔLR) in Fig. 3b under acceleration (Acc), deceleration (Dec), and mixed (Mix) conditions. c Mean and standard deviation of lC among the subjects, where “All” indicates the statistical result from all perturbation conditions. The dots along the bars represent data points (n = 8 for Acc and Dec, n = 6 for Mix, and n = 22 for “All”). d Mean and standard deviation of gL and gR among the subjects.

Discussion

Using time-series data measured during walking under perturbations that disturbed the antiphase relationship between the motion of the left and right legs (Fig. 2), we constructed a function that models the control of interlimb coordination in human walking on the basis of the phase reduction theory and Bayesian inference method. We found that this function has a flat region with values close to 0 in the vicinity of the antiphase condition and steep negative slopes outside this flat region for all subjects studied (Fig. 4). This result indicates that the relative phase between the motion of the left and right legs is not actively controlled until the deviation from the antiphase condition exceeds a certain threshold. In other words, there is a dead zone like that in the case of the steering wheel of an automobile. This differs from our expectation that the relative phase is strictly controlled to maintain the antiphase condition during walking (Fig. 1c). This unexpected finding was obtained for the first time through the quantitative evaluation of the control of interlimb coordination as a function of the relative phase.

Such a forgoing of control implies free for energy cost for the control. This conjectures that the dead zone in the control of interlimb coordination could contribute to energy efficiency during walking. In addition, the forgoing of control induced the relative phase between the legs to be neither exponentially stable nor unstable, but neutral around the antiphase condition (Fig. 3c). Such a neutral stability could enhance maneuverability to change the walking behavior. For example, walking along a curved path requires left-right asymmetry in the walking behavior and a phase shift from the antiphase condition3, which conjectures that the lack of control allows a change in the walking direction. However, large deviations from the antiphase condition could result in a deterioration of gait performance. For this reason, it is conjectured that when the deviation exceeds a certain threshold, control is activated to reduce this deviation. Although the slopes of f(ΔLR) outside the flat region vary among the subjects (Fig. 4d), the length of the flat region, which determines the threshold beyond which the control is exercised, does not (Fig. 4c). This result suggests that the threshold is a universal characteristic inherent in the normal walking of healthy young people.

The present result was obtained by determining the specific shape of the function controlling the interlimb coordination. The method employed here relies to a large extent on recent developments in the identification of the phase equation using big data14. In previous studies23,24, interlimb coordination has also been investigated through use of the phase equation. However, in those works, a simple form (for example a sine function) was assumed for the coupling function in order to determine the function controlling the interlimb coordination. In that case, there appeared no region in which control is not exercised.

The forgoing of control in specific situations in which we expect strict control also appears during quiet standing in humans. Because the upright standing posture is inherently unstable, like an inverted pendulum, it is natural to conjecture that posture is strictly controlled, in order to maintain quiet standing, as we also conjectured that the relative phase of leg motion is strictly controlled in order to maintain the antiphase condition during walking. However, it has been suggested that posture is not actively controlled until deviation from the quiet, upright state exceeds a state-dependent threshold25–28. It is thus conjectured that interlimb coordination during walking and quiet standing have a common motor control strategy in humans.

The leg motion consists of the motion of the thigh, shank, and foot. However, the leg motion is nearly constrained to a two-dimensional space determined by the length and orientation of the limb axis, which connects the hip and toe29. The two-dimensional constraint appears in different speeds30, in other gaits31,32, and in other species33. In addition, some research has reported that there are neural systems that independently represent the length and orientation of the limb axis in cats34,35, suggesting different controls for the length and orientation of the limb axis. We focused on the control of interlimb coordination only about the relationship between the orientations of the left and right legs, and used the measured time series of the limb axis orientation. As shown in Supplementary Note 5, we confirmed that similar results are obtained from multi-dimensional time series of the legs when we focused on the limb axis orientation. We would like to investigate the control of interlimb coordination about not only the limb axis orientation but also the limb axis length in the future.

As mentioned in the previous paragraph, some studies have suggested that the leg motion can be explained in two-dimensional space. However, it is not completely constrained in two-dimensional space and the part outside of the two-dimensional space could affect the phase. Thus, using multipoint measurements may improve the accuracy of phase calculation compared to the method adopted in the present study, enabling us to reveal phase adjustments in detail. There exist some proposed methods to efficiently identify the phase using multipoint measurements20,36. Using such a method, we conducted a preliminary investigation into the control of interlimb coordination and obtained qualitatively similar results, although the flat region was wider than the result shown in this paper. Determining which method is appropriate for investigating the control of leg motion is highly dependent on the control mechanism of human gait and thus remains a topic for future study. However, our key finding—the existence of a dead zone in the interlimb coordination—remains consistent and robust across the methods we examined. This reinforces the validity and significance of our results regardless of the specific method used.

Because each leg is controlled by the corresponding side of the spinal cord, which has the ability to generate rhythm independently37–39, each leg can operate independently. This independence is manifested in the leg movements on a split-belt treadmill with the two belts running at different speeds, in which case the legs exhibit asymmetrical left-right stepping4,40–42. Even though each side of the spinal cord can produce rhythmic movement in each leg separately, the motions of left and right legs are generally coordinated during walking. For example, a reciprocal relationship is maintained between the legs such that only one leg enters the swing phase at a time43. Commissural interneurons contribute to interlimb coordination by mediating interactions between the rhythm-generating locomotor circuits located on each side of the spinal cord44–47. We plan to investigate the neural mechanism governing the dead zone in the control of interlimb coordination in future research. In particular, it would be useful to use the phase reduction method15–18 for a detailed model of the neural system including the musculoskeletal system to derive the two-coupled-oscillator model to evaluate the control of interlimb coordination. The comparison of the results between this theoretical approach and our data-driven approach will provide useful insight into the neural mechanism.

Although the present study focused only on the normal walking of healthy young people, our method has great potential for broader application. For example, because our method can determine the phase coupling functions ΓLR and ΓRL independently, it is applicable to the investigation of the control of interlimb coordination during walking with left-right asymmetry, for example in the case of a curved path3 or a split-belt treadmill4, or in the case of hemiplegic stroke patients9. In addition, the relative phase between the motion of the left and right legs during walking in elderly people and patients with Parkinson’s disease fluctuates significantly5–8. It has been suggested that, unlike young, healthy people, elderly people and patients with Parkinson’s disease do exercise strict control in quiet standing24. We plan to investigate how the forgoing of control in interlimb coordination during walking changes in elderly people and patients with Parkinson’s disease. We believe that our method will contribute to the understanding of not only motor control strategies employed by young, healthy humans, but also motor disorder mechanisms in elderly people and patients with neurological disabilities. Finally, we point out that our method derives results using data measured when the relative phase between the legs recovers after being disturbed by intermittent perturbations of short duration. Although we confirmed that reasonable results were obtained from relatively small amount of measured data in young, healthy humans (see Supplementary Note 6), more data would be necessary in elderly people and patients with neurological disabilities. At present, our approach has the limitation that it does not allow us to increase the number of perturbations within a period and force subjects to walk for a long time. In the near future, we will improve our method to reduce the burden on subjects and develop applications for the diagnosis and treatment of gait irregularities.

Methods

Measurement

To extract the phase coupling functions that govern the control of interlimb coordination as modeled by the phase equation, we use kinematic data measured during walking of eight subjects. However, we frequently encountered difficulty in extracting these functions in the case of steady-state walking under normal conditions, because in such situations, the left and right legs are highly synchronized in the antiphase relationship, and for this reason, the measured data provide little information regarding the interaction between the legs. To overcome this difficulty, we varied the belt speed of the treadmill on which the subjects walked in order to disrupt the antiphase relationship between the legs. In response to such perturbations, the control of interlimb coordination was exercised, and this allowed us to extract the phase coupling functions from the measured data. This study was approved by the Ethics Committee of Doshisha University. Written informed consent was obtained from all subjects after the procedures had been fully explained. All ethical regulations relevant to human research participants were followed. A part of the measured data was used in refs. 48,49 for purposes different from that of the present study.

The subjects walked on a treadmill (ITR3017, BERTEC corporation) with a belt speed of 1.0 m/s. We used a motion capture system (MAC3D Digital RealTime System; NAC Image Technology, Inc.) to measure the kinematic data with a sampling frequency of 500 Hz. Reflective markers were attached to the subjects’ skin at several easily identifiable positions on both the left and right sides: the head, the top of the acromion, the greater trochanter, the lateral condyle of the knee, the lateral malleolus, the head of the second metatarsal, and the heel. We used the markers at the greater trochanter of the hip and the head of the second metatarsal of the foot to calculate the orientations of the limb axis φL(t) and φR(t) in Fig. 1a. We measured the belt speed of the treadmill using a rotary encoder with a sampling frequency of 500 Hz.

We applied short duration intermittent perturbations to the walking behavior by suddenly changing the belt speed. We used two types of perturbations: temporary acceleration and temporary deceleration. The temporary acceleration and deceleration perturbations increased and decreased the belt speed, respectively, by 0.6 m/s over a period of 0.1 s. The belt speed then returned to the original speed of 1.0 m/s in the following 0.1 s. The perturbations began approximately 10 s after the start of the measurement process and were applied approximately every 5 s thereafter. Each trial included ten perturbations and lasted approximately 60 s. We used three types of perturbation conditions for each trial: acceleration, deceleration, and mixed conditions. Under the first two conditions, only temporary acceleration and deceleration perturbations, respectively, were applied, while under the mixed condition, both types of perturbations were applied randomly.

All subjects (Subjects A–H) were healthy men (n = 8; age 21–23 years; weight 50.2–79.5 kg; height 161–182 cm). Although we recruited participants regardless of gender, no women joined, which would be a drawback. Subjects A and B performed 25 trials under both acceleration and deceleration conditions (with the exception that Subject A performed 24 trials under the deceleration condition). The other six subjects (Subjects C–H) performed 15 trials under each of the perturbation conditions. Each trial consisted of between 805 and 1257 walking cycles in the case of acceleration conditions, between 781 and 1248 walking cycles in the case of deceleration conditions, and between 787 and 852 walking cycles in the case of mixed conditions.

Data processing

In preparation for fitting the various quantities appearing in the phase equation with the measured time-series data, we first transformed the time series of the kinematic angle φi(t) to those of the oscillator phase, ϕi(t). (We use the subscript i ∈ (L, R), which indicates the left or right leg, instead of L and R in the subsequent sections.) In addition, we transformed the time series of the external perturbation p(t) to those for the rescaled external perturbation I(t).

Oscillator phase

First, we obtained the time series for the protophase θi(t) from the time series for the kinematic angle φi(t) using the Hilbert transform50 of φi(t), φiH(t), as4 Ai(t)expiθi(t)=φi(t)+iφiH(t),

where Ai(t) is the amplitude. The phase of an uncoupled limit-cycle oscillator is generally defined in phase reduction theory as always increasing at a constant natural frequency (ϕ°=ω for phase ϕ and natural frequency ω)15–18. However, the protophase θi(t) in Eq. (4) does not exhibit such a linear time evolution and is not suitable for the phase equation. To statistically rectify this problem, we transformed θi(t) to ϕi(t) as described in refs. 19,51 such that it tends to evolve linearly in time with the natural frequency through the definition5 ϕi(θi)=2π∫0θigi(θ′)dθ′,

where gi(θi) is the probability density function of θi(t) calculated from the time-series data for all trials of the walking experiments for each subject and each perturbation condition. We used ϕi(t) as the oscillator phase.

External perturbation

To make the phase sensitivity function Zi a dimensionless quantity, we introduced I(t) through the definition I(t)=2πp(t)/(v¯T¯), where v¯ is the average treadmill belt speed, T¯ is the average gait cycle duration, and v¯T¯ is the average stride length for one gait cycle.

Outlier exclusion

Unexpected events during walking, such as stumbling, yields outliers in the derivative of the time series of the oscillator phase. This reduces the accuracy of the fitting. To address this problem, we evaluated the derivative of the oscillator phase at time tτ, viτ, as follows:6 viτ=ϕi(tτ+1)−ϕi(tτ)tτ+1−tτ.

We performed the Smirnov-Grubbs test on the sets of vLτ and vRτ to determine the times at which outliers occur (p < 0.05) from the time-series data of all trials. In this procedure, we excluded time series during which external perturbations were applied (green regions in Fig. 2) from consideration. If an outlier appeared at some time tτ, the data for ϕi(tτ) and I(tτ) were excluded from the original data set. Outliers were often detected during touch-down events, when the belt speed exhibited small, sudden variations (Fig. 2).

Fitting the phase equations

Using the time series of the oscillator phase ϕi(t) and external perturbation I(t), we fit the parameters in the phase equations (Eqs. (1) and (2)) for each subject and perturbation condition. (We obtained 22 sets of results, corresponding to two perturbation conditions for two subjects and three perturbation conditions for six subjects.) Specifically, we obtained values for the natural frequency ωi, the phase coupling function Γij, and the phase sensitivity function Zi on the basis of the Bayesian inference method14,17,21,22.

Phase coupling function

Because the phase coupling function Γij(ψ) is 2π-periodic based on the phase reduction theory15–18, it can be represented by a Fourier series as7 Γij(ψ)=aij(0)+∑m=1Mijaij(m)cosmψ+bij(m)sinmψ,

where aij(m) (m = 0, 1, …, Mij) and bij(m) (m = 1, 2, …, Mij) are coefficients that determine the shape of the phase coupling function, and Mij is the maximum order of the Fourier series. When Mij = 0, only the constant term remains, and we have Γij(ψ)=aij(0). An appropriate value is necessary for Mij, and it is determined through a model selection, as explained below. To avoid redundancy of the two constant terms ωi and aij(0) in Eqs. (1), (2), and (7), we define ω^i:=ωi+aij(0) and Γ^ij(ψ):=Γij(ψ)−aij(0).

Phase sensitivity function

Because the phase sensitivity function Zi(ψ) is also 2π-periodic based on the phase reduction theory15–18, it can be represented by a Fourier series as8 Zi(ψ)=ci(0)+∑m=15ci(m)cosmψ+di(m)sinmψ,

where cij(m) (m = 0, 1, …, 5) and dij(m) (m = 1, 2, …, 5) are coefficients that determine the shape of the phase sensitivity function. Because the phase sensitivity function is often composed of lower-order terms, we used a maximum order of 5. We confirmed that the phase sensitivity functions obtained with this restriction were almost identical to those obtained using higher maximum orders up to 10.

Bayesian inference method

Substituting Eqs. (7) and (8) into the phase equations (1) and (2), we obtain9 ϕ°i(t)=ω^i+∑m=1Mijaij(m)cosmΔij+bij(m)sinmΔij+ci(0)+∑m=15ci(m)cosmϕi+di(m)sinmϕiI(t)+ξi(t),

where Δij = ϕj − ϕi ((i, j) = (L, R), (R, L)). Then, using the time series of the oscillator phase {ϕi(tτ)}i,τ and those of the external perturbation {I(tτ)}τ (τ = 0, 1, …, T), where T is the sample number for each subject and perturbation condition, we obtain the unknown parameters ω^i, {aij(m),bij(m)}i,j,m, {ci(m),di(m)}i,m, and D^i (:=σi2Δt), where Δt is the sampling interval. For simplicity, we define the following:wi:=ω^i,aijT,ciTT,aij:=aij(1),bij(1),aij(2),bij(2),…aij(Mij),bij(Mij)T,ci:=ci(0),ci(1),di(1),ci(2),di(2),…ci(5),di(5)T.

Here, wi, whose dimension is Ni (: = 12 + 2Mij), represents all of the unknown parameters of one phase equation, except for D^i, and aij and ci represent the unknown parameters of the phase coupling function, Γ^ij, and phase sensitivity function, Zi, respectively.

Following the Bayesian estimation method13,14, we define the likelihood function using a Gaussian distribution N as follows:10 P{ϕi(tτ)}i,τwi,D^i,{I(tτ)}τ=∏τ=0T−1Nviτω^i+Zi(ϕi(tτ))I(tτ)+Γ^ij(Δij(tτ)),D^i.

This function represents the probability to reproduce the time series {ϕi(tτ)}i,τ when the parameters wi and D^i and the time series {I(tτ)}τ are given. We adopt a Gaussian-inverse-gamma distribution IG for the conjugate prior distribution as follows:11 Ppriorwi,D^i=Nwiχiold,D^iΣioldIGD^iαiold,βiold.

Here, the hyperparameters χold and Σold determine the mean and covariance matrices, respectively, of the Gaussian distribution N, and the other hyperparameters, αold and βold, determine the shape and scale, respectively, of the Gaussian-inverse-gamma distribution IG, which is given byIGx∣α,β=βαΓ(α)x−α−1exp−βx,

where Γ(⋅) is the gamma function. In particular, the first element of χiold controls ω^i, and the other elements control aij and ci. In this study, the hyperparameters in the prior distribution were set asχioldk:=⟨viτ⟩τk=10k≠1,Σiold:=0.1Ei,αiold=βiold=2,

where Ei is the Ni × Ni identity matrix. According to Bayes’ theorem, we obtain the posterior distribution of wi and D^i from the product of the likelihood function and prior distribution:12 Ppostwi,D^i{ϕi(tτ)}i,τ,{I(tτ)}τ∝P{ϕi(tτ)}i,τwi,D^i,{I(tτ)}τ×Ppriorwi,D^i.

Update rule for hyperparameters

Because of the conjugacy of the prior distribution (Eq. (11)) to the likelihood function (Eq. (10)), the posterior distribution (Eq. (12)) is obtained as a Gaussian-inverse-gamma distribution IG as follows:Ppost(wi,D^i)=Nwiχinew,D^iΣinewIGD^iαinew,βinew.

Therefore, it is characterized by the hyperparameters χinew, Σinew, αinew, and βinew in the same way as the prior distribution (Eq. (11)).

The update rule for the hyperparameters is derived in a simple way in accordance with the Bayesian update rules. First, we define the matrix Fi∈RT×Ni and the vector vi∈RT as follows:Fi:=1gi0Thij0T⋮⋮⋮1giT−1ThijT−1T,vi:=vi0,vi1,…,viT−1T,

where the vectors giτ and hijτ are given bygiτ=I(tτ)1,cosϕi(tτ),sinϕi(tτ),cos2ϕi(tτ),sin2ϕi(tτ),…,cos5ϕi(tτ),sin5ϕi(tτ)T,hijτ=cosΔij(tτ),sinΔij(tτ),cos2Δij(tτ),sin2Δij(tτ),…,cosMijΔij(tτ),sinMijΔij(tτ)T.

Then, the hyperparameters in the posterior distribution are calculated as follows:χinew=ΣinewFiTvi+(Σiold)−1χiold,Σinew=Σiold−1+FiTFi−1,αinew=αiold+T2,βinew=βiold+12viTvi+χioldTΣiold−1χiold−χinewTΣinew−1χinew.

Finally, the unknown parameters wi and D^i of the phase equation (Eq. (9)) are determined from the mean of the posterior distribution.

Model selection

For the phase coupling function (Eq. (7)), we need to choose the model complexity, i.e., an optimal maximum order Mij of the Fourier series. For example, if Mij is too small, there will not be sufficiently many terms to correctly represent the function, whereas if Mij is too large, there will unnecessarily be higher-order harmonic terms involved in representing the function, and this will lead to overfitting. Following standard Bayesian methods13, we determined an optimal value for Mij under the assumption that this value maximizes the model evidence PME for 0≤Mij≤20 as follows:13 Mij=argmax0≤Mγν′≤20PME{ϕi(tτ)}i,τMγν′,

14 PME{ϕi(tτ)}i,τMγν′=P{ϕi(tτ)}i,τwi,D^i,{I(tτ)}τPpriorwi,D^iPpostwi,D^i{ϕi(tτ)}i,τ,{I(tτ)}τ.

The values of ΔLR are distributed within a narrow region near π, and the description provided by the phase coupling function is limited to this region. The lack of data outside this region results in the overfitting. Using this method, the optimal value for Mij was determined in each case to be in the range 11–20.

Statistics and reproducibility

We used one-tailed t-test to confirm the flat region around the antiphase relationship and the steep negative slopes outside the flat region in the control of interlimb coordination by rejecting that the means of the flat region and the slopes were zero, respectively (n = 22; p ≪ 0.01). The data points for the t-test were obtained by approximating the estimated function of the control of interlimb coordination with a piecewise linear function for each subject and perturbation condition.

Supplementary information

Peer Review File

Supplementary Information

Description of Additional Supplementary File

Supplementary Data 1-3

Supplementary information

The online version contains supplementary material available at 10.1038/s42003-024-06843-w.

Acknowledgements

This work was supported in part by the following: JSPS KAKENHI Grant Numbers JP20K21810, JP20H04144, JP20K20520, and JP24H00723; JST CREST Grant Number JPMJCR09U2; JST FOREST Program Grant Number JPMJFR2021; and MEXT KAKENHI Grant Number JP23H04467.

Author contributions

T. Aoyagi and S. Aoi developed the study design. T. Funato performed human experiments. T. Arai and K. Ota analyzed the data in consultation with K. Tsuchiya, T. Aoyagi, and S. Aoi. T. Arai and S. Aoi wrote the manuscript, and T. Aoyagi revised it. All authors reviewed and approved the final version of the manuscript.

Peer review

Peer review information

Communications Biology, thanks Arthur Dewolf, Shai Revzen and Satyajit Ambikefor their contribution to the peer review of this work. Primary Handling Editor: Joao Valente. A peer review file is available.

Data availability

Source data for graphs presented in the main figures can be found in Supplementary Data. Other data that support the findings of this study are available from the authors upon reasonable request.

Code availability

A sample program (Python code) used in this study is available from the authors upon reasonable request.

Competing interests

The authors declare no competing interests.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

These authors contributed equally: Takahiro Arai, Kaiichiro Ota.
==== Refs
References

1. Alexander, R.McN. Principles of animal locomotion. Princeton University Press: Princeton, NJ (2003).
2. Haddad JM van Emmerik REA Whittlesey SN Hamill J Adaptations in interlimb and intralimb coordination to asymmetrical loading in human walking Gait Posture 2006 23 429 434 10.1016/j.gaitpost.2005.05.006 16099160
Haddad, J. M., van Emmerik, R. E. A., Whittlesey, S. N. & Hamill, J. Adaptations in interlimb and intralimb coordination to asymmetrical loading in human walking. Gait Posture 23, 429–434 (2006).16099160
3. Courtine G Schieppati M Human walking along a curved path. II. Gait features and EMG patterns Eur. J. Neurosci. 2003 18 191 205 10.1046/j.1460-9568.2003.02737.x 12859352
Courtine, G. & Schieppati, M. Human walking along a curved path. II. Gait features and EMG patterns. Eur. J. Neurosci. 18, 191–205 (2003).12859352
4. Reisman DS Block HJ Bastian AJ Interlimb coordination during locomotion: What can be adapted and stored? J. Neurophysiol. 2005 94 2403 2415 10.1152/jn.00089.2005 15958603
Reisman, D. S., Block, H. J. & Bastian, A. J. Interlimb coordination during locomotion: What can be adapted and stored? J. Neurophysiol. 94, 2403–2415 (2005).15958603
5. Plotnik M Giladi N Hausdorff JM A new measure for quantifying the bilateral coordination of human gait: Effects of aging and Parkinson’s disease Exp. Brain. Res. 2007 181 561 570 10.1007/s00221-007-0955-7 17503027
Plotnik, M., Giladi, N. & Hausdorff, J. M. A new measure for quantifying the bilateral coordination of human gait: Effects of aging and Parkinson’s disease. Exp. Brain. Res. 181, 561–570 (2007).17503027
6. Plotnik M Giladi N Hausdorff JM Bilateral coordination of walking and freezing of gait in Parkinson’s disease Eur. J. Neurosci. 2008 27 1999 2006 10.1111/j.1460-9568.2008.06167.x 18412621
Plotnik, M., Giladi, N. & Hausdorff, J. M. Bilateral coordination of walking and freezing of gait in Parkinson’s disease. Eur. J. Neurosci. 27, 1999–2006 (2008).18412621
7. Plotnik, M. & Hausdorff, J. M. The role of gait rhythmicity and bilateral coordination of stepping in the pathophysiology of freezing of gait in Parkinson’s disease. Mov. Disord. 23, S444–450 (2008).
8. Tanahashi T Noisy interlimb coordination can be a main cause of freezing of gait in patients with little to no Parkinsonism PLoS One 2013 8 e84423 10.1371/journal.pone.0084423 24391953
Tanahashi, T. et al. Noisy interlimb coordination can be a main cause of freezing of gait in patients with little to no Parkinsonism. PLoS One 8, e84423 (2013).24391953
9. Tseng S-C Morton SM Impaired interlimb coordination of voluntary leg movements in poststroke hemiparesis J. Neurophysiol. 2010 104 248 257 10.1152/jn.00906.2009 20463199
Tseng, S.-C. & Morton, S. M. Impaired interlimb coordination of voluntary leg movements in poststroke hemiparesis. J. Neurophysiol. 104, 248–257 (2010).20463199
10. Aach M Voluntary driven exoskeleton as a new tool for rehabilitation in chronic spinal cord injury: a pilot study Spine J. 2014 14 2847 2853 10.1016/j.spinee.2014.03.042 24704677
Aach, M. et al. Voluntary driven exoskeleton as a new tool for rehabilitation in chronic spinal cord injury: a pilot study. Spine J. 14, 2847–2853 (2014).24704677
11. Bryan GM Franks PW Klein SC Peuchen RJ Collins SH A hip-knee-ankle exoskeleton emulator for studying gait assistance Int. J. Robot. Res. 2021 40 722 746 10.1177/0278364920961452
Bryan, G. M., Franks, P. W., Klein, S. C., Peuchen, R. J. & Collins, S. H. A hip-knee-ankle exoskeleton emulator for studying gait assistance. Int. J. Robot. Res. 40, 722–746 (2021).
12. Grimmer M Riener R Walsh CJ Seyfarth A Mobility related physical and functional losses due to aging and disease - a motivation for lower limb exoskeletons J. NeuroEng. Rehabil. 2019 16 2 10.1186/s12984-018-0458-8 30606194
Grimmer, M., Riener, R., Walsh, C. J. & Seyfarth, A. Mobility related physical and functional losses due to aging and disease - a motivation for lower limb exoskeletons. J. NeuroEng. Rehabil. 16, 2 (2019).30606194
13. Bishop, C.M. Pattern recognition and machine learning. Springer-Verlag: New York (2006).
14. Ota, K. and Aoyagi, T. Direct extraction of phase dynamics from fluctuating rhythmic data based on a Bayesian approach. arXiv:1405.4126 (2014).
15. Ermentrout, G.B. and Terman, D.H. Mathematical foundations of neuroscience. Springer-Verlag: New York (2010).
16. Hoppensteadt, F.C. and Izhikevich, E.M. Weakly connected neural networks. Springer-Verlag: New York (1997).
17. Kuramoto, Y. Chemical oscillations, waves, and turbulences. Springer-Verlag: Berlin (1984).
18. Winfree, A.T. The geometry of biological time. Springer-Verlag: New York (1980).
19. Kralemann B Cimponeriu L Rosenblum M Pikovsky A Mrowka R Phase dynamics of coupled oscillators reconstructed from data Phys. Rev. E 2008 77 066205 10.1103/PhysRevE.77.066205
Kralemann, B., Cimponeriu, L., Rosenblum, M., Pikovsky, A. & Mrowka, R. Phase dynamics of coupled oscillators reconstructed from data. Phys. Rev. E 77, 066205 (2008).
20. Revzen S Guckenheimer JM Estimating the phase of synchronized oscillators Phys. Rev. E 2008 78 051907 10.1103/PhysRevE.78.051907
Revzen, S. & Guckenheimer, J. M. Estimating the phase of synchronized oscillators. Phys. Rev. E 78, 051907 (2008).
21. Onojima T Goto T Mizuhara H Aoyagi T A dynamical systems approach for estimating phase interactions between rhythms of different frequencies from experimental data PLoS Comput. Biol. 2018 14 e1005928 10.1371/journal.pcbi.1005928 29337999
Onojima, T., Goto, T., Mizuhara, H. & Aoyagi, T. A dynamical systems approach for estimating phase interactions between rhythms of different frequencies from experimental data. PLoS Comput. Biol. 14, e1005928 (2018).29337999
22. Ota K Aihara I Aoyagi T Interaction mechanisms quantified from dynamical features of frog choruses R. Soc. Open Sci. 2020 7 191693 10.1098/rsos.191693 32269798
Ota, K., Aihara, I. & Aoyagi, T. Interaction mechanisms quantified from dynamical features of frog choruses. R. Soc. Open Sci. 7, 191693 (2020).32269798
23. Couzin-Fuchs E Kiemel T Gal O Ayali A Holmes P Intersegmental coupling and recovery from perturbations in freely running cockroaches J. Exp. Biol. 2015 218 285 297 10.1242/jeb.112805 25609786
Couzin-Fuchs, E., Kiemel, T., Gal, O., Ayali, A. & Holmes, P. Intersegmental coupling and recovery from perturbations in freely running cockroaches. J. Exp. Biol. 218, 285–297 (2015).25609786
24. Suzuki Y Postural instability via a loss of intermittent control in elderly and patients with Parkinson’s disease: A model-based and data-driven approach Chaos 2020 30 113140 10.1063/5.0022319 33261318
Suzuki, Y. et al. Postural instability via a loss of intermittent control in elderly and patients with Parkinson’s disease: A model-based and data-driven approach. Chaos 30, 113140 (2020).33261318
25. Asai Y A model of postural control in quiet standing: Robust compensation of delay-induced instability using intermittent activation of feedback control PLoS One 2009 4 e6169 10.1371/journal.pone.0006169 19584944
Asai, Y. et al. A model of postural control in quiet standing: Robust compensation of delay-induced instability using intermittent activation of feedback control. PLoS One 4, e6169 (2009).19584944
26. Bottaro A Yasutake Y Nomura T Casadio M Morasso P Bounded stability of the quiet standing posture: An intermittent control model Hum. Mov. Sci. 2008 27 473 495 10.1016/j.humov.2007.11.005 18342382
Bottaro, A., Yasutake, Y., Nomura, T., Casadio, M. & Morasso, P. Bounded stability of the quiet standing posture: An intermittent control model. Hum. Mov. Sci. 27, 473–495 (2008).18342382
27. Nomura T Oshikawa S Suzuki Y Kiyono K Morasso P Modeling human postural sway using an intermittent control and hemodynamic perturbations Math. Biosci. 2013 245 86 95 10.1016/j.mbs.2013.02.002 23435118
Nomura, T., Oshikawa, S., Suzuki, Y., Kiyono, K. & Morasso, P. Modeling human postural sway using an intermittent control and hemodynamic perturbations. Math. Biosci. 245, 86–95 (2013).23435118
28. Suzuki Y Nomura T Casadio M Morasso P Intermittent control with ankle, hip, and mixed strategies during quiet standing: A theoretical proposal based on a double inverted pendulum model J. Theor. Biol. 2012 310 55 79 10.1016/j.jtbi.2012.06.019 22732276
Suzuki, Y., Nomura, T., Casadio, M. & Morasso, P. Intermittent control with ankle, hip, and mixed strategies during quiet standing: A theoretical proposal based on a double inverted pendulum model. J. Theor. Biol. 310, 55–79 (2012).22732276
29. Borghese NA Bianchi L Lacquaniti F Kinematic determinants of human locomotion J. Physiol. 1996 494 863 879 10.1113/jphysiol.1996.sp021539 8865081
Borghese, N. A., Bianchi, L. & Lacquaniti, F. Kinematic determinants of human locomotion. J. Physiol. 494, 863–879 (1996).8865081
30. Bianchi L Angelini D Orani GP Lacquaniti F Kinematic coordination in human gait: relation to mechanical energy cost J. Neurophysiol. 1998 79 2155 2170 10.1152/jn.1998.79.4.2155 9535975
Bianchi, L., Angelini, D., Orani, G. P. & Lacquaniti, F. Kinematic coordination in human gait: relation to mechanical energy cost. J. Neurophysiol. 79, 2155–2170 (1998).9535975
31. Ivanenko YP Cappellini G Dominici N Poppele RE Lacquaniti F Modular control of limb movements during human locomotion J. Neurosci. 2007 27 11149 11161 10.1523/JNEUROSCI.2644-07.2007 17928457
Ivanenko, Y. P., Cappellini, G., Dominici, N., Poppele, R. E. & Lacquaniti, F. Modular control of limb movements during human locomotion. J. Neurosci. 27, 11149–11161 (2007).17928457
32. Ivanenko YP d’Avella A Poppele RE Lacquaniti F On the origin of planar covariation of elevation angles during human locomotion J. Neurophysiol. 2008 99 1890 1898 10.1152/jn.01308.2007 18272871
Ivanenko, Y. P., d’Avella, A., Poppele, R. E. & Lacquaniti, F. On the origin of planar covariation of elevation angles during human locomotion. J. Neurophysiol. 99, 1890–1898 (2008).18272871
33. Catavitello G Ivanenko Y Lacquaniti F A kinematic synergy for terrestrial locomotion shared by mammals and birds eLife 2018 7 e38190 10.7554/eLife.38190 30376448
Catavitello, G., Ivanenko, Y. & Lacquaniti, F. A kinematic synergy for terrestrial locomotion shared by mammals and birds. eLife 7, e38190 (2018).30376448
34. Bosco G Poppele RE Proprioception from a spinocerebellar perspective Physiol. Rev. 2001 81 539 568 10.1152/physrev.2001.81.2.539 11274339
Bosco, G. & Poppele, R. E. Proprioception from a spinocerebellar perspective. Physiol. Rev. 81, 539–568 (2001).11274339
35. Poppele RE Bosco G Rankin AM Independent representations of limb axis length and orientation in spinocerebellar response components J. Neurophysiol. 2002 87 409 422 10.1152/jn.00022.2001 11784759
Poppele, R. E., Bosco, G. & Rankin, A. M. Independent representations of limb axis length and orientation in spinocerebellar response components. J. Neurophysiol. 87, 409–422 (2002).11784759
36. Wilshin, S., Kvalheim, M.D., Scott, C., and Revzen, S. Estimating phase from observed trajectories using the temporal 1-form. arXiv:2203.04498 (2022).
37. Hägglund M Optogenetic dissection reveals multiple rhythmogenic modules underlying locomotion Proc. Natl. Acad. Sci. USA 2013 110 11589 11594 10.1073/pnas.1304365110 23798384
Hägglund, M. et al. Optogenetic dissection reveals multiple rhythmogenic modules underlying locomotion. Proc. Natl. Acad. Sci. USA 110, 11589–11594 (2013).23798384
38. Kudo N Yamada T Morphological and physiological studies of development of the monosynaptic reflex pathway in the rat lumbar spinal cord J. Physiol. 1987 389 441 459 10.1113/jphysiol.1987.sp016665 2824763
Kudo, N. & Yamada, T. Morphological and physiological studies of development of the monosynaptic reflex pathway in the rat lumbar spinal cord. J. Physiol. 389, 441–459 (1987).2824763
39. Soffe SR Roles of glycinergic inhibition and N-methyl-D-aspartate receptor mediated excitation in the locomotor rhythmicity of one half of the Xenopus embryo central nervous system Eur. J. Neurosci. 1989 1 561 571 10.1111/j.1460-9568.1989.tb00363.x 12106115
Soffe, S. R. Roles of glycinergic inhibition and N-methyl-D-aspartate receptor mediated excitation in the locomotor rhythmicity of one half of the Xenopus embryo central nervous system. Eur. J. Neurosci. 1, 561–571 (1989).12106115
40. Dietz V Zijlstra W Duysens J Human neuronal interlimb coordination during split-belt locomotion Exp. Brain Res. 1994 101 513 520 10.1007/BF00227344 7851518
Dietz, V., Zijlstra, W. & Duysens, J. Human neuronal interlimb coordination during split-belt locomotion. Exp. Brain Res. 101, 513–520 (1994).7851518
41. Forssberg H Grillner S Halbertsma J Rossignol S The locomotion of the low spinal cat. II. Interlimb coordination Acta. Physiol. Scand. 1980 108 283 295 10.1111/j.1748-1716.1980.tb06534.x 7376923
Forssberg, H., Grillner, S., Halbertsma, J. & Rossignol, S. The locomotion of the low spinal cat. II. Interlimb coordination. Acta. Physiol. Scand. 108, 283–295 (1980).7376923
42. Prokop T Berger W Zijlstra W Dietz V Adaptational and learning processes during human split-belt locomotion: interaction between central mechanisms and afferent input Exp. Brain Res. 1995 106 449 456 10.1007/BF00231067 8983988
Prokop, T., Berger, W., Zijlstra, W. & Dietz, V. Adaptational and learning processes during human split-belt locomotion: interaction between central mechanisms and afferent input. Exp. Brain Res. 106, 449–456 (1995).8983988
43. Yang JF Lamont EV Pang MYC Split-belt treadmill stepping in infants suggests autonomous pattern generators for the left and right leg in humans J. Neurosci. 2005 25 6869 6876 10.1523/JNEUROSCI.1765-05.2005 16033896
Yang, J. F., Lamont, E. V. & Pang, M. Y. C. Split-belt treadmill stepping in infants suggests autonomous pattern generators for the left and right leg in humans. J. Neurosci. 25, 6869–6876 (2005).16033896
44. Danner SM Spinal V3 interneurons and left-right coordination in mammalian locomotion Front. Cell. Neurosci. 2019 13 516 10.3389/fncel.2019.00516 31824266
Danner, S. M. et al. Spinal V3 interneurons and left-right coordination in mammalian locomotion. Front. Cell. Neurosci. 13, 516 (2019).31824266
45. Laflamme OD Distinct roles of spinal commissural interneurons in transmission of contralateral sensory information Curr. Biol. 2023 33 1 13 10.1016/j.cub.2023.07.014 36446352
Laflamme, O. D. et al. Distinct roles of spinal commissural interneurons in transmission of contralateral sensory information. Curr. Biol. 33, 1–13 (2023).36446352
46. Molkov YI Bacak BJ Talpalar AE Rybak IA Mechanisms of left-right coordination in mammalian locomotor pattern generation circuits: A mathematical modeling view PLoS Comput. Biol. 2015 11 e1004270 10.1371/journal.pcbi.1004270 25970489
Molkov, Y. I., Bacak, B. J., Talpalar, A. E. & Rybak, I. A. Mechanisms of left-right coordination in mammalian locomotor pattern generation circuits: A mathematical modeling view. PLoS Comput. Biol. 11, e1004270 (2015).25970489
47. Zhang H The role of V3 neurons in speed-dependent interlimb coordination during locomotion in mice eLife 2022 11 e73424 10.7554/eLife.73424 35476640
Zhang, H. et al. The role of V3 neurons in speed-dependent interlimb coordination during locomotion in mice. eLife 11, e73424 (2022).35476640
48. Funato T Aoi S Tomita N Tsuchiya K Validating the feedback control of intersegmental coordination by fluctuation analysis of disturbed walking Exp. Brain Res. 2015 233 1421 1432 10.1007/s00221-015-4216-x 25657068
Funato, T., Aoi, S., Tomita, N. & Tsuchiya, K. Validating the feedback control of intersegmental coordination by fluctuation analysis of disturbed walking. Exp. Brain Res. 233, 1421–1432 (2015).25657068
49. Funato T Evaluation of the phase-dependent rhythm control of human walking using phase response curves PLoS Compt. Biol. 2016 12 e1004950 10.1371/journal.pcbi.1004950
Funato, T. et al. Evaluation of the phase-dependent rhythm control of human walking using phase response curves. PLoS Compt. Biol. 12, e1004950 (2016).
50. Pikovsky, A., Rosenblum, M., and Kurths, J. Synchronization: A universal concept in nonlinear sciences. Cambridge University Press: Cambridge (2001).
51. Kralemann B Cimponeriu L Rosenblum M Pikovsky A Mrowka R Uncovering interaction of coupled oscillators from data Phys. Rev. E 2007 76 055201 10.1103/PhysRevE.76.055201
Kralemann, B., Cimponeriu, L., Rosenblum, M., Pikovsky, A. & Mrowka, R. Uncovering interaction of coupled oscillators from data. Phys. Rev. E 76, 055201 (2007).
