Steptostep variations in human running reveal how humans run without falling
Abstract
Humans can run without falling down, usually despite uneven terrain or occasional pushes. Even without such external perturbations, intrinsic sources like sensorimotor noise perturb the running motion incessantly, making each step variable. Here, using simple and generalizable models, we show that even such small steptostep variability contains considerable information about strategies used to run stably. Deviations in the center of mass motion predict the corrective strategies during the next stance, well in advance of foot touchdown. Horizontal motion is stabilized by total leg impulse modulations, whereas the vertical motion is stabilized by differentially modulating the impulse within stance. We implement these humanderived control strategies on a simple computational biped, showing that it runs stably for hundreds of steps despite incessant noiselike perturbations or larger discrete perturbations. This running controller derived from natural variability echoes behaviors observed in previous animal and robot studies.
https://doi.org/10.7554/eLife.38371.001eLife digest
Running at a constant speed seems like a series of repetitive, identical strides, but it is not. There are small variations in each stride. Selfinflicted errors in the forces generated by the muscles, or misperceptions from the senses, may cause these tiny imperfections. Uneven terrain or other outside forces, like a push, can also cause changes in a running stride. People must correct for these small changes as they run to avoid falling down. The only way to correct errors in a stride is by changing the force exerted on the ground by the leg.
Now, Seethapathi and Srinivasan document stepbystep how people correct for small imperfections in their running stride to avoid falling. In the experiments, eight people ran on a treadmill at three different speeds while the motion of their torso and each foot was measured along with the forces of each foot on the treadmill. Seethapathi and Srinivasan found that these runners corrected for minor deviations by changing where each foot lands and how much force each leg applies to the treadmill. The runners placed their foot at a different position on each step and these varying foot positions could be predicted by the errors in the body movement between steps. These errors in body movement could also be used to predict how the runners would change the forces applied by their legs on each step. Imperfections in the stride were mostly corrected within the next step. Errors that would cause the runner to fall to the side were corrected more quickly than errors in forward or backward motion. Seethapathi and Srinivasan incorporated these corrective strategies into a computer simulation of a runner, resulting in a simulated runner that did not fall even when pushed.
These findings may inform the design of robots that run more like humans. They may also help create robotic exoskeletons, prosthetic legs and other assistive devices that help people with disabilities move more fluidly and avoid falling.
https://doi.org/10.7554/eLife.38371.002Introduction
Human running is often modeled as being periodic (Blickhan and Full, 1993; Seyfarth et al., 2002; Srinivasan and Holmes, 2008). But running is not exactly periodic, even on a treadmill at constant speed. Body motion during running varies from step to step (Cavanagh et al., 1977; Belli et al., 1995; Jordan et al., 2007; Jordan and Newell, 2008). This steptostep variability could be due to internal perturbative sources like muscle force noise and sensory noise (Warren et al., 1986; Harris and Wolpert, 1998; Osborne et al., 2005) or small external perturbations (e.g. visual field inhomogeneity, small ground imperfections). To run without falling, the body’s ‘running controller’ must prevent the effects of these small perturbations from growing too large. Here, we provide an experimentally derived lowdimensional characterization of this control that reveals how humans run without falling down.
One classic modeling paradigm for running control assumes that the human leg behaves like a linear spring (Blickhan, 1989; McMahon and Cheng, 1990; Blickhan and Full, 1993). This paradigm has been used to argue how passiveelastic properties may reduce muscle work needed for locomotion (Alexander and Vernon, 1975; Alexander, 1990) and has been useful in examining locomotion in a simplified setting. Variants of these springmass running models have demonstrated stable running (Seyfarth et al., 2002; Seipel and Holmes, 2005; Ghigliazza et al., 2005; Geyer et al., 2006; Srinivasan and Holmes, 2008; Englsberger et al., 2016). These models have been successful in fitting the average center of mass motion during running (Blickhan and Full, 1993; Geyer et al., 2006; Srinivasan and Holmes, 2008). However, understanding running stability requires understanding how deviations from the average motion are controlled. It has been previously recognized that springlike leg mechanics cannot explain how deviations from the average motion are controlled and eventually attenuated (e.g. Ghigliazza et al., 2005; Biewener and Daley, 2007; Maus et al., 2015). Here, we examine the role of active muscle control in running stability, using more general models of human locomotion rooted in Newtonian mechanics (Srinivasan, 2011).
One way of characterizing the running controller is to apply perturbations (for instance, pushes or pulls or sudden changes in terrain) and examine how the body recovers from the perturbations (Van Woensel and Cavanagh, 1992; Daley and Biewener, 2006; Qiao and Jindrich, 2014; Riddick and Kuo, 2016). Instead of such external perturbations, here, we use the naturally occurring steptostep variability (Hurmuzlu and Basdogan, 1994; Maus et al., 2015) to characterize the controller. Previous attempts at examining such variability for controller information focused only on walking (Hurmuzlu and Basdogan, 1994; Wang and Srinivasan, 2012; Wang, 2013; Wang and Srinivasan, 2014) or considered variants of the springmass model (Maus et al., 2015). Here, we directly characterize the control in terms of how humans modulate their leg force magnitude and direction during running. The only way to control the center of mass motion is for the leg to systematically change the forces and the impulses it applies on the ground. We uncover how such center of mass control is achieved. We then implement this humanderived controller on a simple mathematical model of a biped (Srinivasan, 2011), showing that this biped model runs without falling down, despite incessant noiselike perturbations, large external perturbations, and on uneven terrain.
A humanderived controller such as the one proposed here could inform monitoring devices to quantify running stability or fall likelihood (O'Loughlin et al., 1993), or could help understand running movement disorders. Further, implementing such controllers into robotic prostheses and exoskeletons (Dollar and Herr, 2008; Shultz et al., 2015) will allow the human body to interact more ‘naturally’ with the device, rather than having to compensate for an unnatural controller. Some running robots have demonstrated stable running, using a variety of control schemes (Raibert, 1986; Chevallereau et al., 2005; Tajima et al., 2009; Nelson et al., 2019). But these robots fall short of human performance and versatility. Understanding human running may lead to better running robots.
Results
The steptostep variability during running appears superficially random. We show that this variability contains lowdimensional structure, specifically containing information about control strategies involved in running stably. We implement these strategies into a feedback controller, thereby stabilizing a simple mathematical model of a biped, and make further predictions.
We measured body motion and ground reaction forces (GRFs) of humans running for hundreds of steps on a treadmill at three speeds: 2.5, 2.7 and 2.9 m/s (see Materials and methods). Each running step consists of a ‘flight phase’ with neither foot on the ground, and a ‘stance phase’ with one foot on the ground. Figure 1 shows the coordinate system and sign convention: $x$ is sideways, $y$ is foreaft, and $z$ is vertical. The results we present are for the highest running speed and we discuss speeddependence of our results separately. Unless otherwise specified, all quantities and results in equations and figures are nondimensionalized using body mass $m$, acceleration due to gravity $g$, and leg length ${\mathrm{\ell}}_{\mathrm{max}}$. Forces are normalized by $mg$, distances by ${\mathrm{\ell}}_{\mathrm{max}}$, speeds by $\sqrt{g{\mathrm{\ell}}_{\mathrm{max}}}$, time by $\sqrt{{\mathrm{\ell}}_{\mathrm{max}}/g}$, impulses by $m\sqrt{g{\mathrm{\ell}}_{\mathrm{max}}}$, etc.
A hypothesized controller structure for stable running
During flight phase, the body center of mass moves in a nearly parabolic trajectory and the runner has no control over this parabolic motion (as the aerodynamic forces generated by the person are negligible, unlike birds). From Newton’s second law, it follows that the only way to control the center of mass motion is to modulate the total ground reaction force components during stance phase, when the leg is in contact with the ground. However, there are infinitely many ways to modulate the ground reaction forces to control the center of mass motion. Here, we examine how the ground reaction forces are modulated in response to center of mass state deviations during the previous flight apex (Figure 1). A flight apex is defined as when the center of mass height $z$ is maximum. The center of mass position and velocity at flight apex are denoted by $({x}_{a},{y}_{a},{z}_{a})$ and $({\dot{x}}_{a},{\dot{y}}_{a},{\dot{z}}_{a})$, respectively. Because the vertical velocity at flight apex ${\dot{z}}_{a}=0$ by definition, ${\dot{z}}_{a}$ is not considered as an explanatory variable. The absolute horizontal position $({x}_{a},{y}_{a})$ on the treadmill changes with a much slower timescale than other variables. Therefore, our default set of explanatory variables is $({\dot{x}}_{a},{\dot{y}}_{a},{z}_{a})$. We will include the horizontal position $({x}_{a},{y}_{a})$ when we comment later on ‘station keeping’.
Foreaft and sideways impulses independently control center of mass motion
The steptostep variability in the center of mass state at flight apex over hundreds of steps is shown in Figure 2a. To be stable, the runner needs to prevent this motion variability from growing without bound. As noted, the only way to control this motion is by using the ground reaction forces (GRFs). Consequently, the ground reaction force components over the stance phase are also variable (Figure 2b).
The net effect of the ground reaction forces on the center of mass velocity over a stance phase is captured by the force impulse, namely, the integral of the force. The variability in the sideways and foreaft ground reaction impulses over a step (Figure 2c) are wellpredicted by the variability in the center of mass state $({\dot{x}}_{a},{\dot{y}}_{a},{z}_{a})$ at the previous flight apex (Figure 3). Moreover, the sideways impulse depends primarily on the sideways velocity ${\dot{x}}_{a}$ and the foreaft impulse depends primarily on the foreaft velocity ${\dot{y}}_{a}$. Thus, it appears that the control in the foreaft and sideways directions are independent or decoupled. Pooled over all subjects, the bestfit linear model for the sideways impulse ${P}_{x}$ is:
and that for the foreaft GRF impulse ${P}_{y}$ is:
as in Figure 3. All coefficients in Equations (1) and (2) are significant at $p\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{10}^{4}$. Both sideways and foreaft impulses depend negligibly on vertical position deviations, so that including $z}_{a$ in the regression increases the $R}^{2$ values by less than $0.02$.
Almost deadbeat: impulses correct horizontal velocity mostly within a step
The linear models for the foreaft and sideways impulses in Equations (1) and (2) have a simple interpretation. The $\mathrm{\Delta}{\dot{x}}_{a}$ coefficient of about $1$ in Equation (1) (that is, $\mathrm{\Delta}{P}_{x}\approx \mathrm{\Delta}{\dot{x}}_{a}$) implies that sideways velocity deviations are completely corrected in one step, on average (over all steps and all subjects). This correction could have been done over many steps, as would be the case if the coefficient were $0.5$, say. But humans seem to exhibit a ‘onestep deadbeat controller’ on average for sideways velocity deviations (the term deadbeat control refers to when state deviations decay to zero in a finite amount of time). Of course, this singlestep correction is not perfect. An ${R}^{2}$ value of about $0.55$ suggests that the system overcorrects or undercorrects deviations for any given step.
Analogously, the coefficient of $0.72$ in Equation (2) suggests that about 72% of a forward velocity deviation is corrected in a single step, on average. While this is not strictly ‘deadbeat control’, it results in only ${(10.72)}^{2}=0.08$ of a perturbation remaining after two steps, and ${(10.72)}^{3}=0.02$ of a perturbation after three steps, indicating rapid control.
Apextoapex maps also show fast decay of center of mass deviations
We corroborate the above findings regarding perturbation decay with the ‘apextoapex maps’: that is, linear models that describe the relation between deviations in the state at one flight apex and those at the next flight apex. The righttoleft map from the state ${S}_{\mathrm{right}}={[{\dot{x}}_{a},{\dot{y}}_{a},{z}_{a}]}_{\mathrm{right}}$ at an apex preceding a right stance to the state at the next flight apex (preceding a left stance) is, approximately:
where the superscript $}^{\ast$ indicates that the coefficient is not significantly different from zero ($p\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0.05$). The lefttoright matrix ${K}_{\mathrm{L}\to \mathrm{R}}$ is similar to ${K}_{\mathrm{R}\to \mathrm{L}}$, except for the sign changes due to mirrorsymmetry. The matrix product of ${K}_{\mathrm{L}\to \mathrm{R}}$ and ${K}_{\mathrm{R}\to \mathrm{L}}$ — Jacobians of the Poincare map (Hurmuzlu and Basdogan, 1994; Guckenheimer and Holmes, 2013; Maus et al., 2015) — quantify how apex state deviations grow or decay over one stride (two steps). The eigenvalues of this matrix product were all less than one in absolute value, indicating a stable periodic motion. The largest eigenvalue was $0.14$, indicating that at most $14$% of a perturbation remains after a stride on average in any direction.
The low value of ${K}_{\mathrm{R}\to \mathrm{L}}(1,1)$, not significantly different from zero, suggests that a purely sideways velocity perturbation gets corrected essentially over one step on average, consistent with the sideways impulse control (Equation 1). Similarly, the value ${K}_{\mathrm{R}\to \mathrm{L}}(2,2)=0.27$ suggests that 73% of a forward velocity deviation is corrected in one step, consistent with the foreaft impulse control (Equation 2). Finally, the $(3,3)$ element of the step map (Equation 3) suggests that less than 50% of a deviation in vertical position ${z}_{a}$ remains after a step. See (Maus et al., 2015) for a detailed Floquet analysis of human running including more state variables, complementing the simplified version here.
Withinstep vertical impulse modulations control vertical position
The control of vertical position is qualitatively different from that of control in the foreaft and sideways directions, as we cannot use net vertical impulse for vertical position control due to the impulsemomentum considerations below. A flight apex occurs when the center of mass vertical velocity is zero. So, the net vertical impulse between two consecutive flight apexes is also zero (as it equals the change in vertical momentum, according to the impulsemomentum equation). Therefore, changing the net vertical impulse over a stance phase will not accomplish any meaningful control in the vertical direction. However, we will show that by differentially modulating the vertical impulse within one stance phase, we can change the vertical position (${z}_{a}$) from one flight apex to the next, without changing the net impulse.
To show this most simply, consider infinitesimal flight phases and a stance phase from $t=0$ to $t={T}_{\mathrm{step}}$. The total impulse ${P}_{z}$ due to the vertical ground reaction force ${F}_{z}(t)$ equals that due to gravity, which is given by, ${P}_{z}={\int}_{0}^{{T}_{\mathrm{step}}}{F}_{z}(t)\mathit{d}t={\int}_{0}^{{T}_{\mathrm{step}}}mg\mathit{d}t=mg{T}_{\mathrm{step}}$. For a triangular stance force (Figure 4) with peak force ${F}_{\mathrm{peak}}$ at ${t}_{\mathrm{peak}}$, we get ${F}_{\mathrm{peak}}=2mg$. Then, by integrating the vertical acceleration $({F}_{z}/mg)$ twice, the change in vertical position $z({T}_{\mathrm{step}})z(0)$ over a step is given by:
If the step was symmetric about midstance (${t}_{\mathrm{peak}}={T}_{\mathrm{step}}/2$), there is no vertical position change over a step ($z({T}_{\mathrm{step}})=z(0)$). The flight apex vertical position on the next step $z({T}_{\mathrm{step}})$ can be changed by changing ${t}_{\mathrm{peak}}$ relative to ${T}_{\mathrm{step}}/2$ (Figure 4). For example, if $z(0)$ at one flight phase was greater than its nominal value and the runner wishes to reduce it, this simple model predicts that the runner will decrease the firsthalf impulse and increase the secondhalf impulse; doing this is equivalent to delaying ${t}_{\mathrm{peak}}$ relative to ${T}_{\mathrm{step}}/2$ (as in Figure 4). This prediction is in agreement with the following experimentallyderived linear relations for the first half vertical impulse from $t=0$ to ${T}_{\mathrm{step}}/2$, namely ${\mathrm{\Delta}{P}_{z}}_{0}^{{T}_{\mathrm{step}}/2}$, and the second half vertical impulse from $t={T}_{\mathrm{step}}/2$ to ${T}_{\mathrm{step}}$, namely ${\mathrm{\Delta}{P}_{z}}_{{T}_{\mathrm{step}}/2}^{{T}_{\mathrm{step}}}$:
We see that a positive $\mathrm{\Delta}{z}_{a}$ corresponds to a decrease in the firsthalf vertical impulse and an increase in the second half vertical impulse. In addition to the vertical impulse, the landing leg length is also modulated in response to vertical flight apex deviations. Regressing the leg length $\mathrm{\ell}$ at the beginning of stance with the flight apex state, we found that this landing leg length is mostly a function of the vertical position at flight apex:
Thus, a downward position deviation at flight apex would result in landing with a shorter leg length than nominal (e.g. via flexed knee or ankle). A downward position deviation is analogous to a sudden stepup perturbation, so reducing the landing leg length reduces trip likelihood.
Impulse control is achieved by phasedependent force modulations
The linear models above tell us how deviations from nominal motion at flight apex are corrected grossly over the next stance. But they do not tell us how the forces are modified continuously throughout a stance phase. The variability of the GRF components $({F}_{x},{F}_{y},{F}_{z})$ depend on the ‘phase’ of the stride cycle, specifically, the time fraction ${\varphi}_{\mathrm{stance}}$ of stance (Figure 2b). To explain this phasedependent force variability within a single step, we compute the phase dependent sensitivity of $({F}_{x},{F}_{y},{F}_{z})$ to the center of mass state as follows. For each output, say ${F}_{x}$, we divide the stance duration into 20 phases and compute a linear model for ${F}_{x}$ at each of those phases, all with $({\dot{x}}_{a},{\dot{y}}_{a},{z}_{a})$ as inputs. We refer to the coefficients in these linear models as a function of the phase ${\varphi}_{\mathrm{stance}}$ as the phasedependent sensitivities of the GRFs (Figure 5) to the corresponding predictor variable in $({\dot{x}}_{a},{\dot{y}}_{a},{z}_{a})$.
The phasedependent sensitivity of sideways GRF to ${\dot{x}}_{a}$ shows that ${F}_{x}$ is decreased over the whole step to correct a positive sideways velocity deviation at flight and that a majority of this correction occurs during the middle of stance (Figure 5a). Similarly, in response to a positive foreaft velocity perturbation, the foreaft GRF is modulated so that there is a net negative force on the body over the next step (Figure 5b). The sensitivity of the foreaft force ${F}_{y}$ is more in the first half of stance than during the second half of stance, being modulated more during the deceleration phase (roughly ${\varphi}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e}}\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}0.5$) than during the acceleration phase (roughly ${\varphi}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e}}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0.5$).
Foot placement control: step in the direction of the fall
Placing the foot relative to the body allows a runner to modulate the leg force direction and thus the GRF components. The foot position $({x}_{f},{y}_{f})$ relative to center of mass position at the beginning of stance phase $({x}_{s},{y}_{s})$ is predicted by the previous flight apex state as described by the following equations. Specifically, sideways foot placement is described by the following equations:
The foreaft foot placement is described by the following equations:
That is, a sideways velocity perturbation to the body results in the foot being placed further along the direction of the perturbation. So, a rightward perturbation results in a more rightward step. Analogously, a forward velocity perturbation results in the foot being placed further forward relative to the body. As with the impulses, again, there is no significant coupling between sideways and foreaft directions. Foreaft foot placement modulation also depends on vertical position deviations, in a manner that the runner lands with a steeper leg when landing from a higher flight apex ${z}_{a}$. Such dependence of landing leg angle on vertical position is analogous to behavior in terrainchange experiments (Daley and Biewener, 2006; Müller et al., 2012; Qiao and Jindrich, 2012; BirnJeffery and Daley, 2012), as discussed in detail later. We speculate that using foot placement based on center of mass state may be an efficient way to affect the center of mass motion, compared to, say, changing the leg force magnitudes and leg joint torques after the foot touches down (Clark, 2018).
Swing foot repositioning happens during flight, just before foot touchdown
One possibility is that the foot placement deviations are achieved early on during the swing phase and this deviation is preserved during swing until the foot touchdown. However, this does not appear to be the case. Figure 6 shows the fraction of foot placement variance predicted by the swing foot state over the previous step. Less than 10% of the eventual foot placement is predicted by the swing foot at the beginning of flight phase (Figure 6). The explanatory power of the swing foot rises rapidly during the flight phase from less than 10% to a 100% when it becomes the next stance foot, suggesting that most swing foot repositioning may happen during this flight phase.
Center of mass state predicts future foot placement before the foot state does
At the beginning of flight phase (and earlier), the center of mass state is a vastly better predictor of the next foot placement than the swing foot itself (Figure 6). We can predict the foot placement using the center of mass state better than just the relative swing foot state until about 100 ms before foot touchdown. The explanatory power of the center of mass remains flat during flight. This flatness is likely because center of mass state follows a parabolic path during flight and thus accumulates no new variation. This lag between the explanatory power of the center of mass and the foot suggests that the error information in the center of mass state is yet to be fully incorporated into the swing foot repositioning until the flight phase. During the brief flight phase, when the swing foot’s explanatory power increases, information from center of mass state is transferred to the foot, presumably via some mixture of feedback control and feedforward dynamics.
Continuous stance state feedback, stationkeeping, and running speed do not significantly affect stance control
As an alternative to control based on discrete monitoring of deviations at the previous flight apex state, we considered a ‘continuous control’ model. Specifically, we obtained linear models for the GRFs based on the current center of mass state during stance $(\dot{x},\dot{y},z)$. These linear models did not differ significantly in the fraction of GRF variance explained, compared to the apexbased control model ($p=0.94$). In the linear models above, adding the sideways and foreaft apex body position $({x}_{a},{y}_{a})$ to the explanatory variables improves the ${R}^{2}$ values by less than $1.5\%$. Thus, the runners did not prioritize controlling their position relative to the treadmill (stationkeeping). Further, the regression coefficients for $({\dot{x}}_{a},{\dot{y}}_{a},{z}_{a})$ did not vary significantly across the three running speeds ($p\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0.05$, paired ttest).
Approximate leftright symmetry in the control
The running control gains have approximate bilateral (leftright) symmetry. The gains that couple sideways direction variables and either foreaft or vertical direction variables have mirrorsymmetry (see Equations 1,2,8,9). That is, these gains for the left leg’s stance are the negatives of corresponding gains for the right leg’s stance. On the other hand, gains that couple one sideways variable with another sideways variable, or one foreaft variable with another foreaft variable, are the same for the left and right legs without any such sign changes. This mirror symmetry in running control likely follows from the approximate mirror symmetry in body physiology about the sagittal plane and was also found in walking (Wang and Srinivasan, 2014; Ankaralı et al., 2015). This symmetry suggests the lack of a substantially dominant limb for running control, in contrast to the asymmetry and limb role differentiation that occurs in some other tasks (Peters, 1988).
Simple models of running capture GRFs and phasedependent GRF modulations
We now show that the experimentally derived control strategies described above are sufficient to control the running dynamics of a simple mathematical model of a biped. We consider a biped with pointmass upper body and massless telescoping legs capable of generating arbitrary force profiles (unlike a spring). We considered two versions of this biped model (Figure 7a), one with direct control of the leg force and another that produces leg forces via Hilltype muscles (Figure 7b). See Materials and methods for how the nominal running motion and the feedback controllers are specified for the models.
We find that the models’ ground reaction forces are similar to experimental data despite not explicitly matching the curves (Figure 8a). Further, we find that the phasedependent ground reaction force feedback gains for the models are qualitatively similar to the phasedependent gains inferred from experiment (Figure 8b), again, despite not explicitly fitting the shape of these phasedependent gains. This shows that these simple models can not only capture the average motion during running, but also how the runner responds to deviations from the average motion.
Humanderived controller stabilizes a minimal model of bipedal running
The simple models’ running motions are not stable without the controller: an arbitrarily small perturbation makes it diverge from the original running motion. With the foot placement and leg force controller turned on, the running motion is asymptotically stable. Figure 9a–d shows the model recovering from foreaft, sideways, and vertical perturbations at flight apex. It is a mathematical theorem that a stable periodic motion that can reject perturbations at one phase (say, flight apex) can reject perturbations at any phase (Guckenheimer and Holmes, 2013). So, it follows that our model rejects perturbations at any phase.
The inputs to the feedback controller (${\dot{x}}_{a},{\dot{y}}_{a},{z}_{a}$) do not include the absolute sideways and foreaft position $({x}_{a},{y}_{a})$ of the runner. Therefore, the controller does not correct position perturbations (stationkeeping). A sideways or foreaft push to the model results in convergence to the nominal running motion, except for a sideways or foreaft position offset (Figure 9c).
Nonzero leg work for energychanging perturbations
Figure 10 illustrates the leg workloop for the unperturbed run (net zero work) and when positive perturbations are applied to sideways and foreaft velocities, and vertical positions. All such positive perturbations result in net negative work on the first step after the perturbation, reflected in the workloops with net negative area within them. Such net positive or negative leg work is clearly necessary to recover from perturbations that change the total mechanical energy of the runner, as was recognized in prior discussions of the energyneutral springmass model of running (Ghigliazza et al., 2005; Biewener and Daley, 2007; Srinivasan and Holmes, 2008).
Explaining variability: muscledriven running model does not fall despite noise
To simulate the steptostep variability in real human running, we added ‘noise’ to our foot placement and leg forces (for the direct force control model) or muscle activations (for the muscledriven model) and simulated the biped models for a few hundred steps. This noise is meant to model the phenomenon that intended muscle forces tend to deviate from actual muscle forces due to motor noise (Harris and Wolpert, 1998). We find that while the direct leg force control model falls down, the runner with muscles does not fall down for hundreds of steps despite the noise. The stable motion of the center of mass in the presence of noiselike perturbations is shown in Figure 11a. The variability in the center of mass state at flight apex for the model (Figure 11b) as a result of the simulated noisy control is qualitatively similar to the variability found in experiment (Figure 2a). The model is also able to run without falling despite vertical position perturbations at flight apex, which are equivalent to uneven terrain. Thus, even though the model was derived using data on horizontal ground, it is capable of running robustly on uneven terrain. The muscledriven model is robust to motor noise presumably because of the intrinsic stabilizing properties of forcelength and forcevelocity relations (Hogan, 1984; Jindrich and Full, 2002).
Discussion
We have mined the steptostep variability in human running to show how humans modulate leg forces and foot placement to run stably. We then used these dataderived control strategies on a biped model, demonstrating robustness to discrete perturbations and persistent motor noise.
We have shown that humans use foot placement or leg angle control in a manner that they step in the direction of the perturbation, thereby directing the leg force so as to oppose the perturbation. This result provides an empirical basis for ad hoc assumptions about leg angle control made in previous running models (Seyfarth, 2003; Ghigliazza et al., 2005; Peuker et al., 2012). The foot placement controller derived from running data is qualitatively similar to the classic Raibertlike controller used in early running robots (Raibert, 1986) in that the foot placement opposes velocity deviations with no sidewaysforeaft coupling, but differs in that it has a dependence on vertical position perturbations. This makes such robotic controllers inadvertently biomimetic. Humans use similar foot placement control in walking, stepping in the direction of the perturbation (Hof et al., 2010; Wang and Srinivasan, 2014). Previous work had shown that appropriate foot placement is used in running ostriches while turning (Jindrich et al., 2007), running humans in cutting maneuvers (Besier et al., 2003), and turning while walking (Patla et al., 1999).
Some past work on inferring stability from variability focused on kinematic measures of variability such as Floquet multipliers, finitetime Lyapunov exponents (Dingwell et al., 2001) and long term correlations in walking and running variability (Hausdorff et al., 1995; Jordan et al., 2006; Kaipust et al., 2012). Such measures can provide discriminative diagnostic measures (Kaipust et al., 2012), but do not attempt to provide a causal narrative about how locomotion is controlled. Our approach here, rooted in NewtonEuler mechanics, is able to discover potential causal strategies underlying locomotion stability, and by extension, could inform treatment of pathological unstable movements in addition to diagnosis. Other past studies have used variants of the principal component analysis (Cappellini et al., 2006; Maus et al., 2015) to demonstrate that the intrinsic variability in human locomotion may reside in a lower dimensional manifold (Cappellini et al., 2006; Chang et al., 2009; Yen et al., 2009; Dingwell et al., 2010; Maus et al., 2015). Here, by focusing on how the center of mass is controlled through forces, we have implicitly used a physicsbased dimensionality reduction to examine the dominant control strategies.
While our work relies on linear regressions from data, the basic physics relating the inputs and outputs in these models suggest a natural causal account. This causal account, based on simplifying modeling assumptions, ignores the effect of variables not considered here. Our goal here was to delineate the explanatory power of controller descriptions based on center of mass state. To identify the effect of perturbations of other possibly relevant state variables (such as trunk attitude and angular velocity), we may need to either independently perturb these state variables or show that the natural variability in such variables is not significantly correlated with the center of mass state.
The gain relating sideways foot placement and sideways velocity deviation was about 2.5 times greater than the gain relating foreaft foot placement and foreaft velocity deviation; a similar factor of 3 was found in walking (Wang and Srinivasan, 2014), perhaps reflecting the greater sideways instability of a biped without foot placement control (Ghigliazza et al., 2005). Also consistent with lower control authority and a greater fall propensity in the sideways direction, we find that the recovery from a sideways perturbation is faster than from a foreaft perturbation. While stationkeeping was not prioritized over a single step, it may occur on a slower timescale with a multistep controller, not considered here.
The results we have presented have been for data pooled over all subjects. Performing the regressions for data from individual subjects indicates that the dominant terms in the inferred controllers are similar for all subjects; the subjecttosubject variability in the estimated control gains are shown in Figure 12. Figure 12 illustrates how the accuracy of an estimated control gain depends on the number of strides used for regression. For such linear regressions, the error estimate (standard deviation) is expected to decrease with ${N}_{\mathrm{stride}}$ like $1/{N}_{\mathrm{stride}}^{2}$, so that a factor of 10 decrease in error requires a 100fold increase in sample size (Wang and Srinivasan, 2012; Hamilton, 1994). This dependence on ${N}_{\mathrm{stride}}$ may guide selection of sample sizes for future experimental designs.
Our model predicts that when a runner starts at a higherthannormal height at flight apex, or equivalently, encounters a stepdown, the runner lands with a steeper leg angle (Figure 9d). Such behavior has been observed in humans and bipedal running birds running with large unforeseen or visible stepdowns (Daley and Biewener, 2006; Grimmer et al., 2008; Müller et al., 2012; Qiao and Jindrich, 2012). Conversely, stepups decrease touchdown angle, as predicted (BirnJeffery and Daley, 2012). This behavior has been attributed to swing leg retraction just before foot contact (Seyfarth, 2003), but our foot placement controller captures this phenomenon despite not having explicit leg swing dynamics. While the terrain perturbations in the aforementioned experiments were large (5–20 cm), our model is based on data with tiny steptostep deviations (vertical position ${z}_{a}$ s.d. 5 mm). This agreement indicates that humans may use qualitatively similar control strategies for large external perturbations and small intrinsic perturbations. Such foot placement control has also been used to control robots running on uneven terrain (Hodgins and Raibert, 1991).
It is expected that any running controller that achieves asymptotic stability will need to perform net mechanical work in response to perturbations that decrease or increase the body’s mechanical energy (Ghigliazza et al., 2005; Srinivasan and Holmes, 2008; Maus et al., 2015). Our results are consistent with such expectation, as illustrated by the workloops with net mechanical work in Figure 10. Energyconservative springlike leg behavior does not allow such net mechanical work and can achieve only partial asymptotic stability, not being able to handle energychanging perturbations (as noted by (Ghigliazza et al., 2005)). Indeed, it is generally thought that even the springmasslike steady state center of mass motion in running is due to considerable muscle action and has been termed pseudoelastic (Ruina et al., 2005) or pseudocompliant (McN. Alexander, 1997). Remarkably, energyoptimal running movements in models with no leg springs produce similar springmasslike center of mass trajectories (Srinivasan, 2011), with leg muscles performing equal amounts of positive and negative work.
A previous article (Maus et al., 2015) fit running data to variants of the springmass model, allowing the spring stiffness and spring length to change during stance, and showing that constant values for these parameters cannot fit running data. Here, we have used a simpler model to directly describe the control of stance leg force or activation (Figure 8). Such direct control of leg force or activation is perhaps more parsimonious than the simultaneous control of two variables, namely, spring stiffness and length. We have shown that humans modulate GRF continuously over the whole stance phase for control (Figure 5); Maus and colleagues (Maus et al., 2015) assumed, for simplicity, an instantaneous finite energy input at midstance.
The stabilizing responses we have characterized in this study are likely due to a mixture of feedforward dynamics and active neurally mediated feedback control. When we use the term "control" here, we implicitly refer to this mixture. It is hard to rigorously separate the roles of feedforward and feedback control without recording motor neuronal outputs and how these outputs interact with the properties of muscles. Nevertheless, we can determine the feasibility of feedback control by checking whether there is enough time for feedback control, given typical neuromuscular latencies. Our typical flight phase durations are greater than or about roughly equal to the typical short or middle latencies in reflex or feedback loops involving vestibular (Fitzpatrick et al., 1994; Iles et al., 2007) or proprioceptive mechanisms (Pearson and Collins, 1993; Sinkjær et al., 1999). This suggests feasibility of feedback based on flight phase or late stance phase information regarding center of mass state. While we have focused on the control of stance based on the previous flight apex, we have found that equivalent controllers based on the center of mass state at the end of previous stance have similar predictive ability (Figure 6), thus allowing more time for neural feedback. Specifically, the lag between the information in the center of mass state and the swing foot state regarding future foot placement is about 0.1 s for sideways placement and about twice that for foreaft foot placement, suggesting sufficient time for neurally mediated feedback control of foot placement (Figure 6).
Center of mass state or other body state information needed for feedback control could be estimated by the nervous system by integrating sensory signals from vision (Patla, 1997), proprioceptive sensors (especially when the foot is on the ground (Sainburg et al., 1995)), and vestibular sensors (Angelaki and Cullen, 2008), potentially in combination with predictive internal models (Wolpert et al., 1995; Cullen, 2004). In future work, repeating the calculations herein (for instance, Figure 6) for experiments that systematically block or degrade (say, by adding sensory noise) one or more of these sensors may tell us the relative contributions of these sensors to running control. We speculate that most available relevant sensory information is used, perhaps analogous to an optimal state estimator (Kuo, 2005; Srinivasan, 2009), and degrading one sensor may result in sensory reweighting on a slow timescale (Carver et al., 2006; Assländer and Peterka, 2014). Such experiments may also help explicitly distinguish the effects of sensory and motor noise, which we have implicitly combined here into a single residual term in the linear regressions.
In this work, we have obtained a running controller with simplifying assumptions. Because humans have extended feet, nonpointmass upper bodies and legs with masses, the simple pointmass model may not capture all aspects of the running data (Bullimore and Burn, 2006; Srinivasan and Holmes, 2008). Further, we have made simplifying assumptions about muscle architecture, muscle properties (linear forcevelocity relations), and muscle activation, which are meant to capture the main qualitative dynamical features of muscles, rather than model them quantitatively precisely. For instance, we used a linear forcevelocity relation, which may be sufficient to produce dampinglike and stabilizing muscle behavior when activated, but this damping behavior may be accentuated by a more realistic nonlinear forcevelocity relation.
Future work will also involve obtaining controllers for more complex biped models and muscle models with different control architectures, which, for instance, might include feedback control based on not just the center of mass state, but the states of individual body segments. We have focused on linear relations between state deviations and control, as this is naturally suited for small deviations and perturbations that our data contains. In future work, we hope to examine the range of perturbation sizes for which this linear description is accurate by comparing this linear control to responses in experiments with larger perturbations, also inferring nonlinear descriptions should they improve predictive capability. We will also examine other control architectures, for instance, more explicitly incorporating state estimation and considering continuous control of motor outputs based on an estimated state, partly correcting for neural latencies using internal models.
The methods used here are simple and noninvasive: they can be replicated to study running stability and control in other animals, or indeed, other approximately periodic tasks such as flapping flight and swimming. These methods are suitable for analyzing differences in different populations like athletes and nonathletes, the young and the elderly, and adults with and without movement disorders. Once such differences are wellcharacterized, this information could be used, say, in a rehabilitation setting to track progress from a controller in the presence of a movement disorder to a more healthy controller, and to design rehabilitation robots that assist in this progress.
Materials and methods
We collected running data by conducting human subject experiments, obtained linear models to characterize the control strategies hidden in the steptostep variability, and performed dynamic simulations using the inferred controller on a mathematical model of the runner.
Experimental methods
Request a detailed protocolThe protocols were approved by the Ohio State University Institution Review Board and subjects participated with informed consent. Eight subjects, three female and five male (age 25.0 ±5 years, weight 66.8 ±7 kg, height 1.8 ±0.14 m, leg length 1.05 ±0.08 m, mean ± s.d.) ran on a splitbelt treadmill at three constant speeds: 2.5, 2.7 and 2.9 m/s, presented in random order. Each speed had 2075 ± 67 strides across all subjects (one stride = two steps) with subjects running about 3.5 min on average. Subjects wore a loose safety harness that did not constrain their motion. Threedimensional ground reaction forces and moments on each belt of the treadmill were recorded by separate sixaxis load cells (Bertec Inc, 1000 Hz). Body segment motion was measured using markerbased motion capture (Vicon T20, 100 Hz) with four reflective markers on each foot and on the torso.
Calculating input state variables during flight apex
Request a detailed protocolWe define flight apex as when the center of mass velocity reaches its peak height ($\dot{z}=0$). The input to the running controller is drawn from the center of mass state at flight apex, namely position $({x}_{a},{y}_{a},{z}_{a})$ and velocity $({\dot{x}}_{a},{\dot{y}}_{a},{\dot{z}}_{a})$. Unless otherwise specified, we use the flight apex state $({\dot{x}}_{a},{\dot{y}}_{a},{z}_{a})$ as inputs in our linear models. The vertical velocity ${\dot{z}}_{a}$ at flight apex is zero by definition and hence not included as an input. The center of mass velocities are obtained by integrating the center of mass accelerations, that is, the massnormalized net force on the body: $\ddot{x}={F}_{x}/m$, $\ddot{y}={F}_{y}/m$, and $\ddot{z}={F}_{z}/mg$, where ${F}_{x},{F}_{y}$ and ${F}_{z}$ are the measured ground reaction forces on the body. To obtain the integration constants, we assume that the mean velocity and acceleration over the whole trial are zero, because the person does not translate appreciably in the lab frame over a trial. To remove the slow integration drift in the center of mass velocity, we used a highpass filter with a frequency cutoff equal to an eighth of mean step frequency (Luinge and Veltink, 2005; Schepers et al., 2009). Changing this highpass filter cutoff to a twentieth of the step frequency instead, or using a piecewiselinear detrending over 20 steps, do not affect any of this article’s conclusions. This is because the stabilitycritical timescales are much shorter. We ignored air drag here, because including it changed the velocities by less than 10^{5} ms^{1}, which is much smaller than the steptostep variability. We use a weighted mean of four markers, roughly at the sacral level, as an approximation of the center of mass position (Gard et al., 2004; Wang and Srinivasan, 2014; Perry and Srinivasan, 2017).
Calculating the output control variables during stance
Request a detailed protocolWe assume that the following variables are used to control the runner: GRFs, foot placement, and the landing leg length. Stance phases are identified as when the vertical GRF exceeds a threshold value to account for measurement noise (${F}_{\mathrm{z}}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}30$ N). The corresponding stance duration is ${T}_{\mathrm{stance}}$. The GRF impulses $({P}_{x},{P}_{y},{P}_{z})$ for each step are obtained by integrating the GRF components over the stance phase (${P}_{x}={\int}_{0}^{{T}_{\mathrm{stance}}}{F}_{x}\mathit{d}t$, etc). In addition to considering how GRF control occurs grossly over one step, we also consider GRF control throughout stance as a function of stance phase fraction ${\varphi}_{\mathrm{stance}}$. Each stance phase is divided into $n$ bins of duration ${T}_{\mathrm{stance}}/n$. To approximate how the GRFs changes with the stance phase fraction ${\varphi}_{\mathrm{stance}}$, we used the binned averages of the GRF in each of $n=20$ bins.
Linear regressions between the outputs and the inputs
Request a detailed protocolWe compute the mean values of the inputs over all steps in each trial and obtain deviations from these means $(\mathrm{\Delta}{\dot{x}}_{a},\mathrm{\Delta}{\dot{y}}_{a},\mathrm{\Delta}{z}_{a})$. Similarly, we compute the deviations from the means of the output variables $\mathrm{\Delta}F({\varphi}_{\mathrm{stance}})$, $\mathrm{\Delta}P$, and $\mathrm{\Delta}({x}_{f}{x}_{s},{y}_{f}{y}_{s})$. We use ordinary least squares regression to obtain linear models between the inputs and the outputs and report significant coefficients. Specifically, we have $\mathrm{\Delta}\mathrm{Output}=J\cdot \mathrm{\Delta}\mathrm{Input}$, where the Jacobian matrix $J$ represents the matrix of coefficients in the linear model. Each element of the matrix $J$ quantifies the sensitivity of an output variable to small changes in a corresponding input variable, as inferred from the data and subject to the simplifying model assumptions. These sensitivity coefficients could be interpreted as partial derivatives, such as: $\partial {T}_{\mathrm{stance}}/\partial {\dot{x}}_{\mathrm{a}}$, $\partial {F}_{y}({\varphi}_{\mathrm{stance}})/\partial {\dot{y}}_{\mathrm{a}}$, and so on. Unless otherwise specified, the results presented are based on deviations of all subjects pooled together as one dataset, but we find that the models of individual subjects’ data are qualitatively similar (as indicated in Figure 12). The coefficients for the right leg and left leg are computed separately, to accommodate sign changes due to symmetry about the sagittal plane.
Regressions with phasedependent inputs
Request a detailed protocolIn addition to the regressions described above using the flight apex state as the predictor, we used the center of mass state $(\mathrm{\Delta}{\dot{x}}_{a},\mathrm{\Delta}{\dot{y}}_{a},\mathrm{\Delta}{z}_{a})$ at different phases over the previous step to predict each of the stance phase outputs. Specifically, for each stance phase output, we performed $n=20$ regressions, each using the center of mass state at one of the $n=20$ equally spaced gait phases over the previous step, where one full step is defined as starting and ending at a touchdown. This analysis allows us to investigate the predictive ability of the center of mass state at different phases. For these phasedependent regressions, in addition to using the center of mass state as the predictor, we repeated the calculations using the swing foot state (position and velocity relative to the center of mass), so as to compare the different predictive abilities as in Figure 6.
Implementing the dataderived control on a minimal mathematical biped
Request a detailed protocolWe consider two simple models of running, similar in spirit to previous models in terms of simplicity (Blickhan and Full, 1993; Geyer et al., 2006; Srinivasan and Holmes, 2008), but generalized such that the leg forces are not constrained by ad hoc springlikeleg assumptions (Srinivasan, 2011). Instead, the biped controller details are inferred from our experimentally obtained linear models. Both biped models have a pointmass upper body and massless legs (Srinivasan and Ruina, 2006; Srinivasan, 2011), that can change effective leg length during stance by modulating the leg force (Figure 7a). During flight phase, the pointmass body undergoes parabolic free flight. The legs can apply forces on the upper body during stance phase. The two models, dubbed ‘direct force control model’ and ‘muscle control model’ differ in how the leg force is produced and controlled. For the muscle control model, we use a Hill muscle model with forcelength and forcevelocity relations (Figure 7b, c and d). The 3D equations of motion of the pointmass biped are: $m\ddot{x}={F}_{\mathrm{l}\mathrm{e}\mathrm{g}}\cdot (x{x}_{\mathrm{f}\mathrm{o}\mathrm{o}\mathrm{t}})/\ell$, $m\ddot{y}={F}_{\mathrm{l}\mathrm{e}\mathrm{g}}\cdot (y{y}_{\mathrm{f}\mathrm{o}\mathrm{o}\mathrm{t}})/\ell$, and $m\ddot{z}=mg+{F}_{\mathrm{l}\mathrm{e}\mathrm{g}}\cdot (z{z}_{\mathrm{f}\mathrm{o}\mathrm{o}\mathrm{t}})/\ell$, where ${F}_{\mathrm{leg}}$ is the scalar leg force, $({x}_{\mathrm{foot}},{y}_{\mathrm{foot}},{z}_{\mathrm{foot}})$ is the foot position with ${z}_{\mathrm{foot}}=0$ on flat terrain and $\ell =\sqrt{(x{x}_{\mathrm{f}\mathrm{o}\mathrm{o}\mathrm{t}}{)}^{2}+(y{y}_{\mathrm{f}\mathrm{o}\mathrm{o}\mathrm{t}}{)}^{2}+(z{z}_{\mathrm{f}\mathrm{o}\mathrm{o}\mathrm{t}}{)}^{2}}$ is the leg length from body to foot. In the ‘direct force control model’, the object of control is the leg force ${F}_{\mathrm{leg}}$ during stance phase. In the ‘muscle control model’, the object of control is the muscle activation ${a}_{\mathrm{muscle}}$, which is converted to muscle force via the forcelength and forcevelocity equations of Hilltype muscles (Figure 7a–b). See (Zajac, 1989; Srinivasan and Ruina, 2006; Srinivasan, 2011) for more detailed equations of motion and muscle model equations.
Both models have two terms in their control: (1) a feedforward or ‘nominal’ term, that depends only on the average or desired periodic motion and (2) feedback modification of the control in response to state deviations at flight phase. The model’s leg force or muscle activation is modeled as a twoterm sine series of the form ${A}_{1}\mathrm{sin}(2\pi t/2{T}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e}})+{A}_{2}\mathrm{sin}(2\pi t/{T}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e}})$, as shown in Figure 7e. By changing the relative weights of ${A}_{1}$ and ${A}_{2}$, the shape of the leg force profile can be changed from being symmetric about the peak force to being asymmetric, with the peak force preceding or following midstance. We parameterize the running motion using stance duration ${T}_{\mathrm{stance}}$, flight duration ${T}_{\mathrm{flight}}$, 2D foot placement $({x}_{\mathrm{foot}},{y}_{\mathrm{foot}})$, 3D initial conditions for stance $(x(0),y(0),z(0))$, and the coefficients of the twoterm sine series (${A}_{1}$ and ${A}_{2}$). We solve for these variables to obtain a periodic running motion that accurately match the forward speed, step period, step width, and peak leg force from experimental data (Figure 8a) by using an optimization procedure (Srinivasan and Holmes, 2008; Srinivasan, 2011) that enforces a constraint satisfaction tolerance of less than ${10}^{6}$. The runner leaves the ground when it reaches the maximum leg length, but the nominal leg length at landing is assumed to be shorter (95%) than the maximum leg length, as seen in running data (Voloshina and Ferris, 2015). We enforce that left and right stances are mirror symmetric. Unlike previous simple running models, our model’s nominal periodic motion has nonzero step width and a stance phase that is asymmetric about midstance. This asymmetric stance is due to unequal landing and takeoff leg lengths, and the asymmetry of the leg force or muscle activation about midstance.
The foot placement control for the models are based on the experimentally derived control and given by the linear model in Equations 8 and 9. The leg force feedback control based on apex body state, for the direct force control model, has gains as shown in Figure 8b. The muscle control model’s feedback control gains are also shown superimposed in Figure 8b. These control gains were derived for the two models by modifying the Fourier coefficients for the force and muscle activations respectively, so that the linear map from one apex to the next is best matched to that from data (Equation 3). While there are infinitely many controllers, even for this simple biped model, that can approximate the apextoapex map, our simplifying assumptions constrains the space of controllers to produce a unique fit. The leg forces and muscle activations are rectified, so that they never become negative despite feedback control (Blum et al., 2017). The foot placement control and leg force feedback control are activated only when the apex state deviates from nominal.
To obtain a running simulation over many steps, we break up each step into three phases: flight from apex to beginning of stance, the stance phase, and flight from the end of stance to flight apex. The control actions for the next stance are chosen at flight apex. As previously defined for the experimental data, the flight apex is when $\dot{z}$ becomes zero. In some cases, if the vertical velocity is downward when a stance phase ends ($\dot{z}\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}0$), there is no flight ‘apex’ and the controller uses the end of stance state instead of flight apex state as input. The end of flight and thus, the beginning of stance, are determined as the moment when the distance between the body and the target foot position is exactly equal to the landing leg length. The leg length at landing is also controlled based on flight apex state, based on the linear model in Equation 7. At flight apex, if the distance to the next foot position is less than the target landing leg length, the runner immediately goes into stance.
Such a simulation, when started from initial conditions exactly on the nominal periodic motion, results in a perfectly periodic motion when there are no further perturbations. We then resimulated the two pointmass running models for hundreds of steps, in the presence of noisy foot placements and leg forces or muscle activations with steptostep variability. To model the noise in foot placement, we computed the ‘desired’ foot placement based on the center of mass state at flight apex (Equations 89) and then added a deviation drawn from a normal distribution, whose variance equals the foot placement variance not explained by Equation 8. Similarly, we incorporated imprecise control of leg forces or muscle activation in the following manner: for each stance phase, once the leg force $F(t)$ (for model1) or muscle activation $a(t)$ (for model2) is determined based on the center of mass state at the previous flight apex, we ‘corrupt’ these functions by a multiplicative noise term, so that the actual leg force or muscle activation is $F(t)(1+\u03f5)$ or $a(t)(1+\u03f5)$ respectively, where $\u03f5$ is drawn from a normal distribution with variance equal to the unexplained steptostep variability in leg force magnitude. Thus, we use the unexplained variance in the foot placement and leg forces from experimental regressions as a simple model of the intrinsic noise in active control.
Data availability
All anonymized raw human running data (1GB) from this study are available at the following Dryad link: https://doi.org/10.5061/dryad.1nt24m0

