## ABSTRACT

Runners are commonly modeled as spring–mass systems, but the traditional calculations of these models rely on discrete observations during the gait cycle (e.g. maximal vertical force) and simplifying assumptions (e.g. leg length), challenging the predicative capacity and generalizability of observations. We present a method to model runners as spring–mass systems using nonlinear regression (NLR) and the full vertical ground reaction force (vGRF) time series without additional inputs and fewer traditional parameter assumptions. We derived and validated a time-dependent vGRF function characterized by four spring–mass parameters – stiffness, touchdown angle, leg length and contact time – using a sinusoidal approximation. Next, we compared the NLR-estimated spring–mass parameters with traditional calculations in runners. The mixed-effect NLR method (ME NLR) modeled the observed vGRF best (RMSE:155 N) compared with a conventional sinusoid approximation (RMSE: 230 N). Against the conventional methods, its estimations provided similar stiffness approximations (−0.2±0.6 kN m^{−1}) with moderately steeper angles (1.2±0.7 deg), longer legs (+4.2±2.3 cm) and shorter effective contact times (−12±4 ms). Together, these vGRF-driven system parameters more closely approximated the observed vertical impulses (observed: 214.8 N s; ME NLR: 209.0 N s; traditional: 223.6 N s). Finally, we generated spring–mass simulations from traditional and ME NLR parameter estimates to assess the predicative capacity of each method to model stable running systems. In 6/7 subjects, ME NLR parameters generated models that ran with equal or greater stability than traditional estimates. ME NLR modeling of the vGRF in running is therefore a useful tool to assess runners holistically as spring–mass systems with fewer measurement sources or anthropometric assumptions. Furthermore, its utility as statistical framework lends itself to more complex mixed-effects modeling to explore research questions in running.

## INTRODUCTION

Some biomechanical properties are unobservable, latent variables. That is, they represent a phenomenon, feature or behavior of the mechanical system that cannot be measured directly. In spring–mass analyses of running, the body is commonly reduced to such a mechanical system to study its elastic behavior, and the reduced parameters of this model (i.e. a runner's ‘stiffness’) are represented as latent quantities. This stiffness parameter and the spring–mass model draw their theoretical basis from the spring-loaded inverted pendulum (SLIP) model of running (Blickhan, 1989; McMahon and Cheng, 1990). This model represents the body as a mass on a single linear leg-spring (Fig. 1). It is the simplest physical system that captures the salient distinguishing features of a running body: a stance and swing phase of gait with a collision and propulsion of the mass. It is conceptually simple and faithfully describes the dynamics of running with only four deterministic parameters for a given mass (*m*) and velocity (*v*): its leg-spring length (*L*_{0}), a touchdown angle (α_{TD}), a contact time (*t*_{c}) and a spring stiffness (*k*) (Blickhan, 1989). However, the mechanical parameters of this simple system are abstractions from the human runner, and the model itself is characterized by complex nonlinear dynamics. Thus, one is required to make assumptions and estimations of model geometry to facilitate model use in experimental situations (Brughelli and Cronin, 2008).

- BIC
Bayesian information criterion

- COM
center of mass

- COP
center of pressure

*E*energy of the system

*F*_{max}maximal vertical force

*g*gravitational acceleration

*k*spring stiffness

*k*_{leg}leg stiffness

*k*_{vert}vertical stiffness

*L*_{0}leg-spring length

*L*_{00}height of the greater trochanter relative to the ground

*L*_{01}height of the greater trochanter as estimated using the conventional 0.53 ratio of the standing height

*m*mass

- ME
mixed effects

- NLR
nonlinear regression

- PS vGRF
parameterized sinusoidal estimation of the vGRF

- RMSE
root mean-squared error

- SAEM
stochastic approximation expectation maximization

- SLIP
spring-loaded inverted pendulum

*t*_{c}contact time

*T*_{step}step time

*v*velocity

- vGRF
vertical ground reaction force

- VI
vertical impulse

- Δ
*L*linear length of the leg spring

- Δ
*y*COM displacement

- α
_{TD}touchdown angle

- β
duty factor

**θ**parameter vector

*k*

_{vert}) and leg stiffness (

*k*

_{leg}) of the runner to characterize behavior of vertical displacement and leg length changes during the gait cycle under maximal force. As a dynamic SLIP model functions with a single linear elastic spring, the description here will be restricted to that of leg stiffness calculations. This stiffness is defined as the ratio of the maximal vertical force (

*F*

_{max}) to the change in the linear length of the leg spring (Δ

*L*) (McMahon and Cheng, 1990):

*F*

_{max}and

*t*

_{c}from a force plate, and α

_{TD}is approximated as (He et al., 1991):The spring length,

*L*

_{0}, of the model is the distance from the runner's center-of-mass (COM) to the foot's center-of-pressure (COP) at contact – termed the ‘virtual’ leg of the spring–mass system (Blickhan et al., 2018). However, in experimental application it is commonly assumed to be the resting leg length of the runner – termed the ‘effective’ leg of the system. It is generally approximated as the distance from the runner's hip to the foot's COP, either directly measured as the height of the greater trochanter to the ground (Brughelli and Cronin, 2008; He et al., 1991; McMahon and Cheng, 1990) or estimated as a ratio of 0.53 to the standing height of the runner (Morin et al., 2005; Winter, 1979). The vertical force recording is also used to calculate Δ

*y*, the COM displacement, via double-integration of the vertical force (Cavagna et al., 1977). The method of Morin et al. (2005) models the vertical ground reaction force (vGRF) as a sinusoid with a peak of

*F*

_{max}and a half-period of

*t*

_{c}:This method then uses measurements of contact time, flight time and running speed to estimate

*F*

_{max}and Δ

*L*. Similar to McMahon and Cheng (1990), this method assumes the leg spring length to be that of the resting leg length as previously described. Blum et al. (2009) later adapted this sinusoidal method to use the observed

*F*

_{max}and a measured value of α

_{TD}to estimate Δ

*L*rather than using the flight time.

