A spinal network of proprioceptive reflexes can produce a variety of bipedal gaits

TL;DR

A simple spinal reflex system can generate various natural bipedal gaits, suggesting proprioceptive reflexes play a key role in movement control.

Contribution

A fixed-gain reflex controller is shown to produce diverse gaits without rhythmic inputs or state machines.

Findings

The reflex controller generates walking, hopping, and running gaits in multiple directions and speeds.

Natural gaits emerge without high-level control or modulated feedback gains.

The model suggests proprioceptive reflexes are more flexible and essential for locomotion than previously thought.

Abstract

Proprioception is crucial for movement, yet the role of proprioceptive reflexes in legged locomotion is still poorly understood. While previous simulation studies have shown great potential for reflex-based control strategies, these controllers are typically geared to specific gaits, using hand-crafted feedback pathways that are linked to specific gait phases. In this work, we explore the control capabilities of a simple reflex controller that consists of only homonymous and antagonistic length and force feedback pathways with constant gains. This control model can be considered a highly simplified subset of spinal control rather than an attempt to emulate all spinal control functions. Despite its simplicity, we found our control framework capable of producing a wide variety of natural gaits, including walking and hopping, forwards and backwards, and running in different variations and at different velocities – without requiring any rhythmic inputs or high-level state machines modulating the feedback gains. Our work highlights the important role and flexibility of proprioceptive reflexes and suggests a necessary re-evaluation of their role in locomotion. Due to its simplicity and flexibility, our control framework provides a solid basis for the development of higher-level neuromuscular control systems.

A spinal reflex controller with fixed-gain muscle force and length feedback can generate a wide range of natural gaits, offering insight into how proprioceptive reflexes may contribute to locomotor control.

Figures