DryadData from: Steptostep variations in human running reveal how humans run without falling.https://doi.org/10.5061/dryad.1nt24m0
References

Three uses for springs in legged locomotionThe International Journal of Robotics Research 9:53–61.https://doi.org/10.1177/027836499000900205

The mechanics of hopping by kangaroos (Macropodidae)Journal of Zoology 177:265–303.https://doi.org/10.1111/j.14697998.1975.tb05983.x

Vestibular system: the many facets of a multimodal senseAnnual Review of Neuroscience 31:125–150.https://doi.org/10.1146/annurev.neuro.31.060407.125555

Walking dynamics are symmetric (enough)Journal of the Royal Society Interface 12:20150209.https://doi.org/10.1098/rsif.2015.0209

Sensory reweighting dynamics in human postural controlJournal of Neurophysiology 111:1852–1864.https://doi.org/10.1152/jn.00669.2013

Mechanical step variability during treadmill runningEuropean Journal of Applied Physiology and Occupational Physiology 70:510–517.https://doi.org/10.1007/BF00634380

Muscle activation strategies at the knee during running and cutting maneuversMedicine & Science in Sports & Exercise 35:119–127.https://doi.org/10.1097/0000576820030100000019

Unsteady locomotion: integrating muscle function with whole body dynamics and neuromuscular controlJournal of Experimental Biology 210:2949–2960.https://doi.org/10.1242/jeb.005801