While the methods described above have generally been found to have good agreement in their estimation of leg stiffness (Blum et al., 2009; Coleman et al., 2012; Morin et al., 2005), these approaches approximate the additional parameters of the SLIP model, which may misrepresent the fundamental elastic dynamics and spring–mass characteristics of the runner. Assigning the mass of the runner to the SLIP model is a valid assumption, but the other parameters do not translate as analogously as they are commonly assumed. When using either of the above methods, the length of the SLIP model's spring is assumed to be the leg length of the runner – the effective leg. However, to model the runner as a SLIP system, this distance should be that of the runner's COM to the point of contact on the ground – the virtual leg. A human's COM is difficult to determine statically (Clauser et al., 1969), complex to determine dynamically, and not defined by that of the human leg (Kingma et al., 1995; Maus et al., 2011; Naga, 2005; Saini et al., 1998). Clauser et al. (1969) measured the COM-to-height ratio as being 0.58, which would yield a 7–10 cm difference in *L*_{0} estimation from the traditional 0.53 leg length approximation. Similarly, Blum et al. (2009) proposed using a 5% and 10% scaling factor for COM estimation from leg length measurement from Winter (2009). A 10% difference in *L*_{0}, based solely on which COM approximation is adopted, would yield a 7% difference in *k*_{leg} estimation (Morin et al., 2005). Moreover, α_{TD}, which is the angle of the COM relative to the point of contact, is either approximated from Eqn 3 or measured as the angle of the leg at touchdown, requiring additional kinematic measurements and assumptions. Given the complexity and ambiguity of the COM location, these approximations may misrepresent the underlying spring–mass dynamics. Finally, *t*_{c} is commonly assigned as that observed in the runner. While this seems like a valid assumption on the surface, it may also inaccurately model the spring–mass behavior. The final milliseconds of propulsion are often characterized by non-linear elastic dynamics and thus exhibit a marked deviation from SLIP kinetics (Cavagna, 2006). By assigning the observed contact time to a spring–mass model, one would then bias the model towards a contact time longer than what the underlying spring–mass mechanisms would exhibit (da Rosa et al., 2019). Lipfert et al. (2012) demonstrated that by allowing *t*_{c} to vary between model and experiment, they were able to predict COM trajectories of runners more accurately from the estimated model stiffnesses. Even the assumption of the runner's average horizontal velocity introduces systemic error, as the stance velocity (from which the touchdown angle is traditionally estimated) can be 2–3% lower than the step cycle's average velocity, which ultimately challenges the predictive stability of the estimates (Bullimore and Burn, 2007; He et al., 1991).

While these differences may seem small in magnitude, they are experimentally problematic for several reasons. The first reason pertains to the inter-investigation comparability of observations. With traditional kinetic measurement methods, variations in the anatomical accuracy and method of leg length estimation and in the filtering and thresholding for contact time determination will manifest in the calculation of the touchdown angle and the corresponding leg stiffness. When kinematic methods are used for determination of these parameters, additional assumptions are introduced, further challenging direct comparison of results between studies (Brughelli and Cronin, 2008). The second reason pertains to the accuracy of the estimates in representing actual spring–mass system dynamics, necessary for creating and assessing predictive modeling of running systems (Ludwig et al., 2012; Maus et al., 2015). Stable SLIP systems and their physical realizations are highly sensitive to small parameter variations. Changes in the model's spring length by a centimeter or adjustment of the angle by a degree can demand substantial changes in stiffness or the temporal characteristics to maintain stability (Seyfarth et al., 2002). This sensitivity coupled with the approximations described above has challenged the ability to use model estimates in a predictive capacity, as the traditional estimates do not yield symmetric, stable SLIP systems (Bullimore and Burn, 2007). Because *k* is a biomechanically latent variable, and the other spring–mass model characteristics (true *L*_{0}, true α_{TD} and SLIP *t*_{c}) are difficult, if not impossible, to precisely measure in a runner, a method to estimate these parameters simultaneously from the dynamics of a SLIP system and informed by a single high-fidelity data source such as the vGRF may yield more accurate and useful descriptions of the spring–mass characteristics of a runner.

Because of the complex dynamics of SLIP models, direct model–experiment comparisons have been limited, and estimation of best-fit spring–mass parameters from the vGRF has been restricted to iterative simulations, which themselves failed to estimate parameters that produced stable running (Lipfert et al., 2012). Ludwig and colleagues (2012) proposed a method that hybridized kinetic and kinematic measures to estimate spring–mass parameters without a defined biomechanical leg using COM trajectories and best-fit SLIP model simulations. This liberated spring–mass analyses from the ambiguity of leg length and touchdown angle estimates and their corresponding assumptions. However, it requires kinematic analyses with the corresponding motion capture equipment. A method to estimate these measures with kinetic data alone would be experimentally more efficient, invite fewer sources of measurement error, and further facilitate comparison between studies and across species.

Nonlinear regression (NLR) is a candidate for an efficient, functional approach to the problem; it is a numerical method that allows for parameter estimation of a nonlinear function from observed data (Bates and Watts, 1988). If the vGRF of a SLIP system is modeled as a time-varying function, NLR is a promising tool to estimate these SLIP parameters with fewer anatomical or mechanical assumptions and no additional measurement sources. Furthermore, as opposed to alternative optimization techniques for parameter estimation, its functional framework can reveal the nature of interactions among the parameters, allowing a description of the variance structure of the estimates and extending the models to include fixed and random effects on their estimates. Leveraging this technique with a functional form of the SLIP vGRF could yield a comprehensive description of the spring–mass behavior of a runner that is informed by the full shape of the vGRF time series data rather than single extracted values.

The purpose of this investigation was fourfold: (1) to derive a parameterized, time-varying functional form of the SLIP vGRF; (2) to validate the accuracy of the function in describing the actual vGRFs of SLIP models with known parameter combinations; (3) to apply the NLR technique to analyzing the spring–mass characteristics of human runners and to compare it with the conventional methods of spring–mass analyses; and (4) to compare the stability of predicted SLIP running simulations driven by the experimentally observed spring–mass parameters from both the NLR-based and traditional estimates.

## MATERIALS AND METHODS

### A parameterized functional form of the spring–mass vGRF

*F*

_{max}and a period defined by the contact time (

*t*

_{c}) as per Eqn 4 (Blum et al., 2009; Morin et al., 2005; Robilliard and Wilson, 2005). By twice-integrating this sinusoid over

*t*

_{c}, we obtain the vertical displacement of the COM, Δ

*y*, defined by

*F*

_{max}and

*t*

