## Abstract

Running is an essential mode of human locomotion, during which ballistic aerial phases alternate with phases when a single foot contacts the ground. The spring-loaded inverted pendulum (SLIP) provides a starting point for modelling running, and generates ground reaction forces that resemble those of the centre of mass (CoM) of a human runner. Here, we show that while SLIP reproduces within-step kinematics of the CoM in three dimensions, it fails to reproduce stability and predict future motions. We construct SLIP control models using data-driven Floquet analysis, and show how these models may be used to obtain predictive models of human running with six additional states comprising the position and velocity of the swing-leg ankle. Our methods are general, and may be applied to any rhythmic physical system. We provide an approach for identifying an event-driven linear controller that approximates an observed stabilization strategy, and for producing a reduced-state model which closely recovers the observed dynamics.

## 1. Introduction

Running is an essential mode of human locomotion, during which ballistic aerial phases alternate with phases when a single foot contacts the ground. Motion capture can record the kinematics of running in great detail. Models that quantitatively fit these data have broad significance for prevention of falls, improving athletic performance, prosthetics, enhancement of human locomotion and the design of legged robots. The pogo stick, or spring-loaded inverted pendulum (SLIP) provides a starting point for modelling locomotion: McMahon & Cheng [1] and Blickhan [2] showed that the dynamics of an SLIP resemble those of the centre of mass (CoM) of a human runner. Here, we show that SLIP can reproduce experimental observations of CoM positions and velocities, removing at least the following fractions of the variance: 98.4% forward position, 95.5% vertical position, 75.1% lateral position, 46.9% forward velocity, 93.4% vertical velocity and 41.3% for lateral velocity. We show that at the same time as it so closely recovers kinematics, SLIP fails to reproduce the stability properties of humans. We then produce extensions to SLIP which maintain fidelity with the dynamics and also stabilize in ways that reproduce and predict human observations. Our methods are general and establish a paradigm for data-driven modelling of rhythmic phenomena. In particular, we provide an approach for identifying an event-driven linear controller that approximates an observed stabilization strategy. A similar linear controller could be used in robots or prosthetic devices that wish to mimic motion dynamics and control of humans or other study organisms.

We analysed data from individuals running on a treadmill in an attempt to predict their future motions based on the current state of their body. We quantified prediction quality using *relative remaining variance* (rrv)—the variance of prediction residuals divided by variance of the data. An ideal treadmill runner produces identical strides (pairs of consecutive steps) returning to exactly the same position after each. Human runners do not, and some of the variability of their motions represents corrections used to produce stability. In our analysis, we split the deviations from ideal running into parts that correlate with future steps and parts that are uncorrelated with future steps. We analysed the extent and nature of past information that is used. For example, do humans take account of more than one previous step in planning the next? Are there body parts whose movements predict future steps? The dynamics of the SLIP model itself give some insight into these questions.

Carver *et al*. [3] studied the ability of SLIP to execute prescribed trajectories: given a desired sequence of plausible CoM states at apex (the highest point of the ballistic flight trajectory), can parameters be chosen at each step that will execute this sequence? They discovered this is not possible, but that the model can produce all sequences specifying CoM state every stride (two steps). When we instead consider *reduced apex states*—height, forward and lateral velocity of the CoM at apex—we find that parameters can be found such that the model shows the same *reduced apex states* at every single step [3–5]. Inspired by these results, we hypothesize that humans running on treadmills correct deviations from their target trajectory within a single stride. Consequently, state at an apex should not depend on information more than one stride in the past.

Improving our understanding of SLIP has broader implications. Since its appearance in the 1980s, SLIP was shown to be a successful model for ground-reaction forces produced by legged animals of all sizes and leg numbers [6]. The success of SLIP prompted the formulation of the ‘templates and anchors hypotheses' (TAH) [7], which focus on dimensionality reduction as an underlying principle of locomotion dynamics: *animals restrict the dynamics of their many-degree-of-freedom bodies* [*the* ‘anchor’] *to a low-dimensional sub-manifold of their state-space upon which the behavior plays out* [*the* ‘template’]. The low-dimensional *template* models are hypothesized to capture the essence of locomotion across diverse species and scales, even though in each specific animal a more elaborate *anchor* model relates the template dynamics to that particular physical instantiation. The TAH engendered a large body of research in biomechanics, dynamics and robotics [8].