Birds achieve high robustness in uneven terrain through active control of landing conditionsThe Journal of Experimental Biology 215:2117–2127.https://doi.org/10.1242/jeb.065557

The springmass model for running and hoppingJournal of Biomechanics 22:1217–1227.https://doi.org/10.1016/00219290(89)902248

Similarity in multilegged locomotion: bouncing like a monopodeJournal of Comparative Physiology A 173:509–517.https://doi.org/10.1007/BF00197760

Force encoding in muscle spindles during stretch of passive musclePLOS Computational Biology 13:e1005767.https://doi.org/10.1371/journal.pcbi.1005767

Consequences of forward translation of the point of force application for the mechanics of runningJournal of Theoretical Biology 238:211–219.https://doi.org/10.1016/j.jtbi.2005.05.011

Motor patterns in human walking and runningJournal of Neurophysiology 95:3426–3437.https://doi.org/10.1152/jn.00081.2006

Modeling the dynamics of sensory reweightingBiological Cybernetics 95:123–134.https://doi.org/10.1007/s0042200600695

A biomechanical comparison of elite and good distance runnersAnnals of the New York Academy of Sciences 301:328–345.https://doi.org/10.1111/j.17496632.1977.tb38211.x

Whole limb kinematics are preferentially conserved over individual joint kinematics after peripheral nerve injuryJournal of Experimental Biology 212:3511–3521.https://doi.org/10.1242/jeb.033886

