Abstract
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, footstrike 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.0M_{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, rearfoot strike run, forefoot strike run and forefoot 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 straightpath, level running and hopping is not in the horizontal direction of travel (Cavagna et al., 1977). Horizontal force requirements are minimized by an effective steptostep 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 strideaveraged 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 springmass model and multimass models. Models of both types are wellfounded 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 springmass model (Blickhan, 1989; McMahon and Cheng, 1990). The singlemass approach models running and hopping animals as a lumped pointmass mass bouncing on a massless leg spring. This singlemass 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 nonsymmetrical components that the force–time waveforms inevitably contain. These include, but are not limited to, heelstrike 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 multimass models developed from the twomass 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 twomass, stackedspring 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 singlemass models, a primary objective of the multimass models has been to provide detailed explanations of waveform variability, specifically the impact and risingedge variability observed for human joggers (Nigg, 2010; Zadpoor and Nikooyan, 2010). However, the relatively specific objective of the multimass 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 nowelaborate, multimass 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 singlemass models, but considerably simpler than current multimass models. For this purpose, we formulated a twomass 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 twomass model can explain running ground reaction force–time waveforms in their entirety across different speed, footstrike and footwear conditions.
RESULTS AND DISCUSSION
In keeping with our hypothesis, our twomass 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, footstrike 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 multimass 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 footstrike 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 speedspecific foot collisional velocity (Mann and Herman, 1985) (Table 1), the rising edge of the sprint waveform is substantially underpredicted 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 springmass’ modeltype (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 fixedmass value of our lowerlimb 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 wellsupported by the waveformspecific literature values available. Values for the second, the deceleration time of the lower limb upon touchdown, are wellsupported 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 lowerlimb 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 goodnessoffit 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 modelconducive 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 twomass 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 highfrequency positional data from 10 bodily locations. In contrast, the scientific basis of our twomass 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
Because the net vertical displacement of the body over time during steadyspeed, level running is zero, the timeaveraged vertical ground reaction force must equal the body's weight. Thus, the total stanceaveraged vertical force F_{Tavg} can be determined if foot–ground contact time t_{c} and aerial time t_{a} are known: (1) where t_{step} is step time (t_{step}=t_{c}+t_{a}), m is body mass and g is gravitational acceleration (9.8 m s^{−2}).
The ground reaction force waveform represents the instantaneous acceleration of the body's mass. Accordingly, the waveform can be conceptualized as the sum of the instantaneous accelerations of different segments that make up the body's total mass (Bobbert et al., 1991). In our model (Fig. 1), impulse J_{1} results from the acceleration of the lower limb during surface impact, and J_{2} corresponds to the acceleration of the remainder of the body's mass. The total impulse J_{T}, is the sum of J_{1} and J_{2}: (2)
Impulse mass 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: (3) 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: (4) where F_{2avg} is the average force of J_{2} during the interval t_{c}.
Modeled waveforms
The bellshaped force curves F(t) for J_{1} and J_{2} are a result of nonlinear elastic collisions (Cross, 1999) that can be accurately modeled using the raised cosine function: (5) where A is the peak amplitude, B is the center time of the peak and C is the halfwidth time interval. Because of the symmetrical properties of this function, peak amplitude A=2F_{avg}, and the area under the curve is J=AC. The total force waveform F_{T}(t) is the sum of each impulse waveform: (6) 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 onehalf 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 massspecific units.
ACKNOWLEDGEMENTS
The authors thank Lindsay Wohlers, Geoffrey Brown and the two anonymous reviewers for helpful comments on the manuscript.
FOOTNOTES

Author contributions
Each of the three authors, K.P.C., L.J.R. and P.G.W., contributed substantially to the conception of the study, the implementation and evaluation of the model presented, and the writing of the manuscript.

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.

Funding
This work was supported by a US Army Medical Research and Materiel Command award [W81XWH1220013] to P.G.W.

Supplementary material
Supplementary material available online at http://jeb.biologists.org/lookup/suppl/doi:10.1242/jeb.099523//DC1
 © 2014. Published by The Company of Biologists Ltd