The mathematical perspective which underlies the TAH emphasizes asymptotic, long-term dynamics. From this viewpoint, steady locomotion is described by a *limit cycle oscillator*—a dynamical system whose long-term behaviour is a time periodic trajectory describing the motions an animal would perform under perfectly deterministic conditions. The motions we observe seem as those of a limit-cycle oscillator perturbed by the uncertainties of the body and its environment. The Floquet theory [9,10] of oscillators suggests that under such conditions there exist coordinates in which the dynamics near the limit cycle take a simple linear form, which can be estimated using data-driven Floquet analysis (DDFA) techniques [5,11–13]. DDFA plays a central role in our analysis of running.

A mechanical spring-mass hopper is not a complete SLIP model: there is no ‘first principles’ law for choosing model parameters such as leg length, spring stiffness and leg direction. In the context of running, the changes these parameters undergo from step to step are the crux of control, whereas in early SLIP models, the parameters were assumed fixed [6]. Here, we show a collection of models of human running that address the question of parameter selection for control, improving prediction of COM state after one step by a factor of at least 4.3, 5.2 and 23.1 for CoM height, forward and lateral velocity, respectively, compared with the classical, fixed-parameter SLIP model. We also show the fixed-parameter SLIP to be unstable and a poor predictor of human running over more than a single step in our experimental regime. In seeking an improved model for human running we applied DDFA, that is, we reconstructed the linearized dynamics around the estimated limit cycle [5,11–13], and show that this model predicts future motions substantially more accurately than SLIP—unsurprising in the light of the fact that it gives an upper bound on the predictive power achievable by any oscillator-like model with smooth dynamics. We then show a method for distilling predictive factors and state variables from the DDFA model, allowing them to be used to extend SLIP. One extended model is *factor-SLIP* with five extra states, which by construction reproduces the covariance of the prediction residual of the DDFA *full state*. We used *factor-SLIP* to guide the creation of a physically grounded *ankle-SLIP* which uses six extra states corresponding to the state of the swing-leg ankle and yields at least 44% (apex height), 75% (horizontal velocity) and 86% (vertical velocity) of the variance reduction achieved by the 94-dimensional *full state*. In summary, we show the first model of human running that [14] reproduces dynamics (forces, velocities, positions); [2] reproduces stability metrics (return map eigenvalues) and [6] predicts future motion up to the predicted maximum horizon of a full stride.

Our SLIP model is visualized in figure 1. The flight phase is ballistic, and the foot (end of the leg) remains stationary on the ground during stance. *Touchdown* transitions from flight to stance occur when the foot contacts the ground, and *lift-off* transitions from stance to flight occur when the foot loses contacts with the ground. Flight and stance are further broken down into two parts each by *apex* and *nadir* events which occur when the vertical velocity changes sign. Here, we define a step as the interval between two adjacent apices. At the apex transition of the flight phase, that is, at the beginning of a new step, the leg is pointed in a new direction in preparation for the next touchdown and its rest length changes. At nadir, the leg undergoes an instantaneous stiffness change with an accompanying change of rest length such that the leg force remains unchanged, allowing energy to be added or removed from the system [4]. In total, our SLIP model has four continuous domains and four discrete maps that are applied at the transitions between these. The changes in energy, leg stiffness and the location of the toe relative to the CoM in flight (in three dimensions) are five parameters that are reset at each apex (see figure 1 and Methods).

To complete the description of SLIP as a model of running, we must specify how the parameters are reset. There are many variants of these controllers in the literature [3,15–20]. A distinctive aspect of our analysis compared with previous studies is that we primarily used empirical data, rather than simulation, as input for the SLIP controller. For each step of our human runners, we determined SLIP parameters such that the model step matches observed minimum height, energy change, step duration, final apex height and direction [4]. Consequently, the simulated runner will reproduce the entire sequence of human running steps, as defined by the mentioned states, when the corresponding sequence of parameters for all steps and the observed initial apex state for the first step are provided. Given this set of apex states and the step parameters of the subsequent step, we searched for a control law which reliably predicts the observed sequences, thereby completing the SLIP model.

In §2, we describe the details of the SLIP models we investigated. We follow with a description of the experiments and the data reduction techniques, then present our main results followed by discussion of their significance. The final section presents technical aspects of our methods. They may be of interest to investigators who want to perform similar analyses on other locomotion systems or other physical systems which exhibit rhythmic behaviour.