Overview.Our reflex controller can produce a variety of gaits in a human musculoskeletal model. A The planar musculoskeletal model consists of seven segments connected by 6 joints, which are actuated by 18 Hill-type muscles. B Our reflex controller calculates stimulation for each muscle *i* from delayed force $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{F}$$\end{document}$ and length $\documentclass[12pt

Forwards walking.Joint kinematics (first row) and corresponding muscular activation for hip (left), knee (middle), and ankle (right) joints, as well as ground reaction force (bottom right). Data are shown for the mean of all strides of our forward walking solution (red) compared to the model of Geyer and Herr^[28](#CR28)^ (blue), experimental data averaged over all subjects from refs. ^[38](#CR38),[76](#CR76),[77](#CR77)^ (gray area, see also Section “Experimental data''), and of one height and weight matched subject AB06 from Camargo et al.^[38](#CR38)^ (dotted black).Table 2Maximum cross-cor

Running.Joint kinematics (first row) and corresponding muscular activations for hip (left), knee (middle), and ankle (right) joints, as well as ground reaction force (bottom right) for running at different velocities in comparison to experimental data of running at 3.0 m/s from Hamner and Delp^[79](#CR79)^ (gray).

Controller parameters.Parameter values for the five target gaits (left to right) for length (top row) and force (middle row) reflexes, as well as the constant prestimulations (bottom row). Matrix columns represent the source muscle the feedback is coming from, rows show the target muscle receiving a stimulation input from the reflex connection.

Control architecture.Muscle stimulations *u*_*i*_(*t*) for the nine muscles per leg are computed through delayed length ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{L}(t-\Delta t)$$\end{document}$ ) and force ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\odd

Initial states.For each target gait, the simulation starts in a different initial pose (joint angles and velocities). These poses are manually tuned to minimize the initial correction the controller has to make and to direct the optimization towards the target gait. The numeric values of all initial sets can be found in Table S2.

Full text

Introduction

The role of proprioceptive reflexes in human and animal locomotion has been subject to long-standing research and debate^1,2^. In 1906, after observing that decerebrated cats could still produce basic stepping movements, Charles Sherrington proposed that complex movements arise from sequential reflex chains linking sensory input to motor output^3–5^. In 1911, Thomas Graham Brown discovered that rhythmic locomotor activity persisted even after removing sensory input to the spinal cord, indicating the presence of intrinsic spinal circuits—later dubbed central pattern generators or CPGs—capable of producing rhythmic motor output without requiring sensory input signals^2,6,7^. Neuroscientific research often highlights the importance of CPGs^8–10^, yet their exact functioning in the human spinal cord remains unclear^2^. Information from proprioceptive sensors has been shown to be essential for movement generation^11^, the regulation of stepping^12^, and the modulation of CPG output^13,14^. However, the way in which CPGs and reflexes functionally interact remains poorly understood, and in vivo preparations are insufficient to decipher the parameters and structures involved in this interaction, particularly during dynamic locomotion^8,15^.

Predictive neuromuscular simulations are a useful alternative to in vivo preparations, since they allow complete and fine-grained control over the experimental setup and help identify the potential role of individual control primitives^1,16^. In addition, predictive neuromuscular simulations have important potential applications, such as aiding in the design of assistive devices through human-in-the-loop simulation^17–21^, predicting the outcome of medical procedures^22–24^, and predicting gait stability and perturbation responses^25–27^. However, to fulfill these promises, a control framework is needed that is both neurologically grounded and capable of replicating the versatility of human movement.

An important simulation study that highlights the potential of reflexes in locomotion is the work of Geyer and Herr^28^, who developed a locomotion control strategy primarily driven by proprioceptive reflexes. Their work inspired a large number of derivative controllers, including 3D extensions and the addition of higher-level control components^29–33^ and was the basis for several subsequent predictive simulation studies on human gait^34–36^. These controllers produce rhythmic muscle activation by switching between specific reflex pathways based on the current state of each leg (i.e., swing and stance), using hand-crafted state selection logic. Furthermore, reflex pathways are domain-specific, as they are manually picked to mimic human-level ground walking kinematics and muscular activations. This tailoring to level-ground walking and dependence on state-switching limit their ability to generalize to other gaits. Work extending the controller of Geyer and Herr^28^ to other environments or movements mainly focuses on adding a higher-level component^25,30^ (Song and Geyer^29^ also add new reflex pathways in addition to a supraspinal layer) or on further refining and tailoring the states used^32,33^. While this work has provided insights on different parts of the human motor control system, this approach may obscure the potential within reflexes themselves. An important avenue of future research is to extend the capability of these control systems and allow for more complex tasks^37^, which requires a solid spinal foundation.

The purpose of this study is to explore the locomotor capabilities of proprioceptive reflexes in the absence of a high-level state selection or modulation. Instead, we use delayed length and force reflex pathways with constant feedback gains, which are active throughout the simulation, relying solely on reciprocal innervation to produce alternating rhythmic activity (see Fig. 1). For reflex connectivity, our architecture is informed by the basic organizational principles of the mammalian spinal cord^13^ and includes reflex pathways to homonymous and antagonistic muscles. For the sake of simplicity, we have opted to exclude other sensory input modalities, such as muscle fiber velocity, cutaneous sensors, and vestibular sensors. As such, our control model can be considered a highly simplified subset of spinal control rather than a full emulation. Preliminary studies have shown that adding these input modalities did not result in noticeable improvement in the discovered gait patterns.Fig. 1. Overview.Our reflex controller can produce a variety of gaits in a human musculoskeletal model. A The planar musculoskeletal model consists of seven segments connected by 6 joints, which are actuated by 18 Hill-type muscles. B Our reflex controller calculates stimulation for each muscle i from delayed force \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{F}$$\end{document} and length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{L}$$\end{document} feedback with gains κi**j and a constant offset ci. Implemented pathways include homonymous and antagonistic pathways for length and force. Feedback gains and offsets are optimized based on high-level objectives. C The controller can produce a variety of bipedal gaits like forwards and backwards walking, and hopping as well as running.

In spite of its simplicity, we found that these basic reflexes were sufficient to produce a variety of natural bipedal gaits like walking, hopping, and running in different variations and at different velocities—only by selecting different sets of constant feedback parameters. With its versatile functionality, it can be considered a basic starting point on top of which more elaborate control strategies can be developed.

Results

For all five target gaits, we found solutions that did not fall within the maximum simulation time and enter rhythmic cycles (see Video S1 and phase plots Fig. S6). The model reaches the target velocity of 1.2 m/s for forwards and backwards walking and a maximum running velocity of 3.4 m/s. The hopping solutions generate hopping forwards and backwards at a velocity of 1.2 m/s and 1.4 m/s, respectively (compare also Table 1).Table 1. Gait velocityWalk FwdWalk BkwHop FwdHop BkwRunVelocity (m/s)1.20−1.211.22−1.393.42Stride length (m)1.11−0.830.54−0.541.94Stride duration (s)0.920.690.440.380.57Obtained velocities, stride length, and duration for the five target gaits.

Our forwards walking gait matches human kinematics well (see Fig. 2, ϕh, our: R = 0.90, ϕk, our: R = 0.89, ϕa, our: R = 0.79). Even though the knee kinematics are not matched as well as by the reflex controller of Geyer and Herr^28^(ϕk, G**H: R = 0.98), the ankle kinematics show a greater maximum cross-correlation ^28^(ϕa, G**H: R = 0.69). The mismatch in knee kinematics is mostly coming from the leg being straight during stance. The vertical ground reaction force shows the characteristic two peaks found in experimental data (albeit with shifted timings, R = 0.93). Generally, the computed time shifts Δ of the maximum cross-correlation for knee, ankle, and GRF are higher for our model than for the model of Geyer and Herr^28^, but within maximally 3.6% (Δ_GRF). All maximum cross-correlation values and time shifts are given in Table 2.Fig. 2. Forwards walking.Joint kinematics (first row) and corresponding muscular activation for hip (left), knee (middle), and ankle (right) joints, as well as ground reaction force (bottom right). Data are shown for the mean of all strides of our forward walking solution (red) compared to the model of Geyer and Herr^28^ (blue), experimental data averaged over all subjects from refs. ^38,76,77^ (gray area, see also Section “Experimental data''), and of one height and weight matched subject AB06 from Camargo et al.^38^ (dotted black).Table 2. Maximum cross-correlation R and time shift Δ for joint angles and GRFϕhϕkϕa_GRFWalk FwdR**Δ (%)R**Δ (%)R**Δ (%)R**Δ (%)exp-our0.900.00.892.10.793.40.933.6exp-GH0.930.00.980.60.691.00.990.7AB06-our0.910.00.952.80.724.10.942.7AB06-GH0.810.00.980.30.751.70.980.7Runexp-our0.910.00.932.00.25200.788.5Similarity metrics computed between the mean of experimental data (exp) or mean data of subject AB06 (for Walk Fwd only) and simulation data from our model (our) or the model of Geyer and Herr^28^ (GH, Walk Fwd only).

Also, in terms of muscular activations, the controller produces smooth muscular activations that match experimental data reasonably well. Maximum cross-correlation values can be found in Table S3. The straight knee during stance is caused by an absence of VAS activation. Among the variations of forwards walking we found, there is also one gait with a knee flexion during stance and VAS activation, as well as less shifted ground reaction forces (see Fig. S7b for kinematics and activations and Tables S3 and S4 for maximum cross-correlation results). Therefore, this mismatch is likely a result of either the lack of reflex modulation that makes it difficult to differentiate between stance and swing, or the generic cost function, or a combination of both. Note, however, that the gait is still close to human walking as humans exhibit a wide variety of gait kinematics and muscular activations (e.g., height and weight matched subject AB06 from Camargo et al.^38^ (see Fig. 2 and Table 2) also shows very little VAS activation and a relatively straight leg during stance, ϕk, A**B06−our: R = 0.95).

For our running targets, we discovered various solutions with different running velocities. Figure 3 shows the kinematics and muscular activations for running gaits at different velocities ranging from 2.72 to 3.42 m/s. While ankle kinematics are not accurately matched (ϕa: R = 0.25, Δ = 10%), hip and knee angle match experimental data well (ϕh: R = 0.91, ϕk: R = 0.93), even though, as in Walk Fwd, the knee is too straight during stance. Also in Run, ground reaction forces are shifted, and initially ground reaction forces are substantially higher in the model than in the experimental data (R = 0.78, Δ = 8.5%). Except for VAS, which is not activated at all and RF which misses activation during stance, muscular activations show a fair agreement with experimental data (0.8 < R < 0.91, see Table S5), but some muscles saturate at maximum activation (GLU, SOL, ILI, GAS, TA) probably because only speed was maximized and effort was not taken into account.Fig. 3. Running.Joint kinematics (first row) and corresponding muscular activations for hip (left), knee (middle), and ankle (right) joints, as well as ground reaction force (bottom right) for running at different velocities in comparison to experimental data of running at 3.0 m/s from Hamner and Delp^79^ (gray).

The kinematics and muscle activations for Walk Bkw, Hop Fwd, and Hop Bkw are provided in Fig. S10. In an iterative, hand-tuning process, we were furthermore able to find a wide spectrum of solutions for each target gait (see Video S2). Also, a variety of other related gaits, including skipping or backwards running, can be generated. To illustrate the versatility of the reflex controller, we compiled a collection of some other interesting gaits we found (see Video S3). While these gaits are often suboptimal, they still demonstrate a lively appearance that is difficult to quantify.

A comparison of the parameter values of the five target gaits provides some insights into the parameter space of the controller (see Fig. 4). Generally, the solutions are very diverse and differ a lot between the different gaits (see Fig. S1) as well as within gait (see Fig. S8 for parameter values of Walk Fwd 2, a second, manually optimized solution for Walk Fwd). Directly interpreting and comparing prominent reflexes is difficult, as the interference between reflexes can cancel out the influence of a reflex such that its value becomes meaningless (e.g., a large negative prestimulation will cause an inactive muscle even if other reflexes would generate an activation of this muscle).Fig. 4. Controller parameters.Parameter values for the five target gaits (left to right) for length (top row) and force (middle row) reflexes, as well as the constant prestimulations (bottom row). Matrix columns represent the source muscle the feedback is coming from, rows show the target muscle receiving a stimulation input from the reflex connection.

However, more generally, we find that the offsets ci are higher for the gaits without effort minimization (hopping and running) than for the walking gaits, where effort minimization was included. Also, homonymous feedback was not only excitatory for the length reflexes, but also for most of the force reflexes, even though this was not a constraint during the optimization. This supports the idea that positive homonymous force feedback is important for human gait^39,40^. Overall, we find that for all five target gaits, most reflexes are active, with only a few reflex gains set to an absolute value below 0.1 (Walk Fwd: 3, Walk Bkw: 5, Hop Fwd: 1, Hop Bkw: 4, Run: 2). However, as argued before and as the example of VAS during walking and running shows, not all of these reflexes have an influence as they might get canceled out by prestimulations or other reflexes. In fact, the additional walking solution Walk Fwd 2 generates forwards walking with overall smaller absolute reflex gains and 19 reflexes below 0.1 (compare also Fig. S8).

Discussion

Our results demonstrate that a simple set of reflexes with constant feedback gains is enough to produce a wide variety of bipedal gaits—without the need for a state-switching mechanism or central pattern generator. Since our controller does not depend on hand-picked connections or finite-state machines specifically tailored to walking, it can generalize to other natural-looking gaits, exposing the remarkable potential of reflex-based control. We consider the natural versatility in the absence of high-level control to be the main strength of our controller.

The absence of a state switch mechanism and the restriction to homonymous and antagonistic length and force feedback are the main differences between our control architecture and the controller of Geyer and Herr^28^. However, due to its reflexive nature, our controller shares many features with the controller of Geyer and Herr^28^ and subsequent works^29,33^. First, inter-stride variability is normally relatively low as the periodic movement is based on an interplay of sensor signals (details see also Supplementary Material Section 3). Second, the initial state of the model is important and set close to the target motion as the reflex controller is dependent on the sensory feedback to get into equilibrium. In nature, supraspinal mechanisms could drive the system to the specific initial condition where reflexes can take over completely. Third, even though these controllers are not good at tracking kinematics, desired trajectories can be used in the cost function to discover reflex parameters that lead to similar kinematics. On the other hand, kinematics emerge naturally from an interplay between body mechanics, muscles, and control, and generally, muscular activations are smooth and natural without abrupt changes.

Our controller has the ability to naturally converge to human-like activations and walking kinematics, without the use of experimental data during the optimization. Kinematically, the main difference to experimental data is the straight knee during stance. The controller proposed by Geyer and Herr^28^ better reproduces knee angle and ground reaction forces. We suspect that modulation of the VAS reflexes (which Geyer and Herr^28^ approximated through state-switching) is important to properly differentiate between passive VAS during swing and active VAS during stance.

In terms of muscular activation, our results are similar to the ones of Geyer and Herr^28^: for some muscles, our controller better tracks the experimental data, while for others, their controller is more accurate, and both models do show differences to experimental data. Generally, experimental studies show there is a wide range of muscle activity patterns both intra-subject between strides as well as inter-subjects^41,42^. Our controller captures part of this diversity and allows for different forwards walking solutions with different kinematics (e.g., knee flexion during stance), but also differences in muscle activation patterns of several muscles (GLU, VAS, HAM, BF for Walk Fwd compared to Walk Fwd 2). The full spectrum of gait variability of our reflex controller may still be further explored.

Our controller generalizes better to different gaits than the work of Geyer and Herr^28^, which has been tailored to level-ground walking and therefore is not easily extended to other gaits. Song and Geyer^29^ impressively extended the controller to allow for more complex movements like turning movements and running, which required not only additional reflexes and states, but also a supraspinal control layer. They also found the trunk lean to be an important variable to modulate running speed^43^. Our controller can generate running at different speeds (and backwards walking and hopping) without any change to the controller architecture, using only length and force reflexes. The maximum velocity of 3.42 m/s is well below the maximum velocity professional athletes can run at (average velocity for world record for marathon: 5.8 m/s or 5000 m: 6.62 m/s^44^), however, it is close to the mean preferred running speed of 3.7 m/s found in experienced runners^45^. Also, an increase of maximal muscle forces (mimicking training) could likely increase the maximum velocity of our model. We also expect the simplicity of the model and lack of arms to contribute to the slower maximum speed^46^.

All presented results are the starting point for further development and studies to overcome its current limitations. We adopted the current approach to allow studying the broad versatility of the control structure. When focusing on one specific gait, the differences to experimental data could possibly be minimized by a more complex cost function including further optimization criteria like joint pain or kinematic tracking. Kinematic tracking could also provide a way to model subject differences and analyze differences in the resulting reflexes. Even though it might seem counter-intuitive in a forward simulation, for reflex controllers, this provides a neat way to constrain the kinematics, as the reflex controller can not directly follow a trajectory and therefore the kinematics only shape the cost function.

In the same line of studying the broad versatility of the controller, we allowed for all homonymous and antagonistic reflexes, such that the potential of the controller could be explored, and no bias was introduced by selecting reflexes more specifically. This, however, limits the analysis and interpretability of the found solutions, as likely not all reflexes are needed for all gaits, and interaction effects between parameters can lead to reflexes without effect, even if the reflex gain is high. In the future, these limitations could be overcome by a detailed study of the parameter space, including a minimization of active reflexes, an analysis of the connections within as well as between gaits, and an assessment of the initial state dependency. With this, the dependency on the cost function and hyperparameters could be eliminated, providing a more continuous view on the parameter space. A minimization of reflexes could help to obtain a sparser controller structure and a more meaningful set of reflexes that can be better compared to neurophysiological data, and reduce the number of redundant or unnecessary reflexes. The focus of this work was to showcase the versatility and potential of the control structure; however, the additional solution for forwards walking we found by manual tuning and iterative optimization already demonstrates that potentially less reflexes are necessary and might lead to more natural-looking results.

We purposefully designed our controller to be as simple as possible, only including a minimal set of components and sensor signals (specifically, muscle fiber length and force). We recognize that the absence of velocity feedback may appear problematic, but results from our preliminary simulations indicated that including velocity feedback did not result in noticeable improvement of the discovered gait patterns. We believe that further study is required to understand why the effect of including velocity feedback appears limited, compared to using only length feedback from muscle spindles. One hypothesis worth exploring is that additional filtering of the velocity signal is needed, since the velocity of the contractile element in simulated Hill-type muscles can fluctuate rapidly, potentially limiting its suitability as a direct feedback signal.

Our finding that this minimal network of length and force feedback can produce a rich variety of periodic gaits is surprising and encouraging. It is, thus, a functional and versatile subset of spinal control. Still, it deliberately excludes mechanisms shown elsewhere to enhance reflexive adaptability: Rybak et al.^47^ described spinal CPG architectures where multimodal feedback (length + velocity + phase-specific inputs) enhances phase transitions and robustness in cat models. Loeb^48^ and also Prochazka and Gorassini^49^ emphasized how fusimotor drive modulates spindle gain across task conditions. Likewise, recent works have demonstrated that γ-fusimotor modulation nonlinearly shapes spindle output and that anatomical agonist-antagonist definitions may be more functionally fluid than our current rigid reciprocal inhibition network^50,51^. Incorporating these features—velocity feedback, dynamic fusimotor control, and flexible reciprocal inhibition—could significantly improve adaptability, perturbation resistance, and context-sensitive gait modulation.

Our controller can easily be extended in many ways. The inclusion of supraspinal mechanisms to control gait initiation, speed, or transitions is an important direction for future research. This could also solve the dependency on the initial state, as supraspinal mechanisms could take over until a state is reached that allows for fully reflexive control. Also, a change of gait means a complete change of parameters in our current architecture, which is not plausible in humans, even though reflexes do change, e.g., between postural control and locomotion^52^. The proposed controller provides a valuable basis for extensions to more complex and targeted, voluntary movements, as well as different terrains and perturbations demanding a higher-level control structure.

Overall, our proposed methodology is sufficient to demonstrate the power of the spinal reflex control approach without high-level state selection or modulation. To our knowledge, this is the first controller that can generate different gaits using only reflexes and no higher-level control structure, finite-state machine, or CPG. It provides a solid basis for examining the relation between different gaits and incrementally developing more complex neuromuscular controllers^37^.

Reflex controllers have also shown potential in the control of exoskeletons, prosthesis, and robots, as they allow for more realistic human-in-the-loop simulation^20,21,53^. Because our controller does not rely on gait state detection, it may be easier to incorporate in these scenarios.

Conclusion

Our results demonstrate that proprioceptive reflexes are a remarkably powerful control primitive and that state-switching with a finite state machine or a CPG is not needed to generate rhythmic bipedal gaits. Furthermore, our reflex controller allows for a striking versatility of gaits, suggesting that spinal reflexes play a more prominent role in gait than expected. Even though we know that other circuits, such as central pattern generators, play an important role in gait, their importance, relation, and dependencies might need to be re-evaluated. Our proposed controller provides a valuable basis to develop more complex neuromuscular controllers and has possible applications in the control of robots, exoskeletons, prostheses, or ortheses.

Methods

Model

Simulations were performed using a planar musculoskeletal model representing an adult human male of 74.5 kg and 1.8 m height. The model consists of 7 segments (torso and two legs with femur, tibia, and foot) and 9 degrees-of-freedom (DOFs). Hip, knee and ankle joint of each leg are actuated by 9 Hill-type musculotendon units with elastic tendons and muscle fiber damping^54^: gluteus (GLU), iliopsoas (ILI), rectus femoris (RF), hamstrings (HAM), biceps femoris short head (BF), vastus (VAS), gastrocnemius (GAS), tibialis anterior (TA), and soleus (SOL) (see also Fig. 1). Both the mass properties and muscle parameters in our model (optimal fiber length, tendon slack length, pennation at optimal fiber length, maximum isometric force) were taken from Delp et al.^55^, with the exception of the hamstring parameters, for which we used updated data from Arnold et al.^56^. We further adapted the via points of the muscle pathways to match moment arm data from Rajagopal et al.^57^, which uses wrapping geometry to match experimentally measured moment arms. We avoided the use of wrapping geometry in our model for performance reasons, but found that through via-points we could match the moment arms within the margin of error of the source material.

In cases where we changed the muscle geometry, we adjusted the tendon slack lengths in such a way that the muscle fiber lengths remained optimal at the same joint angle (or joint angles for biarticular muscles), thereby minimizing changes in the angle at which passive muscle force commences. The latter is important because our control strategy relies on force feedback, which is sensitive to the passive muscle force that is triggered after a muscle stretches past its optimal fiber length. Following that same rationale, we adjusted the tendon slack length of the ankle plantarflexors in such a way that these muscles better match passive dorsiflexion muscle force measurements found in Le Sant et al.^58^. The curves representing tendon stiffness, as well as the force-length and force-velocity relations, are taken from Millard et al.^54^, including the passive damping (β = 0.1), allowing muscle activation levels to be zero. We use a vmax of 10 m/s, similar to Delp et al.^55^ and Rajagopal et al.^57^. Muscle activation is computed using a first-order dynamics model, with an activation time constant of 0.01 s and a deactivation time constant of 0.04 s^59^.

Ground contacts are modeled using the Hunt-Crossley contact model^60^ with two contact spheres (r = 0.03, stiffness = 5 × 10^6^, dissipation = 1) per foot segment, using the friction model implemented in OpenSim^61^ (static friction = 0.9, dynamic and viscous friction = 0.6). The model is implemented in SCONE^62^ and Hyfydy^63^ and is provided as supplementary material. An overview of the model parameters is provided in Table S1.

Control

Our neural control model is based on proprioceptive feedback from muscle spindles (type Ia and II sensory fibers) and Golgi tendon organs (type Ib sensory fibers)^13,64^. Our aim is not to replicate the full neurophysiological complexity of the spinal cord or to model a full decerebrated preparation, but rather to capture only a minimally functional subset. In our model, muscle spindle receptors signal information about muscle fiber length, and Golgi tendon organ receptors signal information about muscle tension. Even though Ia feedback from muscle spindles includes both fiber length and velocity information, we were unable to see the benefits of including velocity feedback in the results of our preliminary studies (see Discussion for more details). To strike a balance between controller simplicity and capability, and to keep the number of optimization parameters to a minimum, we therefore chose to leave out fiber velocity feedback and focus only on length and force feedback. The investigation of other sensory input modalities, which are also known to be important for human gait^65,66^, is left for future work.

Our model incorporates the following simplified reflex pathways commonly observed in the mammalian spinal cord: homonymous reflex connections to the same muscle and antagonistic reflex connections to all antagonist muscles that cross the same joint^11,13,67^. For length reflexes, homonymous connections are considered monosynaptic and constrained to be excitatory, while antagonistic connections are constrained to be inhibitory, representing connections via Ia inhibitory interneurons. Ib force reflexes are allowed to be either excitatory or inhibitory. While it is well known that Ib afferents have a widespread inhibitory influence, it has been observed directly in cats that homonymous Ib afferents can initiate excitatory responses during locomotor activity^68^, which has been suggested to play an important role in human gait as well^39,40^.

For our model with 9 muscles per leg this leads to connection matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{R}}}}$$\end{document} containing 9 monosynaptic (red, Fig. 5) and 22 antagonistic connections (blue, Fig. 5). All pathways are modeled with neural delays Δt ∈ [10 ms, 35 ms] based on experimental data of H-reflex latency as well as the length of the pathway (i.e., the more distal the more delay), following van der Kruk and Geijtenbeek^69^ (for the delay values see Table S1).Fig. 5. Control architecture.Muscle stimulations ui(t) for the nine muscles per leg are computed through delayed length ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{L}(t-\Delta t)$$\end{document} ) and force ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{F}(t-\Delta t)$$\end{document} ) feedback. The connection matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{R}}}}$$\end{document} (right) contains both homonymous connections innervating the muscle based on its own sensor signal (red in L matrix) and connections between antagonistic muscles (blue in L matrix). Gains are either unconstrained (gray, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\kappa }_{ii}^{F}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\kappa }_{ij}^{F}$$\end{document} ), constrained to be excitatory (i.e., ≥0, red, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\kappa }_{ii}^{L}$$\end{document} ) or inhibitory (i.e., ≤0, blue, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\kappa }_{ij}^{L}$$\end{document} ). Including the constant thresholds ci, the reflex controller contains 71 free parameters.