Asymptotically stable running for a FiveLink, FourActuator, planar bipedal robotThe International Journal of Robotics Research 24:431–464.https://doi.org/10.1177/0278364905054929

ThesisEnergetic efficiency and stability in bipedal locomotion: 3D walking and energyoptimal perturbation rejection. PhD thesis.The Ohio State University.

Sensory signals during active versus passive movementCurrent Opinion in Neurobiology 14:698–706.https://doi.org/10.1016/j.conb.2004.10.002

Local dynamic stability versus kinematic variability of continuous overground and treadmill walkingJournal of Biomechanical Engineering 123:27–32.https://doi.org/10.1115/1.1336798

Do humans optimally exploit redundancy to control step variability in walking?PLOS Computational Biology 6:e1000856.https://doi.org/10.1371/journal.pcbi.1000856

ConferenceDesign of a quasipassive knee exoskeleton to assist runningIEEE/RSJ International Conference on Intelligent Robots and Systems. pp. 747–754.https://doi.org/10.1109/IROS.2008.4651202

Biologically inspired deadbeat control for running: from human analysis to humanoid control and backIEEE Transactions on Robotics 32:854–867.https://doi.org/10.1109/TRO.2016.2581199

Compliant leg behaviour explains basic dynamics of walking and runningProceedings of the Royal Society B: Biological Sciences 273:2861–2867.https://doi.org/10.1098/rspb.2006.3637

