A simple computational method to estimate stance velocity in running

While the approximation in Eqn 1 is appealing in its simplicity and efficiency, its consistent underestimation of the angle may be problematic for three reasons. The first is in propagation of that error through subsequent spring–mass calculations. As leg stiffness estimation relies on that angle (or analogous stance distance), inaccuracies in the estimation will manifest in inaccuracies in estimations of leg-spring compression and leg-spring stiffness. Thus, an underestimation of impact angle will overestimate leg compression, which will subsequently underestimate stiffness (Lipfert et al., 2012; Morin et al., 2005). It is unknown the extent to which these deviations vary with speed. That leads to the second issue, which is that of experimental generalizability. Inaccuracies in measurements will make direct comparison of results (i.e. leg stiffness) challenging between different measurement systems, where kinetic recordings with the above assumptions are systematically biased against other measurement systems, such as kinematic-based (e.g. motion capture) estimations. No two methods will be directly analogous, but additional sources of bias will inevitably obscure interpretations of findings between studies. Finally, the traditional angle approximation may prevent use of the spring–mass model in a predictive fashion. Experimentally derived spring–mass parameters have conventionally failed to produce stable spring–mass model simulations (Lipfert et al., 2012; Ludwig et al., 2012). Bullimore and Burn (2007) hypothesized that the common use of the average gait-cycle velocity in the parameter estimation was one of the primary reasons for the apparent instability of observed spring–mass estimations in runners.

A method to more accurately approximate the stance velocity and the subsequent impact angle for spring–mass running that does not require kinematic recordings may be possible using the energetic relationships of the spring-loaded inverted pendulum (SLIP) dynamics. The first hypothesis of this investigation was that a numerical relationship could be derived that estimates a runner's stance velocity using only the observed contact time and flight time with the average gait cycle velocity and runner anthropometry. The second hypothesis was that this approximation of stance velocity would more accurately approximate the impact angle and leg stiffness in simulated SLIP models across a variety of speeds and sizes. The third hypothesis was that the new stance velocity estimates could be used to estimate the peak horizontal braking and propulsive forces. Morin et al. (2005) demonstrated the accuracy of modeling the vertical ground reaction force of a runner using only the contact time and flight time, but no such simple analog has been demonstrated for the horizontal forces. The fourth and final hypothesis was that the differences between the new method and the conventional methods that were observed in analyzing SLIP simulations would be analogous in human runners, suggesting a more accurate characterization of the runner's spring–mass system.

To assess the accuracy of the proposed method to approximate stance velocity and impact angle in spring–mass systems, the approximations and the conventional angle approximation to simulated SLIP models were compared. One hundred model systems representative of the size and speed of typical runners were generated. The models ranged from 50 to 80 kg with leg-spring lengths of 90–105 cm. Each model was run at 3.5, 4.0, 4.5, 5.0 and 5.5 m s⁻¹ for a total of 500 simulations. Particle-swarm optimizations were used to find stiffness and impact angle values that yielded stable, symmetric systems, where stability was assessed as a system that could exceed 25 passive steps without failure at 4.0–5.5 m s⁻¹ and 10 steps at 3.5 m s⁻¹ (Seyfarth et al., 2002). All models were simulated using the ode45 solver in MatLab (MatLab 2019a, MathWorks, Natick, MA, USA) with the SLIP model equations of motion as described by Eqns 2–6 (Blickhan, 1989).

For each simulation, the initial, final and average stance velocities were recorded for each of the SLIP models along with its actual impact angle, leg stiffness and peak hGRF. Using its contact time and flight time, the model's initial, final and average stance velocities were then approximated, along with the corresponding impact angle, leg stiffness and hGRF with the new method proposed above. For comparison, the impact angle and leg stiffness were also calculated for each simulation model using the conventional approximation from Eqn 1 (He et al., 1991).

