Improving spring–mass parameter estimation in running using nonlinear regression methods

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₀), 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).

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₀ 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₀, 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₀, 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.

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⁻¹ in 0.5 m s⁻¹ 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⁻¹ and 10 steps for speeds of 3.0 and 3.5 m s⁻¹ (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.

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⁻¹. 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.

For each runner, L₀ was recorded using the two common effective leg approximations: as the height of the greater trochanter (L₀₀) relative to the ground and as estimated using the conventional 0.53 ratio of the standing height (L₀₁) (Morin et al., 2005; Winter, 2009). The α_TD for each step was estimated as per Eqn 3 (α_TD1). The traditional stiffness estimation (k₀) 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₁) were generated using the sinusoidal method proposed by Morin and colleagues (2005).

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₀, α_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.

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⁻¹ 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.

Subject characteristics and measured SLIP parameters are provided in Table 2. The two methods using measured and assumed parameters (L₀, α_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⁻¹, 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⁻¹ 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.

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.

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⁻¹ 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.

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).

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.

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₀ (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₀. 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⁻¹, whereas the aforementioned investigation used a maximal speed of 4 m s⁻¹. 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.

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).

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₀, α_TD1, L₀₀ 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.

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.

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.

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).

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₀, α_TD1, L₀₀ and t_c. Bounds were set at lower and upper limits of 5 kN m⁻¹, 63 deg, 80 cm and 0.12 s, and 30 kN m⁻¹, 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.

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₀₀ 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.

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_xi and v_xf derived from solving Eqn A19 tend to perform well across speeds and parameters sets, but the exact solutions for L₀ and α_TD may fall outside of reasonable values in some conditions (e.g. at slower speeds: <3.5 m s⁻¹), where the spring–mass dynamics are less stable. In those instances, it may be prudent to simply fix either α_TD or L₀ 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).

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.

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.