_{c}(Cavagna et al., 1977; Morin et al., 2005):A maximal change in ‘leg’ length (Δ

*L*) of the model at midstance can then be derived by incorporating Eqn 5 into Eqn 2. This relationship includes the resting leg length (

*L*

_{0}) and the touchdown angle (α

_{TD}):The spring constant, or stiffness (

*k*), of the system is defined as the ratio of

*F*

_{max}to Δ

*L*as per Eqn 1 (McMahon and Cheng, 1990). Using this relationship and Eqn 6,

*F*

_{max}can be expressed as a quantity defined by

*k*, α

_{TD},

*L*

_{0}and

*t*

_{c}:Finally, to constrain the function to be non-zero from 0≤

*t*≤

*t*

_{c}, and zero-valued where

*t*>

*t*

_{c}, a logistic multiplier can be added to continuously behave as a Heaviside function:Here,

*f*(

*t*)=1 for

*t*<

*t*

_{c}and

*f*(

*t*)=0 for

*t*≳

*t*

_{c}. Combining Eqns 4, 7 and 8, a parameterized sinusoidal time-varying function is created of the SLIP model's vGRF (PS vGRF), as defined by the four parameters:

*k*,

*L*

_{0}, α

_{TD}and

*t*

_{c}:This defines each vGRF for the SLIP model as a parameter vector

**θ**, where:Therefore, a runner's experimentally observed vGRF curve can be modeled as:where

**y**is an

*n*-by-1 vector of

*n*discrete time points,

**x**is an

*n*-by-1 matrix of the

*n*discrete time points, and

**ε**is an

*n*-by-1 vector of the residual error. NLR can thereby be used to estimate the SLIP parameter θ that minimizes ε for the observed vGRF. Standard estimation approaches, such as least squares or maximum likelihood, can be used for their determination.

### Validation with SLIP simulations

To validate the accuracy of the proposed functional form of the SLIP model vGRF, the PS vGRF time series were compared with the vGRF time series of simulated SLIP models. Stable SLIP simulations were generated using the equations of motion of the sagittal plane SLIP system (Blickhan, 1989; Seyfarth et al., 2002). Models were simulated across seven running speeds from 3 to 6 m s^{−1} in 0.5 m s^{−1} increments. Models were simulated with masses of 50, 60, 70 and 80 kg. Leg lengths of 0.9 and 0.95 m were both used for 50 and 60 kg models, and leg lengths of 1.00 and 1.05 m were used for the 70 and 80 kg models, representative of human runners. Simulations were carried out in MATLAB using the ode45 solver (MATLAB 2019a, MathWorks, Natick, MA, USA) to achieve a parameter set that yielded stability over 25 steps for speeds of 4.0–6.0 m s^{−1} and 10 steps for speeds of 3.0 and 3.5 m s^{−1} (Seyfarth et al., 2002). Single-step vGRF time series were then generated for each model. A time series of vGRF data points was generated using the functional form of the SLIP GRF (PS vGRF; Eqn 9) with the simulated SLIP parameters as direct inputs, and this time series was compared against the actual SLIP model's vGRF by calculating the root mean-squared error (RMSE). Next, the four parameters for each SLIP simulation were estimated from the simulation's vGRF using NLR. The change in each parameter from the known and simulation values was recorded, and the subsequent vGRF time series with the NLR-estimated parameters was generated for each model to compare its RMSE from the SLIP simulation.

### Experimental data collection

To apply the NLR technique to human runners and examine it against traditional measurements, vGRF recordings were used from a public dataset of running biomechanics (Fukuchi et al., 2017). A detailed description of the methods is available from Fukuchi and colleagues (2017), but for the purpose of this study, select files were used from a subset of seven subjects running on an instrumented treadmill at 4.50 m s^{−1}. The vGRF was recorded continuously for 30 s at 300 Hz (FIT, Bertec, Columbus, OH, USA). All vGRF time series recordings were extracted from the database and processed in MATLAB using custom-written algorithms to isolate single step cycles with detection thresholds set at 50 N. The subject's height, weight and foot strike pattern were matched to their coded metadata file. Additionally, each subject's standing leg length was measured as the average height of the left and right legs' anatomical markers corresponding to the greater trochanter as recorded during a standing static calibration relative to the ground. Because the purpose of this investigation was to be demonstrative in nature, analyses and summary statistics for each subject are presented individually.

### Traditional estimation of spring–mass parameters

For each runner, *L*_{0} was recorded using the two common effective leg approximations: as the height of the greater trochanter (*L*_{00}) relative to the ground and as estimated using the conventional 0.53 ratio of the standing height (*L*_{01}) (Morin et al., 2005; Winter, 2009). The α_{TD} for each step was estimated as per Eqn 3 (α_{TD1}). The traditional stiffness estimation (*k*_{0}) was calculated using the kinetic method described by McMahon and Cheng (1990) with the leg compression modeled via the method of He and colleagues (1991). The conventional sinusoidal vGRFs and stiffness estimates (method 1, *k*_{1}) were generated using the sinusoidal method proposed by Morin and colleagues (2005).

### NLR estimation of spring–mass parameters

For each step cycle collected, NLR was used to estimate θ from the vGRF recordings that minimized ε. First, a constrained model was used, where only *k* was estimated, and *L*_{0}, α_{TD} and *t*_{c} were constrained to traditional measurements (method 2). Second, an unconstrained NLR was performed to estimate all four parameters for each step simultaneously, treating each step as an independent observation (method 3). Finally, a mixed-effects (ME) NLR model was estimated for each subject, treating the subject's steps as a random effect (method 4). The results of these estimations were compared with the conventional measurements and stiffness estimations. The Nonlinear Regression tools in MATLAB were used for the NLR analysis (MATLAB 2019a). See Appendix A for details of the implementation; sample code and a demonstration of the methods are available from GitHub (https://git.io/JIwzS).

To compare the quality of each method's approximated vGRF time series, all models were compared with the experimentally observed vGRF on several metrics, including the respective vertical impulses, RMSE of each model's vGRF against the observed curves, and Bayesian information criterion (BIC) of each model. This final measure provided a comparative quality assessment of a model, where a decrease in the value indicates a better fit, as the metric weights by the number of parameters to penalize overfitting.

### Stability of parameter estimations in simulated SLIP models