## 2. Methods

### 2.1. Experimental protocol

Ten subjects participated in the study after giving informed written consent. Subjects were asked to run on an instrumented treadmill at moderate self-selected speed in a consistent manner. We excluded six of the subjects prior to further analysis, because of frequent marker loss (3), switching between forefoot and rearfoot running (1), extreme fatigue (1) and technical failures of the recording equipment (1). The remaining four subjects were male students with amateur running experience (weight: 64.2 ± 4.1 kg, age: 26.0 ± 1.8 years (mean ± s.d.)). After a few minutes of warm-up, subjects ran at self-selected speed on an instrumented treadmill (speed: 2.8 ± 0.2 ms^{–1} ≈ 10 ± 0.6 km h^{–1}). We recorded three-dimensional ground reaction forces at 1000 Hz using an ADAL3D-WR treadmill from TECMACHINE (HEF groupe), 42 166 Andrieux Boutheon, France, with four piezoelectric force sensors (type 9051A, Kistler, Winterthur, Switzerland) mounted below each of two metal plates located directly under the left and right half of the treadmill belt, and four three-dimensional piezoelectric force sensors (type 9077B, Kistler, Winterthur, Switzerland) mounted below the treadmill frame. We also recorded the motion of 31 reflective markers, placed at identical prominent anatomical landmarks (see the electronic supplementary material for details) on the subjects at 250 Hz using v. 10 Qualisys Oqus cameras (Qualisys, Gothenborg, Sweden). We acquired six recordings of 4 min duration from each subject, with a gap of approximately 30 s between them. In total, each subject provided about 1800 strides. We computed CoM motion using an established method based on a complementary filter [21]. We estimated apex times and states by fitting a quadratic model to the CoM height near apex. Using the *Phaser* algorithm [22], we computed an estimate of the dynamical phase of each data sample, allowing the motions to be parametrized by phase instead of time.

### 2.2. Data de-trending

To test stationarity, we computed a 61-stride centred moving average for kinematic data (marker position and velocities, sampled once per stride at apex) and SLIP parameters, for right and left steps. This moving average captures slow trends in the data and is referred to as ‘baseline’. To test if slow trends can be explained by a slow component of a (discretized) Floquet model, i.e. the kind of model we assume throughout our analysis, we compared the standard deviation of the empirical baselines with the baselines of 10 000 realizations of surrogate data.