Running on uneven ground: leg adjustment to vertical steps and selfstabilityJournal of Experimental Biology 211:2989–3000.https://doi.org/10.1242/jeb.014357

BookNonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector FieldsSpringer Science & Business Media.

Is walking a random walk? evidence for longrange correlations in stride interval of human gaitJournal of Applied Physiology 78:349–358.https://doi.org/10.1152/jappl.1995.78.1.349

Adjusting step length for rough terrain locomotionIEEE Transactions on Robotics and Automation 7:289–298.https://doi.org/10.1109/70.88138

Balance responses to lateral perturbations in human treadmill walkingJournal of Experimental Biology 213:2655–2664.https://doi.org/10.1242/jeb.042572

Adaptive control of mechanical impedance by coactivation of antagonist musclesIEEE Transactions on Automatic Control 29:681–690.https://doi.org/10.1109/TAC.1984.1103644

On the measurement of dynamic stability of human locomotionJournal of Biomechanical Engineering 116:30–36.https://doi.org/10.1115/1.2895701

Human standing and walking: comparison of the effects of stimulation of the vestibular systemExperimental Brain Research 178:151–166.https://doi.org/10.1007/s0022100607212

Mechanics of cutting maneuvers by ostriches (Struthio camelus)Journal of Experimental Biology 210:1378–1390.https://doi.org/10.1242/jeb.001545