The conventional impact angle measurement was compared with the new proposed method using the vertical ground reaction force (vGRF) recordings of 28 runners from a public dataset of running biomechanics (mean±s.d. age 34.8±6.7 years). The methods are described in detail by Fukuchi and colleagues (2017). Briefly, the runners ran on an instrumented treadmill (FIT, Bertec, Columbus, OH, USA) at 2.5, 3.5 and 4.5 m s⁻¹, and the ground reaction forces were recorded continuously for 30 s at 300 Hz. For the current study, the vGRF time series recordings were extracted from the database and processed in MatLab using custom algorithms, where single step cycles were isolated with detection thresholds set at 50 N (84±7 steps per subject per speed; 7091 steps total). Each subject's mass was extracted from the metadata (69.6±7.7 kg), and his or her leg length was determined from the average height of anatomical markers on the left and right legs corresponding to the greater trochanter as recorded during a standing static calibration relative to the ground (leg length 92.8±5.0 cm). The contact time and flight time for each step cycle were then recorded for use in the new stance velocity and impact angle calculation, and the conventional impact angle was calculated using Eqn 1. From each of these angle estimates, leg stiffness was calculated from the vGRF recordings as per Eqns 17 and 18 using the method of McMahon and Cheng (1990), where the CoM displacement was determined via double integration of the vertical acceleration (Cavagna, 1975). Peak braking and propulsive hGRF measures were also recorded for each step, and hGRF estimations using the new method were similarly calculated.

To test the accuracy of the velocity and angle estimates provided by the proposed method in simulated spring–mass SLIP model running systems, differences were calculated between the method-estimated average stance velocity and the actual average stance velocity for each simulation. Differences were also calculated between each method's estimated impact angle and leg stiffness and each model's actual impact angle and actual stiffness for the two estimates (i.e. new method and estimate) and the differences between the new and traditional method. Finally, peak hGRF estimates were compared from the new method with the model's actual peak hGRF. A linear regression was conducted on each criterion variable to examine the fixed effects of speed and model anthropometry on the magnitude of any estimate discrepancies. The criterion variable was the respective difference (i.e. average stance velocity, impact angle, leg stiffness and peak hGRF), and the predictors were speed, mass and leg-spring length, respectively. Correlation coefficients were also calculated for agreement between estimates and actual values for descriptive purposes.