Overall, the controller thus produces a stimulation u(t) of each muscle i based on constant-gain delayed length ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{L}$$\end{document} ) and force ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{F}$$\end{document} ) feedback, and a constant offset ci:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${u}_{i}(t)={c}_{i}+{\sum }_{j\in {{{{\mathcal{R}}}}}_{i}}^{}[{\kappa \, }_{ij}^{L}{\tilde{L}}_{j}(t-\Delta t)+{\kappa \, }_{ij}^{F}{\tilde{F}}_{j}(t-\Delta t)]$$\end{document}

Here, ui is saturated to [0,1], ∣κ∣ ≤ 2, ∣c∣ ≤ 1. We use the same pathways and control parameters (κ and c) for both legs. Our controller does not include cross-connections between the left and right legs. The normalized muscle fiber length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{L}$$\end{document} and force \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{F}$$\end{document} are defined as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{L}}_{i}(t)=\frac{{L}_{i}(t)}{{l}_{{\mbox{opt}},i}}-{l}_{0}{{\mbox{and}}}\,{\tilde{F}}_{i}(t)=\frac{{F}_{i}(t)}{{F}_{{\mbox{max}},i}}$$\end{document}

with lopt,i the optimal fiber length and Fmax,i the maximal isometric force of muscle i (see Table S1). We set l0 = 0.5 to ensure \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{L}}_{i}(t)$$\end{document} is within a realistic range while avoiding length inputs to become negative. All reflex gains and offsets remain constant for the duration of the simulation; there is no modulation or switching based on different states. Since the offsets represented by ci also remain constant, they do not constitute a feedforward control signal and are functionally identical to the length offsets represented by l0.

