Running performance, energy requirements and musculoskeletal stresses are directly related to the action–reaction forces between the limb and the ground. For human runners, the force–time patterns from individual footfalls can vary considerably across speed, foot-strike and footwear conditions. Here, we used four human footfalls with distinctly different vertical force–time waveform patterns to evaluate whether a basic mechanical model might explain all of them. Our model partitions the body's total mass (1.0*M*_{b}) into two invariant mass fractions (lower limb=0.08, remaining body mass=0.92) and allows the instantaneous collisional velocities of the former to vary. The best fits achieved (*R*^{2} range=0.95–0.98, mean=0.97±0.01) indicate that the model is capable of accounting for nearly all of the variability observed in the four waveform types tested: barefoot jog, rear-foot strike run, fore-foot strike run and fore-foot strike sprint. We conclude that different running ground reaction force–time patterns may have the same mechanical basis.

## INTRODUCTION

The bodily motion of terrestrial animals that use bouncing gaits is determined by the action–reaction forces between the limbs and the ground. However, the predominant orientation of these forces during straight-path, level running and hopping is not in the horizontal direction of travel (Cavagna et al., 1977). Horizontal force requirements are minimized by an effective step-to-step maintenance of forward momentum once an animal is up to speed. Vertical force requirements, in contrast, can exceed body weight by a factor of two or more during periods of limb–ground contact (Weyand et al., 2000). Large vertical forces result from two factors: the need for stride-averaged vertical forces to equal the body's weight, and limb–ground contact periods that comprise only a fraction of the total stride time. Consequently, the vertically oriented ground reaction forces exceed horizontal forces by a factor of five or more, and lateral forces by greater margins.

The vertical force versus time waveforms of individual running and hopping footfalls can vary considerably in duration, amplitude and shape. This variation has been documented for a variety of species (Cavagna et al., 1977) and most comprehensively for humans (Bobbert et al., 1991; Munro et al., 1987). At present, several factors are known to introduce the shape variation that occurs predominantly in the initial portion of these force–time waveforms. These include: running speed (Bobbert et al., 1991; Kuitunen et al., 2002; Munro et al., 1987; Weyand et al., 2009; Weyand et al., 2010), the portion of the foot that initially contacts the running surface (Cavanagh, 1987; Chi and Schmitt, 2005; Dickinson et al., 1985; Ker et al., 1989; Lieberman et al., 2010; Nigg et al., 1987) and footwear (Liu and Nigg, 2000; Ly et al., 2010; Nigg et al., 1987; Nigg and Liu, 1999; Zadpoor and Nikooyan, 2010). Current understanding rests heavily on the two types of models most frequently used to interpret these waveforms: the spring-mass model and multi-mass models. Models of both types are well-founded and have undergone extensive evaluation. However, neither was formulated to explain these waveforms in full.

The most basic treatment of the vertical force–time waveforms is provided by the classic spring-mass model (Blickhan, 1989; McMahon and Cheng, 1990). The single-mass approach models running and hopping animals as a lumped point-mass mass bouncing on a massless leg spring. This single-mass model explains many aspects of running and hopping gaits with remarkable accuracy given its mechanical simplicity (Bullimore and Burn, 2007; Farley et al., 1993; Ferris and Farley, 1997; McMahon and Cheng, 1990). However, this classic model was formulated largely for broad evaluative purposes, not specific quantitative ones. Accordingly, the perfectly symmetrical force–time waveforms the model predicts (Bullimore and Burn, 2007; Robilliard and Wilson, 2005) cannot account for the non-symmetrical components that the force–time waveforms inevitably contain. These include, but are not limited to, heel-strike impacts at slow speeds and extremely rapid rising edges at faster ones (Kuitunen et al., 2002; Weyand et al., 2009; Weyand et al., 2010).