To assess the performance of the estimates in the human runners, a similar analysis was conducted, where differences between the new method and the traditional method were calculated at each speed for impact angle and leg stiffness. Then, a linear mixed-effects regression analysis was conducted using the difference between these two methods as the criterion variable. Speed and leg length were treated as fixed effects, and the subjects were modeled with random intercepts. For the horizontal force measurements, a similar analysis was conducted using the difference between the estimated and observed peak hGRF in both braking and propulsion. All analyses were assessed with a Type I error control rate of P<0.05, and P-values were adjusted for multiple comparisons within each family of hypotheses (i.e. simulations and experimental analyses) using the Benjamini–Hochberg method for false discovery rate control (Benjamini and Hochberg, 1995). All model distributions were verified for normality via their residual and quantile–quantile plots. All data analyses were conducted in R v4.0.2 (http://www.R-project.org/).

Across speeds and models, the stance velocity estimate demonstrated near-perfect agreement with the simulation models' actual stance velocity (r>0.99). The mean difference between the estimated and actual velocities within models was −0.009 m s⁻¹ or −0.2% (95% confidence interval, CI: −0.008, −0.009 m s⁻¹; limits of agreement: −0.023, +0.005 m s⁻¹). These relationships are shown in Fig. 2. The actual average stance velocity for all models was 0.06 m s⁻¹ lower (1.35%) than the gait cycle average velocity (‘average running speed’). While small in relative magnitude, the new method's bias was related to both speed and leg length, where it decreased slightly at faster speeds and increased slightly at slower speeds. There was also a small interaction of these effects. The magnitudes of these relationships are provided in Table 1 and visualized in Fig. 2.

The new method's stance velocity estimate yielded more accurate impact angle estimates across all speeds and models. Bias for the new method was −0.05±0.04 deg, whereas the conventional angle approximation carried an average bias of −0.36±0.03 deg. Correspondingly, within models, the new method estimated the impact angle to be 0.31±0.03 deg steeper than with the traditional method. In both methods, the bias from the actual angle increased in magnitude at faster speeds and with shorter legs with a further significant interaction between the two. The difference between the methods followed that same pattern, where it was augmented at faster speeds and in models with shorter legs. Fig. 2 compares the new method against the simulation's known values, and Fig. 3 demonstrates the relationships between the new and traditional estimates. Table 2 describes the linear models of the bias for each method.

The new method also more accurately estimated the stiffness of each of the models across speeds and geometries. The conventional method underestimated the stiffness by 218±36 N m⁻¹ (0.34±0.06 BW/L₀), whereas the new method only underestimated it by 32±25 N m⁻¹ (0.05±0.04 BW/L₀). Within models, the average difference between the two methods was 186±19 N m⁻¹ (0.29±0.02 BW/L₀). Analogous to the angle estimates, the bias of each method was related to both speed and leg length with significant interactions, though the magnitude of these relationships was also small relative to the absolute values. Similarly, Fig. 2 presents the new method against the known simulated values, and Fig. 3 visualizes the new and traditional estimates. Table 3 also provides analyses of the bias.

The horizontal force estimates from the new method were highly correlated with the true simulation model forces (hGRF: r=0.82, BW: r=0.73). The average bias was −27.1 N (−0.05 BW, −8.5%), with 95% limits of agreement of −91.4–37.3 N (−0.15–0.06 BW). Similar to the other estimates, the bias was significantly related to the speed and leg length of the models as well as to the interaction of these effects. The magnitude of these relationships was more substantial, with the bias decreasing by 42 N (−0.07 BW, −13%) for every 1 m s⁻¹ increase in speed, and increasing by 1.4 N (0.004 BW, 0.8%) for every 1 cm increase in leg length. This is visualized in Fig. 4, and the results are summarized in Table 4.

In human runners, the difference in impact angle estimates between the new method and the conventional method followed the same patterns observed in the simulations. The new method estimated angles that were 0.40±0.07 deg steeper than the conventional estimate. Similar to the simulations, that difference had a small but significant relationship to both speed and leg length, where it decreased at faster speeds and increased slightly with longer legs. The difference between the leg stiffness estimates from the two angle estimation methods also followed the simulation patterns. However, the magnitudes of the differences were larger and distributed non-normally. The strong rightward kurtosis suggested a logarithmic distribution, so we used a log₁₀ transformation on the differences prior to analysis. The average magnitude of the differences after transformation was 0.59 kN m⁻¹. These results are detailed in Table 5 and presented in Fig. 5.

The hGRF estimation did not perform as well on the human runners as it did on the simulation models. In both braking and propulsion, the horizontal force estimates were inconsistent, with correlations of r=0.33 and r=0.50, respectively. The average bias and standard deviation at 2.5, 3.5 and 4.5 m s⁻¹ in braking were +72±40 N (+0.10 BW, 38%), −11±54 N (−0.02 BW, −1%) and −75±80 N (−0.11 BW, −17%). In propulsion, the average bias at each of the three speeds was +123±32 N (+0.18 BW, 82%), +53±31 N (+0.08 BW, 24%) and +13±37 N (+0.02 BW, 5%). Similarly, these results are described in Table 5 and visualized in Fig. 5.

Using the energetic states of a simple spring–mass running system, we present a method to approximate stance velocity using the contact times and flight times of a runner. The new method provided accurate average stance velocities across a range of running speeds. That estimation further provided more accurate estimations of impact angle and leg stiffness in spring–mass models. The precision of the velocity estimates, however, did not translate as analogously to the estimations of horizontal acceleration and corresponding hGRF patterns in the simulations. When applied to human runners, the new method further provided values for the impact angle and leg stiffness that were similarly distinct from traditional estimates as compared with the differences observed in the simulations, suggesting the new method may provide accurate stance velocity estimations and more accurate spring–mass parameters without the need for kinematic recording systems.

The accuracy of the method on simulated spring–mass systems was tested where the exact stance velocity, impact angle and leg stiffness were known. The first two hypotheses were supported, as the method provided accurate estimations of stance velocity, and those velocity estimations further provided more accurate impact angle and leg stiffness approximations in the models. The new method's estimations were highly consistent with the known values of the simulations (r>0.99 for all), and the bias was small in relative magnitude (all <0.4%). Across all measures, the bias was significantly impacted by the speed and leg length of the running model. For example, in the stance velocity estimates, the bias was negligible at 3.5 m s⁻¹, but increased modestly across speeds to −0.017 m s⁻¹ at 5.5 m s⁻¹ (−0.31%). Longer-legged models had smaller bias compared with the shorter models but, similarly, the bias was small in relative magnitude (e.g. −0.29% for 90 cm models versus −0.11% for 105 cm models). When assessed as dimensionless quantities (e.g. against the Froude number instead of absolute speed), similar small factor-dependent bias was present. Therefore, it is likely that these small dependencies in both the traditional methods and the new estimation may be related to the nuanced implications of the model assumptions, such as the linear decrease in velocity through stance.

From a practical perspective, the magnitudes of the speed and length dependencies are negligible, as the amount that the bias changed across speed and leg length was consistently and considerably smaller than the bias from using traditional methods to approximate stance velocity, impact angle or leg stiffness (Fig. 3). Again, considering the stance velocity estimate, the traditional method that assumes the average gait cycle velocity overestimated the stance velocity by 0.06 m s⁻¹, or 1.35% (Fig. 3 and Table 1), with similar speed and length dependencies. The new method only overestimated it by 0.009 m s⁻¹, or 0.18% (Fig. 3 and Table 1).

The hypothesis pertaining to the horizontal force estimations was not entirely supported. In the aggregate, the consistency of the new method's estimate among the simulations was good, with correlations of 0.82 and 0.73 for force and acceleration values, respectively. The bias indicated a modest underestimation of the force (−27 N, 0.05 BW, −8.5%), but the performance of the estimate and the magnitude of the bias were strongly influenced by the speed of the model. Across the five speeds, the bias increased from +15 N (+6.3%) at 3.5 m s⁻¹ to −69 N (−19%) at 5.5 m s⁻¹.

This is likely due in part to the difficulty in approximating the horizontal force progression of a SLIP system with a sinusoid. The vertical force progression of a spring–mass system has been accurately approximated as a sinusoid across a range of speeds and model geometries (Burns et al., 2021), but no similar investigation has tested that approximation across a breadth of distinct models for the horizontal force progressions. Robilliard and Wilson (2005) explored their sinusoidal approximation with several numerical SLIP simulations across a range of impact angles in a single model, and the horizontal acceleration predictions were similarly challenged.

The SLIP hGRF progression is not perfectly symmetric within the braking and propulsive phases. The influence of the spring in the horizontal deceleration of the mass is at its peak at initial contact, but it falls to zero as it pendularly rotates to a vertical orientation at midstance. As such, the peak braking force in a SLIP system occurs prior to the midpoint of the braking phase, as opposed to a sinusoidal progression, which would result in the peak occurring exactly at the midpoint. This earlier peak (and correspondingly later peak in the propulsive phase) will likely result in a modestly greater amplitude than in a symmetric sinusoid as well. That trend towards an under-approximation of the peak hGRF via the sinusoid was observed, as all speeds associated with robustly stable SLIP systems (>4.0 m s⁻¹) elicited lower estimations of the hGRF. From a practical perspective, this may challenge the direct implementation of this method for hGRF estimates. However, the bias was strongly predicted by the model speed and geometry, which suggested an opportunity to use a corrective algorithm to better predict the peak hGRF values in practice, which is described later.

The steeper touchdown angle and higher leg stiffness values of the new method's estimations were also observed in humans. The magnitude of the angle bias against the traditional estimation was similar in the human runners to that observed in the simulations: +0.33 deg in humans versus +0.40 deg in simulations at 3.5 m s⁻¹. In the simulation analyses, the steeper angles and higher leg stiffnesses from the new method were closer to the true spring–mass system values, which would support the notion that those same patterns in the human estimations may be better modeling the ‘actual’ spring–mass characteristics of the runner via more accurate stance velocity approximations. That was expected, as adopting the gait cycle average velocity instead of the stance velocity will necessarily underestimate the impact angle and correspondingly overestimate the stiffness.

While the impact angle and leg stiffness estimations in humans corresponded to the trends observed in the simulations, the horizontal force estimations were not as consistent. The average bias at 3.5 m s⁻¹ was only −0.02 BW in braking, but the spread was large (i.e. standard deviation of 0.08 BW). However, this overestimated the actual braking hGRF at slower speeds and underestimated it at faster speeds. In propulsion, it was most accurate at the fastest speed, with an average bias of +0.02 BW, though the spread was similarly large with a standard deviation of 0.05 BW. These discrepancies in accuracy in braking and propulsion were likely due to three factors: (1) asymmetries in the energetic balance through stance (Cavagna, 2006), (2) variability in braking patterns within human runners (Munro et al., 1987), and (3) the instability of SLIP model systems and their parameters at slow speeds (i.e. 2.5 m s⁻¹) (Seyfarth et al., 2002).

In the simulated models, the bias in the peak horizontal force estimations was significantly predicted by both speed and leg length, with the relationship being most pronounced in the acceleration estimates (Fig. 4). This presented an opportunity to test a correction transformation on the human force estimates using the linear model coefficients from the simulations as provided in Table 4. By transforming the force estimations with corrections for speed (−0.065 BW per m s⁻¹ with an intercept of −0.046 BW at 4.42 m s⁻¹), the force estimations improved across speeds (r=0.69 and 0.85 for braking and propulsive forces, respectively), and the speed dependency of the bias improved. An additional leg length correction did not substantively improve the estimate, and the individual variation following the speed correction remained large (Fig. 6). This suggests the method may need further refinement for consistent prediction of horizontal forces in running bodies.

One area to realize this refinement in the estimates for running subjects may be in the initial leg length assumption. The spring length in a runner modeled as a SLIP system is that of the CoM to the CoP. We adopted the common assumption that the resting length L₀ corresponded to the distance of the runner's greater trochanter (GrT) to the CoP while standing. In humans, the CoM is higher than the GrT while standing, consequently underestimating the CoM–CoP distance. This discrepancy is justified by the notion that a runner lands with a flexed knee at contact, suggesting that the resting GrT–CoP is a reasonable estimate for the initial CoM–CoP distance while running (Bullimore and Burn, 2007). However, a longer spring will necessitate steeper impact angles and higher leg stiffnesses, and the choice of this definition will affect experimental values accordingly (Müller et al., 2016). Blum et al. (2009) proposed using a sex-specific correction factor of 1.05 and 1.10 for women and men, respectively, to reconcile the difference of the GrT and CoM lengths, and Burns et al. (2021) found that a nonlinear regression estimation of L₀ from a vGRF time series approximated the distance as 1.05 times that of the GrT measurement. In the analyses of the simulation models here, the resting spring length as an input was known. Consequently, none of the observed variation in the simulation estimates was due to error in that approximation and, correspondingly, the absolute differences in the estimates and the true values were quite small. However, in the human runners, this quantity was unknown and approximated via the aforementioned GrT measurement. Therefore, any discrepancy in these anthropometric assumptions may have contributed to the greater individual variation in both the stiffness estimates and the hGRF approximations in human runners.

The new method accurately estimated stance velocity and impact angle in running using only measurements of contact time, flight time and average running speed, with mass and leg length as inputs. Thus, it can be easily implemented in experimental settings where high-resolution kinematic recordings are infeasible or unavailable. That includes common kinetic investigations using overground runways with timing gates or instrumented treadmills where the belt speed is used. Furthermore, it can be used in field-based investigations where average running speed is captured from a GPS watch and temporal measures from accelerometers, inertial measurement units (IMUs) or portable high-speed cameras. This complements the work of Morin and colleagues (2005), who presented means to estimate peak vertical force, vertical oscillation, and vertical and leg stiffness using these same inputs. Furthermore, it built on the work of Robilliard and Wilson (2005), who presented a means to approximate spring–mass dynamics with temporal and horizontal velocity measurements. The increase in accessible and accurate temporal recording technology, such as commercial IMUs and smartphone-integrated high-speed cameras (Balsalobre-Fernández et al., 2017), creates new opportunities to integrate and apply these methods for more accurate, generalizable spring–mass analyses in a variety of ecological settings.

The new technique for stance velocity and angle approximations was limited by its requirement of a numerical equation solver. It can be easily implemented in any numerical computing language (e.g. MatLab, R or Python), but that requires familiarity with one of those tools. Both a MatLab and an R function have been provided (see Supplementary Materials and Methods or https://git.io/JmxRB) to implement the new method, where the stance velocity, impact angle, hGRF estimate and corrected hGRF estimate are provided. If a researcher or practitioner is not proficient with one of these tools, a simple correction for the stance velocity would provide an improvement over the conventional use of the average step cycle velocity. Here, the simulation models' stance velocity was 98.65±0.09% of their average step cycle velocity, so a correction of 0.987 would be a reasonable adjustment in the absence of a numerical equation solver. This, however, will inevitably be less accurate than the method's estimate given the subject- and speed-specific energetic relationships, so the numerical estimation is strongly recommended. Furthermore, the method's derivation assumes a level ground through the gait cycle, restricting its application to flat-terrain running. An extension of the estimate to include energetic changes related to changes in ground height between steps may be useful, as runners modulate temporal parameters and spring–mass characteristics as they navigate gradients and uneven terrain (Iversen and McMahon, 1992; Müller and Blickhan, 2010; Müller et al., 2012). Finally, the application of the method for hGRF estimations should be approached with caution. The aforementioned correction improves the estimation, but the estimate still exhibits enough variance to challenge its use for contexts where a high sensitivity in hGRF measurement is required.

A novel method was developed and presented to estimate the average stance velocity in runners using only their contact time, flight time, average speed, mass and leg length or height as inputs. This allows for accurate estimation of stance velocity without kinematic recordings, facilitating a more accurate description of spring–mass characteristics of runners. The new method accurately estimated the stance velocity of SLIP running system simulations across a range of speeds and sizes, and its corresponding impact angle and leg stiffness estimates were closer to each system's actual angle than the conventional approximation. This further resulted in more accurate leg stiffness estimations. It was hypothesized that the average stance velocity could be used to estimate the peak hGRF forces in a spring–mass system, but that hypothesis was not strongly supported, as there was modest agreement between the estimates and the actual values but large speed- and length-dependent variations in the estimate's bias. In human runners, the estimates for the impact angle and leg stiffness followed the same patterns as those observed in the simulations when compared against the traditional estimates, suggesting more accurate modeling of the runners' spring–mass behavior. The new method can be applied to laboratory and field-based investigations alike to accurately estimate stance velocity and improve spring–mass analyses.

The authors would like to thank all the researchers – many of whom can be found in the in the References section – whose work informed this project. They would also like to thank the two anonymous scientists who provided a review of the manuscript.