Optimization

Our reflex controller contains 71 free parameters (9 ci, 31 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\kappa }_{ij}^{L}$$\end{document} and 31 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\kappa }_{ij}^{F}$$\end{document} ), which dictate the behavior produced by the controller. We explore the capabilities of the reflex controller by defining five target gaits: forward walking (Walk Fwd) and hopping (Hop Fwd), backwards walking (Walk Bkw) and hopping (Hop Bkw), and running (Run).

For the five target gaits, based on the initial gaits we found, we extracted and manually tuned the initial state of the model (torso and joint angles and velocities) to roughly reflect the target gait (see Fig. 6 and Table S2). To produce a periodic motion, the reflex controller needs to find an equilibrium where the interplay of reflexes produces cyclically reoccurring states and sensor signals. At the beginning of the simulation, the system first has to initialize the feedback signals and enter the limit cycle. As such, the initial state can be regarded as a perturbation from which the controller needs to recover before entering the cyclic motion. Details on this initial behavior and, in general, the variability of the gait cycles are shown in detail in the Supplementary Material (Section 2). A well-chosen initial state minimizes the initial correction the controller must perform and helps to guide the optimization towards the intended target gait.Fig. 6. Initial states.For each target gait, the simulation starts in a different initial pose (joint angles and velocities). These poses are manually tuned to minimize the initial correction the controller has to make and to direct the optimization towards the target gait. The numeric values of all initial sets can be found in Table S2.