A second, more complex variety of multi-mass models developed from the two-mass ideas initially put forward by McMahon (McMahon et al., 1987) and Alexander (Alexander, 1988). These models have evolved in their complexity, largely by building upon Alexander's two-mass, stacked-spring model (Alexander, 1988; Alexander, 1990; Derrick et al., 2000; Ker et al., 1989). Contemporary versions include at least four masses and more than a dozen spring, mass and damping elements (Liu and Nigg, 2000; Ly et al., 2010; Nigg and Liu, 1999; Nikooyan and Zadpoor, 2011; Zadpoor and Nikooyan, 2010). In contrast to the single-mass models, a primary objective of the multi-mass models has been to provide detailed explanations of waveform variability, specifically the impact and rising-edge variability observed for human joggers (Nigg, 2010; Zadpoor and Nikooyan, 2010). However, the relatively specific objective of the multi-mass models has limited the breadth of their application. Evaluations typically ignore the descending edge of the waveforms and have been limited to jogging speeds. Accordingly, the ability of the now-elaborate, multi-mass models to explain either the falling edge of jogging waveforms or the entirety of the waveforms from intermediate and fast running speeds is not known.

Here, we seek to explain running ground reaction forces in full with an approach that is slightly more complex than the single-mass models, but considerably simpler than current multi-mass models. For this purpose, we formulated a two-mass model that theorizes that running vertical force–time waveforms consist of two components, each corresponding to the motion of a discrete portion of the body's mass. A smaller component (*m*_{1}) corresponds to the impact of the lower limb with the running surface while a larger component (*m*_{2}) corresponds to the accelerations of the remainder of the body's mass (Fig. 1A). We hypothesize that our two-mass model can explain running ground reaction force–time waveforms in their entirety across different speed, foot-strike and footwear conditions.

## RESULTS AND DISCUSSION

In keeping with our hypothesis, our two-mass model was able to account for virtually all of the duration, amplitude and force–time pattern variability present in the vertical ground reaction force waveforms analyzed. Despite the large differences in waveform characteristics introduced by different speed, foot-strike and footwear conditions, our model accounted for an average of 97% of the individual force–time relationships (mean *R*^{2}=0.97±0.01) and a minimum of 95% (Fig. 1). The accuracy of these fits across the heterogeneous waveforms tested suggests that two mechanical phenomena, acting in parallel, are sufficient to explain running ground reaction forces: (1) the collision of the lower limb with the running surface, and (2) the motion of the remainder of the body's mass throughout the stance phase.

The accuracy of the fits achieved using a model with only two mass components and with mass component values held constant across conditions differs from prevailing paradigms in several respects. First, while the sequential additions of third, fourth and fifth mass components to multi-mass models over the last two decades (Liu and Nigg, 2000; Ly et al., 2010; Nigg and Liu, 1999; Nikooyan and Zadpoor, 2011; Zadpoor and Nikooyan, 2010) may describe physical and mechanical reality as theorized (Nigg, 2010; Zadpoor and Nikooyan, 2010), these additional masses may also be unnecessary for waveform prediction. Second, the conclusion that the mass quantity decelerated upon foot–ground impact differs substantially for rear foot versus forefoot impacts (Lieberman et al., 2010; Nigg, 2010) should be reconsidered. The close fits we report here using a constant value of 8.0% of the body's mass across all foot-strike conditions indicates that a variable ‘effective mass’ may be unnecessary for accurate modeling and could be mechanically incorrect. For example, if we predict the sprint running waveform analyzed here (Fig. 1G,H) using the effective mass proportions suggested for a forefoot impact [*m*_{1}=1.7% and *m*_{2}=98.3% of total body mass *M*_{b} (Lieberman et al., 2010)] with a speed-specific foot collisional velocity (Mann and Herman, 1985) (Table 1), the rising edge of the sprint waveform is substantially under-predicted and the overall goodness of fit is considerably reduced (*R*^{2}=0.95 to 0.82; see supplementary material Fig. S1). Third, the model's general features and simplifying assumptions permit impulse *J*_{1} and *J*_{2} durations and forces to be independent. In contrast, the dual ‘stacked spring-mass’ model-type (Alexander, 1990; Derrick et al., 2000; Ker et al., 1989) that Alexander originally introduced (Alexander, 1988) uses a serial, coupled configuration that may be incapable of predicting the brief simultaneous impulses responsible for the characteristic pattern of sprint running waveforms.