^{−1}using the four parameter values from both the traditional estimates and the ME NLR estimates. They were compared for stability on the basis of how many steps each model could run before becoming unstable and ‘falling’ (Lipfert et al., 2012; Seyfarth et al., 2002). The SLIP simulations were generated using the equations of motion of the sagittal plane SLIP system (Blickhan, 1989; Seyfarth et al., 2002), where the model was assigned the stiffness, touchdown angle and leg length from each estimation method. The system was assigned a horizontal flight velocity of 4.59 m s

^{−1}, which produced average step cycle velocities within 1% of 4.5 m s

^{−1}. The apex height above the touchdown height was assigned using a calculation assuming simple ballistic motion during flight. The flight time was calculated using the difference of the observed step time (

*T*

_{step}) and the observed (traditional) or estimated (NLR) contact time:The simulated spring–mass system for each subject then passively ‘ran’ until its point mass crossed the ground line, indicating a fall. The simulations were terminated if the model reached 1000 steps, which was defined as a stable model.

## RESULTS

### Simulation and validation

For each of the mass–length combinations, 5 simulations with distinct parameters were generated at each of the 7 speeds, providing 280 unique SLIP models. The parameter set for each model was used to directly generate a vGRF time series, and NLR was used to estimate parameters from each SLIP model and generate a NLR-adjusted vGRF time series. The PS vGRF function (Eqn 9) provided an accurate estimation of actual SLIP vGRFs at moderate and faster speeds on its own, and the use of NLR provided small adjustments of parameters that yielded more accurate estimation of actual SLIP vGRFs across all speeds. The RMSE of each method is summarized in Table 1, as well as the magnitude and percentage change of the NLR-adjusted parameters. The PS vGRF of Eqn 9 was an excellent estimation of the spring–mass vGRF at moderate and faster speeds (4.0 m s^{−1} and faster). NLR adjustment of the model parameters improved the fit of Eqn 9 to be excellent across all speeds, and the adjustment of these parameters was small in magnitude (0.01–0.41%). Fig. 2 provides sample vGRFs of three SLIP models at three speeds with the direct-input and NLR-adjusted PS vGRF functional estimations of the actual vGRF overlaid.

### NLR estimation of spring–mass parameters in runners

Subject characteristics and measured SLIP parameters are provided in Table 2. The two methods using measured and assumed parameters (*L*_{0}, α_{TD} and *t*_{c}) to estimate stiffness – a conventional sinusoidal approximation (method 1) and the single-parameter constrained NLR model (method 2) – yielded the poorest fits to the observed data in the vGRFs generated from the parameters (RMSE: 230.2 and 221.2 N, respectively). The full NLR method estimating all four parameters and treating each step as an independent observation (method 3) improved the fit (RMSE: 170.6 N), and the full NLR method treating each step as a random effect yielded the best fit (RMSE: 155.3 N). Correspondingly, the vertical impulses (VI) followed the same pattern, with the traditional approximation overestimating the VI (223.6 N s) compared with the observation (214.8 N s), and the ME NLR models more closely matching it (209.0 N s). The summary of parameter estimates and model fit for each method and each subject are summarized in Table 3.

The stiffness estimates from the full ME NLR model yielded values consistent with traditional kinetic approximations (−0.2±0.6 kN m^{−1}, mean±s.d.). The estimated leg lengths of the model (i.e. virtual spring–mass leg) tended to be longer than the measured leg length and conventional estimate (i.e. effective legs) by 4.2±2.3 and 4.1±2.3 cm, respectively. The touchdown angles were slightly steeper than the conventional estimates by 1.2±0.7 deg. The effective contact times of this model were 12±4 ms shorter than observed. The stiffness estimates from the full independent NLR model were consistent with the other methods (−0.2±0.6 kN m^{−1} against the conventional method). The leg lengths estimated via this approach were similar to the measured lengths (−0.7±0.4 cm), and the touchdown angles were similar to conventional approximations as well (0.4±0.1 deg). The effective contact times of this model were similarly 12±4 ms shorter than observed. The subject-specific differences among parameters are summarized in Table 4. A sample sequence of steps with vGRFs of methods 1, 3 and 4 fitted to the observed GRF is shown in Fig. 3.

Within five of the seven subjects, the ME NLR model indicated that step-to-step adjustments in leg length and touchdown angle were highly covaried (0.99), and in several subjects, adjustments in stiffness were moderately covaried with estimated leg lengths (0.03–0.84). Changes in contact time tended to be uncorrelated with adjustments in any other parameters. Correlation matrices for the four parameters for each subject are compiled in Table S1. All steps for a single subject with their individual models from the random-effects parameter estimates are shown in Fig. 4.

In addition to the improved fit of the vGRF curves, the spring–mass parameters from the fixed effects of the ME NLR method yielded more stable SLIP models that ‘ran’ further before failing in the majority of subjects as compared with simulations generated using traditional estimates. For four of the seven subjects, their NLR models outperformed the traditional models, and in one of the seven subjects, the traditional model outperformed the NLR model. In two of the subjects, both methods produced stable dynamics. The vGRF for each of subject A's simulated SLIP models can be seen in Fig. 5. The results for each subject are provided in Table 5.

## DISCUSSION

We have presented a novel method to use NLR approaches to estimate the spring–mass parameters of a runner using only the observed vGRF time series as the input. First, we derived and validated a parameterized, time-dependent functional form of the spring-loaded inverted pendulum's vGRF using the sinusoidal approximation. Then, we applied this to a group of runners and demonstrated the similarities and differences to a conventional estimation of spring–mass behavior. The NLR technique provided stiffness estimations that were consistent with traditional methods, but it more accurately modeled the runner's vGRF, more closely approximated the observed vertical impulse of the runner, and ultimately yielded spring–mass parameter sets that produced more stable simulated running. This was due to further adjustments in leg length, touchdown angle and contact time – values modeled with uncertainty via NLR but typically constrained in conventional estimates of stiffness.

### Validity of parameterized sinusoidal GRF function

The first contribution of this work was the derivation and validation of the parameterized sinusoidal vGRF function (PS vGRF). The vGRF of a runner has commonly been modeled as a sinusoid (Cross, 1999) and adapted for spring–mass analyses (Blum et al., 2009; Morin et al., 2005; Robilliard and Wilson, 2005), yet the function has never been systematically validated against the SLIP model whose behavior it was approximating. Robilliard and Wilson (2005) compared their sinusoidal approximation of the vGRF with a numerical SLIP simulation, but their analysis was restricted to a single model across several angles, and the model itself was assigned parameters simulating a horse. To our knowledge, no systematic validation of the sinusoid as a SLIP approximation has been carried out across speeds and geometries or in models representing human runners.