We do not optimize for specific motions through tracking objectives, but instead minimize a high-level cost function J, which consists of a velocity target and a term for minimizing muscle activation as a measure of muscle fatigue^70^:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J={{\mbox{w}}}_{{{{\rm{vel}}}}}\cdot {J}_{{{{\rm{vel}}}}}+{{\mbox{w}}}_{{{{\rm{act}}}}}\cdot {J}_{{{{\rm{act}}}}}$$\end{document}

and mainly dictates the direction and speed of the gait. Another included optimization criterion is continuous gait without falling for the duration of the simulation tsim. Simulations are terminated after tmax = 30 s or if the center of mass (CoM) of the model goes below 0.5 m. A gait fulfills this criterion if it does not fall within the maximum simulation time, i.e., if tsim = tmax This criterion is integrated into Jvel.

To optimize for the five different gaits with the same workflow, we furthermore introduce a flag s_v_ into Jvel, which allows to select between the different optimization conditions of reaching a target velocity vtgt (flag sv = 0) or maximizing the absolute velocity of the model (flag sv = 1 or  −1 for forwards and backwards movements, respectively):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J}_{{{{\rm{vel}}}}}=\left\{\begin{array}{ll}\frac{d}{{t}_{{{{\rm{max}}}}}}\hfill\quad &\,{{{\rm{if}}}} \, {{{\rm{sv}}}}\,=-1\hfill\\ | 1-\frac{\bar{v}}{{v}_{{{{\rm{tgt}}}}}}| +1-\frac{{t}_{{{{\rm{sim}}}}}}{{t}_{{{{\rm{max}}}}}}\quad &\,{{{\rm{if}}}} \, {{{\rm{sv}}}}\,=0\hfill\\ -\frac{d}{{t}_{{{{\rm{max}}}}}}\hfill\quad &\,{{{\rm{if}}}} \, {{{\rm{sv}}}=1}\hfill\end{array}\right.{{{\rm{with}}}}\,\bar{v}=\frac{d}{\,{\mbox{max}}({t}_{{{{\rm{sim}}}}},2)}$$\end{document}

Here, d is the completed distance, which corresponds to the horizontal displacement of the CoM of the segment closest to the origin at the end of the simulation (tsim). Jvel also encourages stability: if sv ≠ 0, Jvel is calculated based on the maximum simulation time tmax (set to 30 s), penalizing early falls where with the same velocity the reached distance d is smaller (and thus Jvel is smaller) the earlier the model falls. If the optimization target is to reach a target velocity (sv = 0), the average velocity of the model \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{v}$$\end{document} has to be calculated based on the actual simulation time tsim to be compared to the target velocity. To foster solutions that do not fall, a penalty for early terminations is added in this case ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1-\frac{{t}_{{{{\rm{sim}}}}}}{{t}_{{{{\rm{max}}}}}}$$\end{document} ). Furthermore, we only calculate the velocity based on tsim if tsim > 2 s, as otherwise the optimization gets caught in local minima where the model falls within the first 2 s at a velocity close to the target velocity.