Dynamic stabilization of rapid hexapedal locomotionThe Journal of Experimental Biology 205:2803–2823.

The structure of variability in human walking and running is SpeedDependentExercise and Sport Sciences Reviews 36:200–204.https://doi.org/10.1097/JES.0b013e3181877d71

An optimal state estimation model of sensory integration in human postural balanceJournal of Neural Engineering 2:S235–S249.https://doi.org/10.1088/17412560/2/3/S07

Measuring orientation of human body segments using miniature gyroscopes and accelerometersMedical & Biological Engineering & Computing 43:273–282.https://doi.org/10.1007/BF02345966

Constructing predictive models of human runningJournal of the Royal Society Interface 12:20140899.https://doi.org/10.1098/rsif.2014.0899

The mechanics of running: how does stiffness couple with speed?Journal of Biomechanics 23:65–78.https://doi.org/10.1016/00219290(90)900422

Invited editorial on “Interaction of leg stiffness and surface stiffness during human hopping”Journal of Applied Physiology 82:13–14.https://doi.org/10.1152/jappl.1997.82.1.13

Leg adjustments during running across visible and camouflaged incidental changes in ground levelJournal of Experimental Biology 215:3072–3079.https://doi.org/10.1242/jeb.072314

BookThe PETMAN and Atlas Robots at Boston DynamicsIn: Goswami PV A, editors. Humanoid Robotics: A Reference. Springer. pp. 169–186.