The PS vGRF approximation alone proved to be a valid representation of the SLIP vGRF at moderate and faster speeds (4.0 m s^{−1} and faster), and with minor adjustments of the parameters via NLR, at slower speeds. Direct input of the model parameters yielded poor approximations of the SLIP vGRF at these lower speeds. That was likely due to several features of the system. First, the SLIP model was inherently unstable at these lower speeds, with parameter configurations that were increasingly more constrained and even infeasible (Seyfarth et al., 2002; Seyfarth et al., 2003). Second, and likely most notably, the derivation of the PS vGRF's amplitude relied on the double-integration of the underlying sinusoid to approximate the vertical oscillation (Eqn 5). So, if its similarity to the SLIP model decreased at lower speeds by a fractional amount, that error would have propagated according to a power law through the integrations and yielded much greater discrepancies in the positional approximation and the corresponding sinusoidal amplitude approximation here.

By using NLR to estimate the spring–mass parameters from the ‘observed’ SLIP numerical simulation rather than directly inputting the simulation's parameters into the PS vGRF, however, we obtained parameters that were very close (all <0.5% different) to the simulation's actual values with fits that were excellent across all speeds. That was likely due to the aforementioned error propagation from the double integration, with the small parameter adjustments ‘correcting’ any underlying deviations or rounding errors. Given that the magnitudes of the adjustments were negligible and within the common reporting sensitivity of the values, the simulation comparisons suggested that the sinusoidal approximation alone was valid at moderate and faster speeds, and that NLR estimation of SLIP parameters with the PS vGRF was a valid technique to estimate spring–mass parameters across all speeds.

### NLR-estimated spring–mass parameters yield more accurate GRFs

While conventional analyses are informed by discrete kinetic data points or spatiotemporal values (e.g. maximal force and/or contact time) and are constrained by geometric approximations of the runner's anthropometry (e.g. leg length and/or touchdown angle), this method allows the spring–mass parameter approximation to be informed by full vGRF time series and allows for uncertainly in the otherwise assumed parameters. Here, these estimates resulted in vGRFs that more closely modeled those produced by the runners, with the average error (RMSE) of the ME NLR method being 155 N versus 230 N in the traditional sinusoidal method. The BIC of the NLR models, a measure of fit that penalizes the addition of parameters and overfitting, paralleled the RMSE patterns, as it decreased with each model as compared with the conventional sinusoid and was lowest in the ME NLR models (average BIC of 59,073 versus 63,204 for the ME NLR versus convention). That further resulted in VI values that more closely matched the observation: the observed and ME NLR-estimated average VI here was 214.8 N s and 209.0 N s, respectively, while the traditional estimate was 223.6 N s. The conventional methods that constrain the spring–mass estimates, especially to the explicitly observed *t*_{c}, thus yielded a parameter set that over-estimated the runner's total vertical force–time relationship.

Both Morin et al. (2005) and Blum et al. (2009) observed this limitation of the sinusoidal approximation in their respective explorations, observing VI biases of 5.3% and 10.5%, respectively. Morin et al. (2005) used a duty factor relationship to further estimate the peak vGRF (rather than explicitly constraining it), which resulted in 6.9% lower peak force, and thus the lower bias. Blum et al. (2009) mitigated the discrepancy by applying a correction factor to the sinusoid's amplitude, defined as the ratio of the observed and modeled VI, and found it to be similar to the duty factor correction. However, both of these methods simply attenuate the peak force and not the underlying temporal dynamics of the spring–mass behavior. They are thus biased by any small but significant deviations from the modeled behavior, such as the distinct nonlinear elasticity in the final moments of propulsion, where the magnitude of the vGRF forces is small, but the foot is still in contact with the ground (Cavagna, 2006; Cavagna et al., 2008). This is conceptually similar to using the ‘effective’ *t*_{c} – the time for which the runner exceeds body weight during stance – which has been shown to be more sensitive in discriminating spring–mass parameters between runners of varying ability (da Rosa et al., 2019).

### NLR-estimated spring–mass parameters differ from conventional measurements

The ME NLR modeling provided stiffness values consistent with those of the traditional estimate, but it revealed differing geometries and temporal relationships among the parameters. The hypothesis for longer leg lengths than those approximated by traditional leg length- or height-based measurements was supported. The average difference within the subjects between the traditional leg-length assumption and the ME NLR-estimated length was 4.2 cm. This suggested that the effective COM of the runners was better modeled as more distal to the ground than typically assumed. The static COM for a human standing is higher than that of the leg height, with Clauser et al. (1969) measuring it as 0.58 of the standing height. Interestingly, Blum et al. (2009) used a correction factor of 1.10 of the measured greater trochanter height (determined from subjects lying on a force plate), which would have corresponded to a ‘virtual’ spring-leg length 9.4 cm higher than that of the biological leg. Our approximations similarly predicted longer legs, and the estimation fell roughly halfway between the measured leg length and the static COM. The NLR method thus allows for subject-specific estimation of a best-fit effective COM location without the segmental assumptions, balance plates or kinematic markers otherwise required (Lafond et al., 2004; Winter, 2009).

In addition to the longer spring legs, the NLR estimation suggested that the runners tended to run with touchdown angles that were 1.2 deg steeper than the conventional estimate. Given the sensitivity of spring–mass systems to this angle, inaccurate approximation of this value would mischaracterize the system and could characterize an otherwise infeasible combination of spring–mass parameters (Seyfarth et al., 2002). The traditional method of approximating the touchdown angle (Eqn 3) necessarily underestimates a SLIP system's actual touchdown angle (see Appendix B), but it is also dependent on the assumed leg length and contact time, further confounding its accuracy. Finally, as described above, the ‘effective’ spring–mass contact times of the runners were lower than the observed values by an average of 12 ms. This would also imply a corresponding reduction in the modeled flight time and an increase in the duty factor of the runner's spring–mass system for a given step time. This temporal difference in the contact phase was likely due to the aforementioned nonlinearity in the moments prior to toe-off. Here, the magnitude of the vGRF was small relative to the rest of the time series, but the foot nevertheless remained in contact with the ground and thus extended the contact time. Clark et al. (2017) used a cosine bell-curve to capture that nonlinear elasticity fully, but the shape was not informed by the spring–mass parameters per se. Also, the selected event thresholds and filtering parameters of raw vGRF data influence the precise estimations of heel contact and toe off, and thus can significantly alter contact time estimations (Tirosh and Sparrow, 2003). The NLR method presented here resolved that sensitivity by using the entire vGRF curve to estimate the contact time of a spring–mass system that best described the systemic dynamics of the runner, rather than assigning a fixed value to that system.