The muscle activation cost term Jact consists of integrated cubed activation during the time interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[\frac{1}{3},\frac{2}{3}]{t}_{\max }$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J}_{{{{\rm{act}}}}}=\frac{1}{d}\left(\int_{t = \frac{1}{3}{t}_{\max }}^{\frac{2}{3}{t}_{\max }}{\sum }_{i}{a}_{i}{(t)}^{3}\right).$$\end{document}

Cubed muscle activation has been suggested as a measure of muscle fatigue^70,71^ based on experimental data indicating a cubic relationship between muscle force and endurance^72,73^. The restriction to the time interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=[\frac{1}{3},\frac{2}{3}]\,{t}_{\max }$$\end{document} helps to initially encourage finding gaits without falls before minimizing fatigue ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t < \frac{1}{3}\,{t}_{\max }$$\end{document} ) and prevents falling towards the end of the simulation due to micro-optimization of activations ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t > \frac{2}{3}\,{t}_{\max }$$\end{document} ). As gains are constant throughout the simulation and the gaits emerge from the interplay of the reflexes, an adaptation to minimize the activation during \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=[\frac{1}{3},\frac{2}{3}]\,{t}_{\max }$$\end{document} automatically influences the rest of the simulation as well. We do find that the gaits where activation minimization is active show some more inter-stride variability than the ones without; however, the strides are not different depending on which time interval they are in (compare Fig. S7a).