Surrogate data were obtained by first computing a linear model of the (zero-mean) data at step *i* + 1 using the (zero-mean) data at the previous step *i* as predictor (see section below), and simulating this model with dynamical uncorrelated noise as input (i.e. a random process similar to the hypothesized form of the runner's dynamics). Next, each output dimension was rescaled such that it has the same standard deviation as the measured data after subtracting the mean. Without the presence of unmodelled slow-timescale dynamics, we would expect that 61-stride centred moving averages of these surrogate data frequently have larger standard deviation than the empirical baselines. However, we found that for the majority of kinematic states and SLIP parameters fewer than 10 of 10 000 surrogate datasets have a baseline with larger standard deviation than the original baseline, leading us to conclude that a statistically significant (*p* ≤ 10^{−3}) trend exists in the time series at this timescale.

Because we found significant trends in the majority of kinematic data and empirical SLIP parameters, we applied a 61 stride centred moving average de-trending for all data prior to any further analyses.

### 2.3. Empirical data

From the recorded data of the remaining four subjects, we computed all marker positions and velocities relative to the CoM, and excluded the absolute position of the runner on the treadmill as well as some redundant marker information from the state space. We resampled the data with respect to phase instead of time, using phase estimates obtained from Revzen & Guckenheimer [22], allowing strides to be compared with each other as the sets of 50 samples of 94-dimensional kinematic data. We removed slow trends in the data (see above). All subsequent analysis was performed on the detrended data; for further details of the data preparation and analysis, see primary sources [4,21,22].

Empirical data are available at the Dryad Digital Repository [23].

We used linear regression to obtain linear mappings from the state (typically kinematic state of the human runner, or parts thereof) at some phase to the state at later phases, or to SLIP parameters of the upcoming step, together with a bootstrap procedure to estimate the distribution of relative remaining variance (rrv) following application of the regression.

### 2.4. Bootstrap confidence intervals

The bootstrap procedure [24] is a parameter-free statistical method to obtain a distribution of a quantity of interest (e.g. a regression matrix or an rrv value, see §2.5) from a single dataset. First, a random subsample of the dataset is drawn. Then, the quantity of interest is computed on this subsample. These steps are repeated until enough realizations of the quantity of interest are obtained (typically 100–200 repetitions). If the quantity of interest is scalar (e.g. an rrv value), a confidence interval can be obtained from its distribution.

### 2.5. Linear regression and prediction

To estimate the relationship of kinematic states, kinematic states at later phases and SLIP parameters, we built corresponding linear models. Here, to obtain a linear model for a (possibly) vector-valued quantity *Y* (e.g. SLIP parameters, kinematic state at some phase *φ*_{2}) from a predictor *X* (e.g. kinematic state at phase *φ*_{1} < *φ*_{2}), we used ordinary least-squares linear regression on a set of corresponding observations of *Y* and *X*
where *k* denotes the *k*th iteration of the bootstrap procedure, and *k _{i}* the

*i*th sample of the

*k*th bootstrapped dataset (which becomes the

*i*th sample of the original dataset if the bootstrap procedure is omitted). The subscript ‘F’ denotes the Frobenius norm which is the square root of the sum of squares of all matrix elements.

In matrix notation, this reads
with the solution . Here, denotes the Moore–Penrose (right) pseudo-inverse of *X _{k}*, which we obtained by the scipy.linalg.pinv function from SciPy [25].

We assessed the quality of this model on the out-of-sample data , that is, those pairs of data that were not used for determining . We used the ratio of variance of to the variance of as the ‘relative remaining variance’ (rrv) measure reported in our results. We performed several (typically 100–200) bootstrap iterations on the *X*,*Y* pairs to obtain confidence intervals for the rrv numbers reported.

Contrary to naive intuition, an rrv value 1 still indicates some prediction ability, because finite sample effects affect the prediction matrix and allow some of the predictor noise to add to the pre-existing uncertainty in the data. rrv values larger than 1 are possible.

For each subject, we performed independent regressions and predictions.

### 2.6. Modified spring-loaded inverted pendulum model

The SLIP model we used is a hybrid dynamical system representing the movements of a point mass in three dimensions. Its equations of motion are given in figure 1 and were integrated numerically using a Dormand–Prince ODE integrator [26]. The model is described by a mass *m* at position *r* on top of a massless telescopic linear spring leg with rest length *l*_{0} and stiffness *k*. During flight, the leg's orientation is described by the inclination angle *α* and the azimuthal angle *β*, with *α* = *β* = 0 the forward direction on the treadmill (the *x*-direction). An auxiliary parameter Δ*E* controls the net work done by the leg while in contact with the ground. This energy change is introduced at nadir () by updating the leg rest length *l*_{0} and stiffness *k* according to
2.1which leaves the leg force unchanged.

Each subject's mass was estimated from the measured ground reaction forces separately for each trial to account for weight changes owing to sweating and drinking (fluctuations typically around 0.4 kg ≈ 0.6% body mass). We normalized the leg stiffness and energy change to the mass for each trial, because the mass appears only as a constant factor in those terms. That way, varying subject mass does not directly affect SLIP parameters.

The stance dynamics end when the leg length *l* reaches , and transition to flight.

At apex, new values for the parameter vector are chosen using a specified linear control law. Note, however, that the changes in (*k*, *l*_{0}) are applied at the subsequent nadir rather than at apex.

### 2.7. Empirical spring-loaded inverted pendulum parameter estimation

In our SLIP model, a step is the interval between two adjacent apices. Using the method of Ludwig *et al*. [4] on each step of a human runner, we can find SLIP parameters such that in the simulation model they result in an exact match of the final apex velocity and height, step duration and nadir height. It is important to note that the new SLIP parameters, together with initial conditions, contain all information about the final CoM state at apex.

### 2.8. Controlled spring-loaded inverted pendulum models

Controlled SLIP models are hybrid dynamical system models which propagate the CoM state in continuous time using a forward simulation of SLIP. Non-CoM states are propagated from apex to apex using discrete linear maps. Initial conditions are taken from experimental data or from the last simulated step, respectively. For each step *i*, SLIP parameter *p ^{i}* updates are performed at initial apex using an affine control law
2.2The reference parameters (set points of the affine control law) for right and left steps (subscripts R and L) correspond to a periodic orbit with the same apex states, step durations and nadir heights as the average experimental steps for the selected subject. The parameter prediction matrix (the ‘feedback law’) was obtained via a least-squares regression of the corresponding model state

*x*(CoM states together with non-CoM states off their respective baselines) against the predictor variables for that model.

### 2.9. Return map eigenvalues

To obtain return maps for kinematic data, we performed linear regression on the difference between the kinematics and the baseline kinematics, regressing them against the same data two steps (one stride) ahead.

To test plausibility of the results, we additionally computed return maps in a different way, using multiple Poincaré sections. Multiple sections take more empirical data into account, but in theory, the recovered eigenvalues should not be affected by the number of Poincaré sections taken into account. In practice, other groups have reported advantages of this approach [27]. We selected *N* () Poincaré sections with equal phase difference, and computed bootstrapped sets of maps *B _{i}*

_{,k}mapping from one Poincaré section

*i*to the next Poincaré section

*i*+ 1. We computed the return maps

*C*as a product of the maps : . For each set of return maps—one for every number of Poincaré sections used—we computed the corresponding set of eigenvalues and plotted the distribution of those eigenvalues for each set in figure 9.

_{k}Eigenvalues for reduced-state space definitions like *ankle-SLIP* were computed in an analogous way, replacing the full kinematic state with the corresponding reduced-state space.

For controlled SLIP models, the corresponding eigenvalues are the eigenvalues of the Jacobian around the periodic solution. To estimate these Jacobians, we used two-sided finite differences of ±0.01 for each model dimension. We verified that the recovered eigenvalues did not depend on the size of the differences, within the precision we report.

## 3. Results

The primary goal of the work we present here is to predict CoM motion in human runners. One of our contributions is providing a standard against which running models can be tested by applying the recently developed techniques of DDFA [9,10,12,13,28] to the dynamics of the CoM, which are presumed to evolve in the stability basin of a stable limit cycle. We assume that the small magnitude of deviations from this hypothetical limit cycle justifies the linear approximations we used here. We constructed the linearized dynamics of states near the limit cycle to form a state-evolution model (the *full state* model in the sequel), which includes estimates of the linearized return map (stride map; Poincaré map [10]) commonly used for stability analysis.

The predictive performance of this *full state* model establishes a reference ‘gold standard’ for evaluating smaller template models of running. In the *full state* model, we applied ordinary least-squares regression to the complete 94-dimensional state of motion capture measurements and CoM, using those as predictors for that state at later phases. Similarly, we computed a linear map using *full state*'s state to predict SLIP parameters for the upcoming step. Because *full state* has access to the richest collection of predictors of the models we examine, it is expected to have the best predictive ability.

Our analysis compared two classes of models—simulated models and linear models. In the simulated models, we took CoM motions during a step as evolving according to numerical integration of the SLIP equations from the initial apex conditions. In the linear models, each future state was computed as the image of a phase-dependent linear map of the initial apex state.

In the *full state* SLIP, we assumed that the update law for SLIP parameters (the *controller*) is given by a linear map from the complete set of 93-dimensional kinematic data (excluding vertical velocity which vanishes at apex) to SLIP parameters. We determined the coefficients of this discrete linear controller by regression, minimizing the residual between the predicted and observed SLIP parameters of the upcoming step. This controller serves as a benchmark upper bound for prediction performance for other controllers examined in this paper, as it is free to act with a different linear map at every phase while having access to the most state information data.

The *full state* SLIP is not an autonomous system: future model states are not given purely as a function of past states; instead, they depend on the totality of measurements we collected. To create a self-contained model of running, we first examined the simplest alternative—the *mean parameter-SLIP* model—which is closest to the traditional SLIP models of the 1990s [8,29,30]. In the *mean parameter-SLIP* model, we held the SLIP parameters of left and right steps constant, with values that correspond to a periodic motion with the same step durations, nadir heights and apex states as the average empirical left and right steps. We found the parameter values by solving for a periodic trajectory were within 0.1 standard deviations of the mean parameter values for left and right steps, respectively. We obtained predictions for individual running steps from *mean parameter-SLIP* by starting the simulation with initial conditions (here: CoM state) taken from individual running steps in the empirical data. We found that *mean parameter-SLIP* does many times worse than *full state* when it comes to prediction of human responses. It *increases* the residual of apex height, forward velocity and lateral velocity by at least a factor of 4.3, 5.2 and 23.1, respectively, across subjects.

After establishing that a large gap in predictive ability lies between *mean parameter-SLIP* and *full state* SLIP, we set out to construct autonomous low-dimensional models whose predictive performance comes close to that of the *full state* model. One strategy uses linear maps similar to the *full state* model, but first reduces the full state by projection onto a subspace. The model's state on that subspace is then predicted for the subsequent apex using a (discrete time) linear map. Each low-dimensional model thus contains a linear map that predicts SLIP parameters, and another linear map that predicts the non-CoM states of the model.

The *full state* map from state to parameters is rank 5, as there are only five parameters. Inspired by factor analysis computations in statistical inference, we chose the subspace orthogonal to the null space of this rank 5 map (the *factors*) together with reduced CoM state at apex, giving the *factor-SLIP* model's state. In this model, each set of SLIP parameters close to the periodic orbit corresponds to a unique linear combination of five computationally obtained factors. In addition to predicting the values of SLIP parameters for the next step, we used *factor-SLIP*'s state to predict the factor values at the next apex, thereby allowing for a self-contained model to be produced if the data support it. By construction, this reduced-state discrete time linear controller produces the same output covariance as the original *full state* control law for SLIP parameters for the first step (see the electronic supplementary material for details). It does not necessarily produce the same output covariance for its own state at next apex, leading potentially to degraded prediction performance in subsequent steps. We test the ability of *factor-SLIP* to predict its own state at next apex in figure 2, and find that *factor-SLIP* can comparatively well predict its own future state compared with the reference given by the *full state* model, qualifying *factor-SLIP* as a self-contained reduced model.

Figure 3 illustrates the 5 × 93 projection matrix of the *factor-SLIP* model. Note that the components of this matrix that come from ankle position and velocity make the largest contributions to the prediction of SLIP parameters. This motivated our second choice of model: *ankle-SLIP*, whose state consists of the CoM state and the ankle position and velocity of the next stance leg. The *ankle-SLIP* uses the mechanical state (position and velocity) of the CoM, and the mechanical state of an additional mass—the protracting swing leg—to generate a prediction for both SLIP parameters and ankle states at the next apex. Unlike *factor-SLIP*, the additional discrete time states added in *ankle-SLIP* have a direct physical meaning. We performed the same analysis of autonomy for *ankle-SLIP* as we did for *factor-SLIP* (figure 2), and found that this subsystem is not fully autonomous, but predicts most of the change in its own future states.

Our final simplified controller explores the use of *augmentation*, whereby parameters of a dynamical system are turned into state variables. We used the SLIP parameters themselves as discrete-time state variables and computed linear control laws for left and right steps, in each of which the state of the SLIP at apex is augmented with the values of all SLIP parameters from the previous step. The resulting *augmented-SLIP* model takes a linear map of the CoM apex state and the previous step's SLIP parameters into a new set of values for the SLIP parameters. The *augmented-SLIP* model gives an upper bound on what a linearized controller with access only to SLIP state information could accomplish.

We compared the ability of each of these models to predict the SLIP parameters humans used for the next step (figure 4). *Factor-SLIP* performs as well at prediction as the reference *full state. Ankle-SLIP* shows slightly poorer, but comparable performance, whereas *augmented-SLIP* has much poorer performance. As these SLIP parameters together with CoM state constitute the state of *augmented-SLIP*, the poor performance implies that no affine feedback on the *augmented-SLIP* apex state reproduces the SLIP parameter choices that humans use. By comparison, *full state* and even *ankle-SLIP* perform much better prediction. We conclude that linear prediction is in fact possible, but that the necessary information is unavailable in the *augmented-SLIP* state. Furthermore, these results imply that accurate prediction of human CoM state requires non-CoM state information.

The analysis of Carver *et al*. [3] suggests that runners should return to their periodic motions within a stride following disturbances. Consequently, we hypothesized that deviations from the limit cycle remain predictable (i.e. retain information) only during a single stride, and tested for the expected loss of prediction ability as a function of phase. In figure 5, we present the ability of phase-dependent linear models with the various predictors to predict state variables over a horizon of one stride into the future. Our results demonstrate that nearly 80% of the variance was not predicted one stride into the future in any of our choices of state variable and prediction algorithm, whereas for some combinations (e.g. CoM lateral velocity as predicted by the *full state* model) 80% of the variance remained predictable when looking only one step into the future. Overall, the *full state* provided a benchmark against which the reduced models were evaluated. Note that *augmented-SLIP* does poorly at prediction across the board, whereas *factor-SLIP* does nearly as well as *full state*. The more physically grounded *ankle-SLIP* did only a bit worse than *factor-SLIP*, except for its greatly reduced ability to predict CoM height beyond the next lift-off. Possibly, this was owing to left–right asymmetry. *Ankle-SLIP* has a single set of state variables used for both ankles, with independently determined controllers for left and right steps. It is possible that an *ankle-SLIP*-like controller with separate states for left and right ankles could reduce the performance gap.

As anticipated by the theoretical predictions of Carver *et al*. [3], adding information from previous steps did not lead to improved prediction performance (figure 6). We found no evidence in our data that human gait control uses any form of multi-stride memory, neural or mechanical.

The models presented to this point have been phase-dependent linear functions, derived by regression to optimize prediction of SLIP parameters. An alternative approach we explored consists of using the corresponding SLIP models in forward simulations for prediction. The models differed only in the control law which acted at the apex event. In figure 7, we compared the predictions for linear models with their simulation counterparts.

SLIP does not exactly reproduce the average CoM motion during a step (figure 7*d–f*), but can reproduce each desired reduced apex state. This results in a systematic prediction error during the step, that vanishes at apex. Figure 7*b,c* shows that at phases where this systematic error is small, namely in the flight phase (around the middle of each stride) for horizontal and lateral velocity, we see good agreement between the simulated CoM states and the corresponding predictions from linear models. In fact, the simulated and linear models have almost the same predictive performance for a single step at apex, and still perform similarly after a stride (figure 8). This shows that the choice of the state space is more important for prediction performance at apex than the choice of the model type (simulation or linear).

As an additional check for dynamic similarity between different models, we compared their return map eigenvalues in the complex plane. The eigenvalues express how quickly nearby trajectories approach the limit cycle. They are invariant under smooth coordinate changes, so differences in the eigenvalues give a measure of the dynamic dissimilarity of two systems. We computed eigenvalues estimated with different numbers of Poincaré sections (*N* = 1, 2, 3, 5), and computed the mapping between subsequent sections via least-squares regression for the experimental kinematic data, taking the return map as the product of these maps (see Methods for details). The observed eigenvalues depended on the number of sections used (figure 9). It has previously been noted that empirically computed eigenvalues depend on the sections at which they are computed, an apparent contradiction to their theoretical properties [11]. A more complete understanding of the limitations of eigenvalue estimation with one or more sections remains a topic for further investigation. We then computed the return maps for our different simulated SLIP models and the corresponding linear models which share the same state space. Our results show a correspondence in eigenvalues between the simulated SLIP models and their linear counterparts, providing evidence for the dynamic similarity of SLIP and linear models. This reinforces our finding that the choice of state space is more important for prediction of apex states than the choice of the model type (SLIP or linear). Furthermore, the magnitude of the reduced models' eigenvalues is in good agreement with those of the kinematic data which correspond to the *full state* SLIP.

## 4. Discussion

Our analysis leads us to propose that CoM motions of human running can be modelled by *ankle-SLIP*, a hybrid dynamical model with states comprised of the CoM and the position and velocity of the swing-leg ankle. The future values of the five SLIP parameters, and six swing-leg ankle states (three each of positions and velocities) are approximately predicted by the current apex state and swing-leg ankle state. The reduction to *ankle-SLIP* might in part be attributable to the small number of subjects, or to the constraints of the experiment—subjects had to maintain the same average velocity and direction to stay on the treadmill. A larger subject population or less constrained running might yield a different choice of reduced model. Nevertheless, our methods could still be applied to produce appropriately reduced models which achieve close correspondence with the observations.

Note that our model only maintains the state for a single ankle even though that ankle alternates between the left and right leg. One interpretation is that during stance any dynamical state of the stance ankle is lost, being entirely encoded in the CoM state, thereby admitting a description that has merely one ankle instead of two. This form of dimensionality reduction is unique to hybrid oscillators, and was only described recently [31].

Runners are not, generally speaking, bilaterally symmetric. In our results, this is apparent in that the right apex map of subjects is not only different from the mirrored left apex map, abut it also predicts significantly different outcomes (SLIP parameter, CoM and ankle states); details of our results regarding gait control asymmetry are relegated to future work.

The fact that we find ankle states to be an important predictor of CoM control does not imply the existence of an active controller in the nervous system that uses ankle states as input. Ankle states might merely be an output of some control process that we do not observe, and this process need not be neural in nature. For example, leg retraction right before touchdown has been shown to stabilize SLIP [19] and does not necessarily require feedback from the nervous system. The dynamics of the body may modulate the swing ankle mechanically, thereby modifying the touchdown event and influencing stability. While swing leg retraction shows the importance of fore–aft motions of the ankle, our analysis has also shown ankle position relative to CoM to be an important predictor. Hence, *ankle-SLIP* contains the full three-dimensional ankle state.

The choice of using the system's state at apex for prediction is motivated by literature, not by physiology. The SLIP model is insensitive to the orientation and length of the leg during flight: these become relevant only when the foot contacts the ground. We do not assert that the use of apex data to compare human runners with our model is optimal or superior to comparisons made at other phases of the gait cycle, but SLIP controllers have used apex states as feedback.

Under our experimental conditions with comparatively slow running speed, a constant parameter SLIP does not reproduce the qualitative behaviour of the observed perturbation response. We could further show that SLIP needs an extension of the state space to be able to display human-like reactions. However, at higher velocities, a constant parameter SLIP can produce stable periodic gaits [18] which humans might rely on at faster running speeds. The methods presented here allow one to readily test such hypotheses, either directly by comparing eigenvalues and eigenvectors to model predictions, or through perturbation analysis. Even at high speeds, the constant parameter SLIP has Floquet multipliers far from zero—it is emphatically not a deadbeat controller—and thus perturbation recoveries are expected to persist over multiple strides. In general, measuring transient responses to stimuli that perturb the periodicity of running motions is a powerful technique for probing the structure of running dynamics and control [32–35].

Our results show that SLIP, as a parametric model of human CoM trajectories during a step, can be used in conjunction with a simple discrete time linear control law to obtain a prediction of human running motions that is nearly as good as allowed by theoretical predictions [3]. The study of control exercised at isolated instants, triggered by time or sensory events, is a quickly burgeoning field in control theory [14,36].

The Poincaré maps of the periodic solutions of SLIP in the region relevant to our data are smooth, and thus linearizable at the fixed point. Equipped with an arbitrary linear controller, it is not surprising that SLIP can be made to match the eigenvalues of observed human dynamics. Our contribution is not in asserting that such a control law can exist, but in estimating several such control laws and comparing their predictive ability.

With larger disturbances, control of running may become notably nonlinear. Our approach can explore and evaluate the capabilities of several low-dimensional models which also show nonlinearities, and potentially come up with a model that captures empirical nonlinearities. Future work based on our methods has the potential to peer over the limits of linear predictive models, by using the rich set of tools available for nonlinear analyses.

The simplicity of the control law we identified suggests that discrete time robot controllers can create autonomous human-like runners and improved prosthetic aids to human locomotion. The question of producing a human-like linear control law based on some parameters of the individual, prosthetic or robot remains largely open. Replication of our work with additional individuals could help elucidate the factors governing such differences. In our experimental set-up, we find that ankle kinematics, which can easily be measured outside of laboratory environments using accelerometers, are particularly important in the control of running. We think that further investigation of ankle motion and associated sensory feedback could provide important insights into how humans run. However, under different conditions, other kinematic states might provide better reduced order models than those obtained with ankle data.

Our computational procedures, based on DDFA, can find broad application: the study of rhythmic phenomena in the physical sciences, and the design and control of devices that operate in periodic rather than steady state regimes.

## Funding statement

H.M.M., C.L. and A.S. were partly supported by German Research Foundation (DFG) under grant nos. SE 1042/6-1 and STR 533/7-1. S.R. was supported by ARO Young Investigator Programme award W911NF-12-1-0284. J.G. was supported by NSF grant no. 1006272.

- Received August 13, 2014.
- Accepted November 19, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.