### NLR modeling yields parameters that produce more stable running models

Traditional spring–mass calculations fail to produce parameters that elicit stable SLIP running models with symmetric bounces (Bullimore and Burn, 2007). This has been resolved by iterating on one or more of the calculated parameters to produce a stable, symmetric simulation, such as *k* (McMahon and Cheng, 1990), α_{TD} (Bullimore and Burn, 2007) or *L*_{0} (Lipfert et al., 2012). Bullimore and Burn (2007) identified three reasons for this model–experiment mismatch: experimental error in measured values, mechanical differences between subjects (human or otherwise) and the model, and the use of the overall mean velocity as opposed to the stance-phase velocity. Rather than making mechanical assumptions to superimpose the model on the subject, these issues may be resolved by imposing the subject on the model – that is, determining the spring–mass parameters that produce a model that behaves like the subject does. Ludwig and colleagues (2012) demonstrated a means to do this by fitting SLIP simulations to experimentally observed COM trajectories, and their method managed to produce symmetric systems over single step cycles. However, their method required both kinetic and kinematic measurements from subjects as well as iterative model simulations.

Our method attempted to resolve those three issues identified by Bullimore and Burn (2007) using only kinetic force recordings and within a statistical framework for computational and analytical efficiency. The use of a single, high-fidelity measurement source – vGRF – improved on their first point, and the modeling of uncertainty on the mechanical parameters of the system resolved the second. The derivation of the PS vGRF (Eqn 9) and the corresponding length–angle approximation (Appendix B), which approximates the flight and stance velocities, resolved the third point, as the PS vGRF equation itself was not determined by the stance velocity. Correspondingly, the NLR method produced SLIP models that demonstrated symmetric bounces and were generally as stable or more stable than models with traditional parameters. The SLIP models with NLR-estimated parameters for four of the seven subjects ran further, and two more of the seven achieved stability with both methods. As discussed previously, the parameter adjustments to realize these stability improvements were often not large in magnitude, and they were supportive of biomechanical hypotheses (e.g. slightly longer model legs). Seyfarth and colleagues (2002) demonstrated that small changes in these parameters can affect model stability, which further suggested that these mechanical parameters could be more safely modeled with uncertainty.

It was surprising that several of the models exhibited stability with the traditional method, with one outperforming the NLR-based estimates. Lipfert et al. (2012) performed a similar exercise by iteratively fitting force–length curves to experimental observations to find a best-fit *L*_{0}. However, their parameters failed to produce stable SLIP models, which they resolved by allowing the models to vary the touchdown angles from the observed estimate. We would hypothesize that our traditional models achieving stability did so for two reasons. First, we generated our simulations with horizontal flight velocities higher than the average velocity so that the simulation's step cycle velocity matched the experimental velocity, improving simulation–experiment matching. Second, we used a higher actual average running velocity itself. The parameter space for stable spring–mass running increases exponentially with speed (Seyfarth et al., 2002). The subjects used in the current study were running at 4.5 m s^{−1}, whereas the aforementioned investigation used a maximal speed of 4 m s^{−1}. It remains to be explored whether the stability of the estimates persists at lower speeds and whether the degree to which they outperform increases or decreases as the running speeds vary.

### NLR modeling facilitates efficient model–experiment comparisons of spring–mass systems in running

Another advantage of our approach is that it facilitates comparison of a runner with a spring–mass system. Because of the complexity of the system's dynamics, there is no closed-form analytical solution that describes its mechanics. Previous attempts to compare runners with the system have been limited, restricted to iterating simulation of SLIP behavior (Geyer et al., 2005; Lipfert et al., 2012). That approach is computationally intensive, preventing more complex or comprehensive analyses of many steps within a runner or a cohort. The current study demonstrated that the PS vGRF function provided a robust approximation of the SLIP vGRF. Therefore, comparison of experimental observation with best-fit curves can provide a metric of how closely a runner behaves like the spring–mass system. For example, subject C had dynamics that more closely resembled a simple SLIP model than did subject B (RMSE: 112 N versus 164 N).

### ME NLR modeling reveals correlation patterns in step-to-step spring–mass parameter adjustments

In applying the NLR technique to estimate the spring–mass parameters, the ME modeling improved the fit of the model beyond independently modeling each step. That benefit of the ME models likely stemmed from three characteristics of the modeling. First, the random effect term on each parameter came from a conditional distribution that maximized the observed likelihood function of the parameters, rather than an independent, discrete value (Bates and Watts, 1988; Feodor Nielsen, 2000). Thus, the solution was more robust to the complex interactions of the parameters and less sensitive to the starting estimates. Second, it gave the model fewer overall parameters and more degrees of freedom, as evidenced by the reduction in the BIC across models, with the ME NLR model having the lowest BIC and fewest defining parameters. Third, the sensitivity of the independent models to the starting estimates may be compounded by any errors in the solution process across the 80 steps. To minimize solver bias, we used each step's conventional parameters as starting estimates for its model (each step's *k*_{0}, α_{TD1}, *L*_{00} and *t*_{c1}), but that cannot rule out any bias or local optima from the gradient-based solution. The stochastic approximation expectation maximization (SAEM) algorithm used in the ME estimation uses global optimization tools and thus may have been more robust to this issue.

The ME NLR model provided further utility in revealing a variance–covariance structure among parameters. This allows for hypothesis testing among the parameters when fixed-effect terms are introduced. Moreover, it reveals the covariance of the parameters in the model estimates. Here, the leg length and touchdown angle were highly covaried in some of the runners, suggesting that these terms may have been better modeled together as a single geometric term (see Appendix B for suggested methodology) or with one as a fixed parameter. Furthermore, both stiffness and contact time had little covariance with the other parameters in most of the subjects, suggesting that those parameters maintained independence within the system.