Incidence of and risk factors for falls and injurious falls among the Communitydwelling elderlyAmerican Journal of Epidemiology 137:342–354.https://doi.org/10.1093/oxfordjournals.aje.a116681

Online steering: coordination and control of body center of mass, head and body reorientationExperimental Brain Research 129:0629–0634.https://doi.org/10.1007/s002210050932

Legadjustment strategies for stable running in three dimensionsBioinspiration & Biomimetics 7:036002.https://doi.org/10.1088/17483182/7/3/036002

Compensations during unsteady locomotionIntegrative and Comparative Biology 54:1109–1121.https://doi.org/10.1093/icb/icu058

ConferenceShortcomings of thinking like a spring loaded inverted pendulum (SLIP) model: mechanical compensation to foreaft perturbations in human runningDynamic Walking Conference 2016.

Control of limb dynamics in normal subjects and patients without proprioceptionJournal of Neurophysiology 73:820–835.https://doi.org/10.1152/jn.1995.73.2.820

Ambulatory estimation of center of mass displacement during walkingIEEE Transactions on Biomedical Engineering 56:1189–1195.https://doi.org/10.1109/TBME.2008.2011059

Running in three dimensions: analysis of a Pointmass Sprungleg modelThe International Journal of Robotics Research 24:657–674.https://doi.org/10.1177/0278364905056194