With the two cost function terms Jvel and Jact we then define the optimization parameters for the five target gaits. Humans typically walk at a preferred walking velocity, possibly minimizing energy expenditure^74^ and muscle fatigue^70^. The energetically optimal walking speed found by Ralston^74^ is at 1.23 m/s, and the subjects of the experimental dataset we use for comparison show a mean preferred speed of 1.17 m/s^38^. We therefore set a target velocity of 1.2 m/s for forwards walking. For simplicity, we set the same target velocity for backwards walking as well. Furthermore, we include activation minimization normalized to distance as a cost function term to take muscle fatigue into account^70^. For hopping and running, we optimize for maximum velocity and therefore deactivate the activation cost term. An overview of our optimization parameters can be found in Table 3. As discussed before, with this generic approach, the initial state has a large influence on the resulting gait, as e.g., hopping forward and running share the same optimization parameters and only differ in their initial state.Table 3. Optimization parametersParameterWalk FwdWalk BkwHop FwdHop BkwRunw_vel_100100100100100 –vtgt (m/s)1.2-1.2––– –s_v_001-11w_act_10001000000Chosen values of the optimization parameters for each of the five target gaits.

We optimize the parameters using the Covariance Matrix Adaptation Evolution Strategy (CMA-ES)^75^. The initial parameter mean values are rough estimates obtained through experimentation and are provided as part of the supplementary material. To allow sufficient exploration to cover all gaits and avoid local minima, we set the standard deviation to 0.2. We use the same mean values and standard deviations for all optimizations. Optimizations are terminated when the average relative improvement over the last 500 iterations is less than 10^−5^ per iteration.