The ability to characterize step-to-step parameter adjustments and their underlying variance structure presents new opportunities to study the variability patterns within an individual's gait. Current methods of characterizing gait variability assess temporal relationships of single parameters (e.g. stride length: Jordan et al., 2006; or COM excursion: Schutte et al., 2015) or compare the phasic relationships of specific joint segments (e.g. thigh–shank and shank–foot: Hafer et al., 2016). Here, we used the NLR method to assess systemic behavior with the spring–mass template across the entirety of a vGRF sequence. This provides both a means to quantify the variability in the parameters and a tractable physical realization of those adjustments. These within-subject system-level adjustments were brought into relief with the NLR analysis here, and its statistical framework provides a means to further explore their dynamics with the application of more advanced ME models.

### Limitations

While the NLR method for spring–mass analyses provided the aforementioned advantages, it had several limitations in scope and application. First, it carried several assumptions of traditional spring–mass analyses in that it was restricted to level-ground, stable-speed running. Second, the analysis was restricted to the vGRF of the SLIP system, as the sinusoidal approximation necessary for the functional analysis does not apply to the horizontal force progressions of the SLIP system. While the sinusoidal approximation of the vGRF was found to be a valid assumption across speeds and geometries, a sinusoidal approximation of the horizontal GRF would likely not map as analogously because of the greater relative pendular variation in the horizontal dynamics. In a single subject across a wide range of touchdown angles, Robilliard and Wilson (2005) demonstrated that a sinusoidal approximation of the horizontal GRF deviated from a SLIP simulation by 10–20% in the peak values in angles that would be expected from human runners. Third, when considering all four parameters, the PS vGRF does not necessarily have a unique solution with respect to spring length and touchdown angle. That can be resolved by analyzing multiple steps with a ME approach (method 4) or by using a three-parameter PS vGRF with the length–angle determination presented in Appendix B (Eqn A3). However, that only provides a measure of variance for the single ‘geometric’ parameter (Eqn A2) and not the length and angle independently. If analyzing single steps, one can also use the three-parameter PS vGRF and adopt one of the traditional approximation methods for the leg length or angle described above while still using NLR to determine best-fit stiffness and effective contact time values. Fourth, when analyzing independent steps, the method can be sensitive to starting parameters. We standardized this by using the conventional measure for each step as the starting estimate and the runner's average of those for all steps in the ME model (see Appendix A). This sensitivity was also further resolved as more steps are included in the analyses, incentivizing the researcher to analyze multiple steps per subject. Fifth, the subjects analyzed in this demonstration were relatively homogeneous in their physical characteristics, and they were observed at a single speed. While we did include runners of distinct footstrike types, we sought to control greater variation in our observed conditions to bring characteristics of the methodology into relief. This may limit the generalizability of the absolute parameter comparisons (e.g. similar stiffness or longer legs), but it establishes an effective framework to explore these and other hypotheses in future investigations.

Finally, the method was more computationally intensive than the traditional approaches. The appeal of the traditional sinusoidal method lies in its simplicity and its field-based inputs (i.e. contact time and flight time). While NLR provides a parameter set with a more accurate vGRF representation, it requires a force platform or instrumented treadmill and more demanding computation. Of note, the methods described here can be extended to model the effects of autocorrelated error often associated with force plate measurements (e.g. drift, resonance, or filtering artifacts) using approaches such as a weighted residual or directly modeling the relevant covariances (Vonesh, 1992). Furthermore, Appendix C provides a derivation of a functional approximation of the vertical displacement time series similarly characterized by the runner's spring–mass parameters. While this has yet to be explored in an experimental setting, this may present an opportunity to extend the NLR methods described here to field-based vertical displacement recordings via accelerometry and inertial measurement systems.

### Conclusion

We have presented a viable method to functionally generate a vGRF time series for a runner that is characterized by four spring–mass input variables and models the vGRF of a running SLIP system. The PS vGRF matched the vGRF of simulated SLIP models, which provided the first systematic validation of the sinusoidal vGRF as a model of the SLIP vGRF and supported its application in running. We then described the means to use that function with NLR to estimate a runner's stiffness, as well as his or her virtual leg length, touchdown angle and contact time, using only the observed vGRF time series; that is, using the holistic shape of the runner's vertical force curve to inform the spring–mass model estimates. This liberates stiffness estimates from assumed geometric and temporal constraints, facilitating more efficient and better-fitting spring–mass approximations. When used across many steps with a mixed model, the NLR technique yielded stiffness estimations that were consistent with traditional estimates. The virtual leg length approximations were longer than the traditional leg length measurements, and the effective contact times were shorter than the observed values. Together, these NLR-estimated spring–mass characteristics yielded vGRFs that more closely simulated the observed vGRFs. In addition to its fidelity as an analytical technique, this method has broad application in modeling more complex research questions with both fixed and random effects, such as including multiple runners in a multilevel ME NLR regression that would allow testing of various runner covariates and cohort parameter differences.

### APPENDIX A

#### Details related to NLR implementation

All analyses were conducted in MATLAB using the Nonlinear Regression tools within the Statistics and Machine Learning toolbox (MATLAB 2019a). Sample code and a demonstration of the methods are available from GitHub (https://git.io/JIwzS).

#### Starting values

Custom initial parameter values for the NLR models were assigned to each subject as a parameter set that minimized the sum of squared errors against all steps together via nonlinear least-squares optimization seeded with values corresponding to the mean value of the subject's conventional parameters: *k*_{0}, α_{TD1}, *L*_{00} and *t*_{c}. Bounds were set at lower and upper limits of 5 kN m^{−1}, 63 deg, 80 cm and 0.12 s, and 30 kN m^{−1}, 74.5 deg, 120 cm and 0.40 s for the four parameters, respectively. These were applied to allow a liberal range of possible values that could be expected from a human runner at this speed (Seyfarth et al., 2002), while also restricting the solver from exploring unreasonable or infeasible parameters.

#### Parameter constraints

^{10}versus 10

^{30}). Note that these logistic terms are distinct from the multiplier used in Eqn 9 to ensure that the modeled force was zero for

*t*>

*t*

_{c}.

#### ME model estimation