A movement criterion for runningJournal of Biomechanics 35:649–655.https://doi.org/10.1016/S00219290(01)002457

Swingleg retraction: a simple control model for stable runningJournal of Experimental Biology 206:2547–2555.https://doi.org/10.1242/jeb.00463

Running with a powered knee and ankle prosthesisIEEE Transactions on Neural Systems and Rehabilitation Engineering : A Publication of the IEEE Engineering in Medicine and Biology Society 23:403–412.

Soleus longlatency stretch reflexes during walking in healthy and spastic humansClinical Neurophysiology 110:951–959.https://doi.org/10.1016/S13882457(99)000346

Optimal speeds for walking and running, and walking on a moving walkwayChaos: An Interdisciplinary Journal of Nonlinear Science 19:026112.https://doi.org/10.1063/1.3141428

Fifteen observations on the structure of energyminimizing gaits in many simple biped modelsJournal of the Royal Society Interface 8:74–98.https://doi.org/10.1098/rsif.2009.0544

How well can springmasslike telescoping leg models fit multipedal sagittalplane locomotion data?Journal of Theoretical Biology 255:1–7.https://doi.org/10.1016/j.jtbi.2008.06.034

ConferenceFast running experiments involving a humanoid robot2009 IEEE International Conference on Robotics and Automation. pp. 1571–1576.

A perturbation study of lower extremity motion during runningInternational Journal of Sport Biomechanics 8:30–47.https://doi.org/10.1123/ijsb.8.1.30

Biomechanics and energetics of running on uneven terrainJournal of Experimental Biology 218:711–719.https://doi.org/10.1242/jeb.106518

ThesisSystem identification around periodic orbits with application to steady state human walking. PhD thesisThe Ohio State University.

ConferenceSystem identification and stability analysis of steady human walking and the swing leg dynamicsASME 2012 5th Annual Dynamic Systems and Control Conference Joint with the JSME 2012 11th Motion and Vibration Conference. pp. 19–23.https://doi.org/10.1115/DSCC2012MOVIC20128663

Visual control of step length during running over irregular terrainJournal of Experimental Psychology: Human Perception and Performance 12:259–266.https://doi.org/10.1037/00961523.12.3.259

An internal model for sensorimotor integrationScience 269:1880–1882.https://doi.org/10.1126/science.7569931

Jointlevel kinetic redundancy is exploited to control limblevel forces during human hoppingExperimental Brain Research 196:439–451.https://doi.org/10.1007/s0022100918684

Muscle and tendon properties models scaling and application to biomechanics and motorCritical Reviews in Biomedical Engineering 17:359–411.
Article and author information
Author details
Funding
National Science Foundation (NSF CMMI grant 1254842)
 Manoj Srinivasan
Schlumberger Foundation
 Nidhi Seethapathi
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Andy Ruina for useful comments on an early draft. This work was supported by NSF CMMI grant 1254842 and a Schlumberger Foundation Faculty for the Future fellowship.
Ethics
Human subjects: The protocols were approved by the Ohio State University Institutional Review Board under protocol number 2012H0032. All subjects participated with informed consent.
Version history
 Received: May 15, 2018
 Accepted: February 6, 2019
 Version of Record published: March 19, 2019 (version 1)
Copyright
© 2019, Seethapathi and Srinivasan
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 4,767
 views

 428
 downloads

 27
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
 Medicine
Erectile dysfunction (ED) affects a significant proportion of men aged 40–70 and is caused by cavernous tissue dysfunction. Presently, the most common treatment for ED is phosphodiesterase 5 inhibitors; however, this is less effective in patients with severe vascular disease such as diabetic ED. Therefore, there is a need for development of new treatment, which requires a better understanding of the cavernous microenvironment and cellcell communications under diabetic condition. Pericytes are vital in penile erection; however, their dysfunction due to diabetes remains unclear. In this study, we performed singlecell RNA sequencing to understand the cellular landscape of cavernous tissues and cell typespecific transcriptional changes in diabetic ED. We found a decreased expression of genes associated with collagen or extracellular matrix organization and angiogenesis in diabetic fibroblasts, chondrocytes, myofibroblasts, valverelated lymphatic endothelial cells, and pericytes. Moreover, the newly identified pericytespecific marker, Limb BudHeart (Lbh), in mouse and human cavernous tissues, clearly distinguishing pericytes from smooth muscle cells. Cellcell interaction analysis revealed that pericytes are involved in angiogenesis, adhesion, and migration by communicating with other cell types in the corpus cavernosum; however, these interactions were highly reduced under diabetic conditions. Lbh expression is low in diabetic pericytes, and overexpression of LBH prevents erectile function by regulating neurovascular regeneration. Furthermore, the LBHinteracting proteins (Crystallin Alpha B and Vimentin) were identified in mouse cavernous pericytes through LCMS/MS analysis, indicating that their interactions were critical for maintaining pericyte function. Thus, our study reveals novel targets and insights into the pathogenesis of ED in patients with diabetes.

 Computational and Systems Biology
Largescale microbiome studies are progressively utilizing multiomics designs, which include the collection of microbiome samples together with host genomics and metabolomics data. Despite the increasing number of data sources, there remains a bottleneck in understanding the relationships between different data modalities due to the limited number of statistical and computational methods for analyzing such data. Furthermore, little is known about the portability of general methods to the metagenomic setting and few specialized techniques have been developed. In this review, we summarize and implement some of the commonly used methods. We apply these methods to real data sets where shotgun metagenomic sequencing and metabolomics data are available for microbiome multiomics data integration analysis. We compare results across methods, highlight strengths and limitations of each, and discuss areas where statistical and computational innovation is needed.