Indeed, the model's design was essential for achieving close fits to waveforms with variable rising edges, smooth falling edges and significantly different durations. Given the fixed-mass value of our lower-limb mass component, the close fits to the variable rising edges were achieved predominantly via the two model inputs (Table 1) responsible for the shape of collisional impulse *J*_{1} (Fig. 1). Values for the first of the two, the vertical velocity of the lower limb at touchdown, are well-supported by the waveform-specific literature values available. Values for the second, the deceleration time of the lower limb upon touchdown, are well-supported by the detailed analysis of Nigg et al. (Nigg et al., 1987) for waveform 2, but are lacking for the other three. Fits along the smoother falling edges depend directly upon impulse *J*_{2} because of the early conclusion of the *J*_{1} collisional event. Given a known physical basis for determining total impulse *J*_{T} from contact and step times (Eqn 1; see Materials and methods), correctly quantifying impulse *J*_{2} depends solely on the quantity subtracted for impulse *J*_{1} (Eqn 2). While empirical validation clearly remains for several elements of our model, the fits achieved using anatomical mass inputs, realistic lower-limb velocities, and one mechanical explanation across conditions raise the possibility that the running force–motion relationship may be more general than previously recognized.

An additional factor in the accuracy of the fits we report was undoubtedly the model evaluation method adopted. The method chosen allowed us to assess a greater variety of waveforms than would have been possible via direct experimentation, but also involved two potential limitations. First, because the model fits were generated by varying the inputs, the goodness-of-fit values obtained should be regarded as the upper performance limits of the model. Second, the digitizing process enabling our approach might have transformed the literature waveforms into more model-conducive shapes. We were able to evaluate this second possibility empirically by applying the inputs used to fit two of the digitized waveforms (3 and 4) to the original waveform data. This process yielded fits that were the same or slightly greater for the original (respective *R*^{2} values of 0.98 and 0.96) versus digitized versions because the original waveforms were so closely reproduced by digitizing (see supplementary material Tables S1–S4).

We close by providing respective, illustrative examples of the basic and applied advances made possible by the concise physical basis of our two-mass model. One basic insight provided by the framework of the model is the identification of a mechanical strategy that runners can adopt to achieve faster speeds. By simply increasing the lower limb's velocity prior to touchdown, and reducing deceleration time during impact, runners can elevate the collisional impulse (*J*_{1}) and total ground reaction forces as needed to attain faster speeds (Weyand et al., 2000; Weyand et al., 2009; Weyand et al., 2010). Both the existing literature data (Table 1) and our modeling results (Fig. 1) are consistent with this being a primary mechanism by which faster human runners do, in fact, attain faster sprint running speeds.

In application, the conciseness of the model could translate into practical techniques for determining ground reaction forces indirectly. At present, the lone indirect assessment method available (Bobbert et al., 1991) is scientifically rigorous, but impractical for broad usage. The existing technique involves the instantaneous summation of the accelerations of seven body segments based on high-frequency positional data from 10 bodily locations. In contrast, the scientific basis of our two-mass model (Eqns 1, 2, 3, 4, 5 and 6) reduces the data needed for indirect force determinations to three basic variables: aerial time, contact time and the vertical velocity of the lower limb. Thus, our model may allow video and other motion capture techniques to become practical tools for determining vertical ground reaction forces without direct measurement.

## MATERIALS AND METHODS

### Model formulation

*F*

_{Tavg}can be determined if foot–ground contact time

*t*

_{c}and aerial time

*t*

_{a}are known: where

*t*

_{step}is step time (

*t*

_{step}=

*t*

_{c}

*+t*

_{a}),

*m*is body mass and

**is gravitational acceleration (9.8 m s**

*g*^{−2}).

*J*