Due to the complexity and discontinuity of our objective function (which involves a full musculoskeletal simulation in which the model can fall), our optimization problem is subject to many local minima. As a result, the optimization result highly depends on the stochastic sequence drawn by the CMA-ES optimizer. Due to these limitations, we run 20 optimizations for each target gait, each with different random seeds (which determine the stochastic sequence used by the optimizer), and pick the best result(s). While these results are unlikely to represent the absolute global optimum, we expect them to be close candidates.

Experimental data

We compare the kinematics and muscular activations of our model to experimental data for walking and running. While humans are capable of performing hopping forward and backwards, unfortunately, there is a lack of experimental data for comparison.

The experimental data of human walking is mostly extracted from Camargo et al.^38^. We use the treadmill trials of all participants and extract ten strides at the velocity matching the walking velocity of 1.2 m/s of our model. For each participant, we then average over these strides. Next, we compute the mean and the mean absolute deviation from the mean and display this interval. To showcase the diversity in human gait, we also compare to data of one participant, who is closest to our model in terms of weight and height (AB06, 1.8 m, 74.8 kg). As Camargo et al.^38^ does not contain electromyography (EMG) data for ILI and BF, we furthermore extract data for these two muscles from other sources. ILI is compared to data taken from Perry^76^ as used in Geyer and Herr^28^ and BF data comes from Blazkiewicz^77^. As Perry^76^ and Blazkiewicz^77^ do not provide data of multiple subjects, we directly use the extracted data and do not display an interval for ILI and BF. Note, that the EMG data from Camargo et al.^38^ was normalized for each subject with respect to the average amplitude at a speed of 1.35 m/s, while data for BF from Blazkiewicz^77^ is given in μV and ILI data from Perry ^76^ was normalized to the maximum manual muscle test value. Therefore, a comparison of magnitude to our simulation data, normalized to maximum activation, is difficult.

Backward walking data was digitized from van Deursen et al.^78^ who normalized EMG data for each subject to the percent of maximal voluntary contraction. For running, we digitized data from Hamner and Delp^79^ of experienced long-distance runners running at 3 m/s. EMG normalization in Hamner and Delp^79^ was performed for each muscle using the maximum voltage across all trials for each subject.

To enable a quantitative comparison between simulation data and experimental data, we follow the approach of Geyer and Herr^28^ and compute the maximum cross-correlation R of the mean signal. We also extract the corresponding time shifts Δ expressed in percent of the gait cycle (maximum shift \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Delta }_{\max }=20 \%$$\end{document} ). Again, interpretation of the fit to EMG should be done carefully, as it has been argued that cross-correlation might not be suitable for comparing EMG values (between different subjects)^80^.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Supplementary information

Supplementary material Description of Additional Supplementary Files Video S1 Target Gaits Video S2 Within-Gait Variability Video S3 Other Gaits Video S4 Video Abstract Reporting Summary