For the full ME NLR parameter estimation, the SAEM algorithm was used to estimate a random-effects model for each subject, treating an individual's steps as random effects (Feodor Nielsen, 2000). Each subject's model was seeded with an initial random-effects covariance matrix with the diagonal equal to 2 times the observed variance in the conventional parameters across steps, with the variance in *L*_{00} as estimated by the variance in α_{TD0} relative to the mean *t*_{c}. Five burn-in iterations were used, and five chains were simulated. For each of the three phases of the SAEM algorithm, 300, 200 and 100 iterations were performed, respectively.

### APPENDIX B

#### A singular length–angle approximation

*A*such that:To enable use of NLR to estimate spring–mass parameters for isolated steps, this term can be incorporated into Eqn 9 as a three-parameter model:Using NLR, this method can provide an estimate for

*A*. With that estimation, properties of the spring-loaded inverted pendulum's dynamics with a horizontal velocity approximation can be used to isolate

*L*

_{0}or α

_{TD}and solve for their values with

*A*.

*y*

_{0}is the vertical position at landing and takeoff,

*y*

_{f}is the peak vertical excursion during flight, and

*v*

_{x}_{i}is the horizontal velocity during flight and at the initial point of stance:The energy of the system can also be considered at midstance where the COM is at its lowest point. Here, Δ

*y*is the vertical oscillation during stance defined by Eqn 5, Δ

*L*is the change in leg length during stance as defined by Eqn 6, and

*v*

_{x}_{f}is the horizontal velocity at midstance:Without use of additional measurement equipment (i.e. markers or accelerometers), the exact values of the horizontal velocity throughout the gait cycle are unknown. The common approximation of the relationship of the leg length to the touchdown angle as described in Eqn 3 is inappropriate, as it will underestimate the values given that the average velocity during stance will be necessarily less than the average velocity of the full gait cycle. If the change in horizontal velocity is approximated during stance as a linear decrease from

*v*

_{x}_{i}to

*v*

_{x}_{f}, the velocity values can be related with the duty factor, β:

*A*known, Eqns A2 and A19 can therefore be solved for unique identities of

*L*

_{0}and α

_{TD}. Eqn A19 can be expressed with known quantities:

**,**

*g**m*,

*v*,

*k*,

*A*and

*t*

_{c}, and solved numerically for

*L*

_{0}. With

*L*

_{0}, Eqn A2 can be used to then solve for α

_{TD.}

As a caution, it should be noted that Eqns A4–A7 rely on data derived from accelerations and are themselves approximations of velocities. Therefore, any errors or deviations from absolute dynamics are propagated when integrating to displacement values for the system. So, the approximations for *v _{x}*

_{i}and

*v*

_{x}_{f}derived from solving Eqn A19 tend to perform well across speeds and parameters sets, but the exact solutions for

*L*

_{0}and α

_{TD}may fall outside of reasonable values in some conditions (e.g. at slower speeds: <3.5 m s

^{−1}), where the spring–mass dynamics are less stable. In those instances, it may be prudent to simply fix either α

_{TD}or

*L*

_{0}with a conventional assumption (e.g. Eqn 3 or a height approximation), and fit a three-parameter model for use with independent, single steps.

A MATLAB function to implement this approximation is available with the code from GitHub (https://git.io/JIwzS).

### APPENDIX C

#### Application of NLR to COM trajectories

The methods described in the current study are for the use of NLR applied to vGRF recordings. This constrains the technique to experimental conditions where force recordings are possible (e.g. running over a force plate or on an instrumented treadmill). To facilitate exploration of these techniques in field-based settings, it may be possible to use the runner's COM displacement time series as captured by accelerometers or inertial measurement units.

*y*

_{0}) taken as

*L*

_{0}sinα

_{TD}, a functional form of the vertical displacement during stance as defined by the four spring–mass parameters can be written as:To facilitate use of NLR, the logistic function from Eqn 8 can be modified and used, or a state-shift can be employed to describe the switch to simple ballistic motion after toe-off:Combining these, a continuous description of the runner's vertical displacement can be characterized through stance and flight:This functional form of the vertical displacement time series can be applied to kinematic or accelerometry-derived recordings of runners and may enable application of NLR methods. This would allow for estimation of spring–mass parameters using the computational methods described herein without the requirement of force recordings. While this is a presentation of the theoretical functional form required for NLR, the utility of its application with field-based equipment and measurement systems has yet to be explored.

## Acknowledgements

We would like to thank the two anonymous scientists who reviewed this article. Their suggestions and critiques improved the early iterations of the manuscript, and we are grateful for the time, insights and intellect they contributed to the process.

## Footnotes

**Author contributions**

Conceptualization: G.T.B., R.G.; Software: G.T.B.; Validation: G.T.B.; Formal analysis: G.T.B., R.G.; Investigation: G.T.B., R.G., R.F.Z.; Writing - original draft: G.T.B., R.G.; Project administration: R.F.Z.

**Funding**

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

**Data availability**

Sample data and code for a demonstration of the NLR implementation in MATLAB are available from GitHub: https://git.io/JIwzS.

## References

*Nonlinear Regression Analysis and Its Applications*

*J. Biomech.*

*J. Exp. Biol.*

*J. Biomech.*

*Scand. J. Med. Sci. Sports*

*J. Theor. Biol.*

*J. Appl. Physiol.*

*J. Exp. Biol.*

*Am. J. Physiol.*

*J. Exp. Biol.*

*J. Exp. Biol.*

*Weight, Volume, and Center of Mass of Segments of the Human Body*

*J. Biomech.*

*Am. J. Phys.*

*Front. Physiol.*

*Bernoulli*

*PeerJ*

*J. Theor. Biol.*

*J. Sports Sci.*

*J. Appl. Physiol.*

*Gait Posture*

*J. Biomech.*

*J. Biomech.*

*J. Theor. Biol.*

*J. Biomech.*

*J. Exp. Biol.*

*J. R. Soc. Interface*

*J. Biomechanics*

*J. Appl. Biomech.*

*J. Exp. Biol.*

*J. Biomech. Eng.*

*PLoS ONE*

*J. Biomech.*

*J. Exp. Biol.*

*J. Appl. Biomech.*

*Stat. Med.*

*Biomechanics and Motor Control of Human Movement*

**Competing interests**

The authors declare no competing or financial interests.