_{1}results from the acceleration of the lower limb during surface impact, and

*J*corresponds to the acceleration of the remainder of the body's mass. The total impulse

_{2}*J*

_{T}, is the sum of

*J*

_{1}and

*J*

_{2}:

*m*

_{1}is the 8.0% of the body's total mass attributed to the lower limb, while impulse mass

*m*

_{2}is the remaining 92.0%. Impulse

*J*

_{1}is quantified from the deceleration of

*m*

_{1}during surface impact: where Δ

*t*

_{1}is the time interval between touchdown and vertical velocity of

*m*

_{1}slowing to zero, Δ

*v*

_{1}is the change in vertical velocity of

*m*

_{1}during Δ

*t*

_{1}, and

*F*

_{1avg}is the average force during the total time interval (2Δ

*t*

_{1}) of impulse

*J*

_{1}. After the

*J*

_{1}time interval, the model assumes

*F*

_{1avg}=0.

*J*

_{2}is determined from

*J*

_{1}and total impulse

*J*

_{T}as: where

*F*

_{2avg}is the average force of

*J*

_{2}during the interval

*t*

_{c}.

### Modeled waveforms

*F*(

*t*) for

*J*

_{1}and

*J*

_{2}are a result of non-linear elastic collisions (Cross, 1999) that can be accurately modeled using the raised cosine function: where

*A*is the peak amplitude,

*B*is the center time of the peak and

*C*is the half-width time interval. Because of the symmetrical properties of this function, peak amplitude

*A*=2

*F*

_{avg}, and the area under the curve is

*J*=

*AC*. The total force waveform

*F*

_{T}(

*t*) is the sum of each impulse waveform:

*A*

_{1}is calculated from

*F*

_{1avg}using the Δ

*v*

_{1}and Δ

*t*

_{1}terms in Eqn 3, and

*B*

_{1}and

*C*

_{1}equal the time Δ

*t*

_{1}after touchdown for the vertical velocity of

*m*

_{1}to reach zero.

*A*

_{2}is calculated from

*F*

_{2avg}in Eqn 4, and

*B*

_{2}and

*C*

_{2}equal one-half the contact time

*t*

_{c}.

### Modeled versus actual waveforms

We digitized (Engauge, version 4.1) four published waveforms that varied in duration, amplitude and shape (Table 1). Model fits of the four digitized waveforms (Fig. 1) were performed via a manual iterative process that constrained the inputs for Δ*t*_{1} and Δ*v*_{1} to values deemed realistic on the basis of existing literature. Inputs for *t*_{c} and subsequent *t*_{a} were determined from the waveforms using a threshold of 60 N. In two cases (waveforms 3 and 4), goodness of fit between modeled and original data waveforms were determined to supplement the evaluation of the digitized versions.

Model fits were quantified in two ways: (1) in force units standardized to the body's weight (*W*_{b}) using the root mean square statistic (RMSE), and (2) for goodness of fit using the *R*^{2} statistic. Digitized waveforms were interpolated as needed to provide force data on a per millisecond basis for these analyses. We hypothesized that the model would explain 90% or more (i.e. *R*^{2}≥0.90) of the force–time variation present in each of the four waveforms analyzed. Data for all digitized, modeled and original waveforms used in the analysis are provided in supplementary material Tables S1–S4.

All variables are presented in SI units, but, per convention, force waveforms are illustrated in mass-specific units.

## Acknowledgements

The authors thank Lindsay Wohlers, Geoffrey Brown and the two anonymous reviewers for helpful comments on the manuscript.

## FOOTNOTES

**Funding**

This work was supported by a US Army Medical Research and Materiel Command award [W81XWH-12-2-0013] to P.G.W.

## References

**Competing interests**

The authors declare competing financial interests. Peter Weyand, Laurence Ryan and Kenneth Clark are the inventors of US Patent #8363891 which is owned by Southern Methodist University and contains scientific content related to that presented in the paper. The patent is licensed to SoleForce LLC in which the three aforementioned individuals are equity partners.