## ABSTRACT

The relationship between gait mechanics and running ground reaction forces is widely regarded as complex. This viewpoint has evolved primarily via efforts to explain the rising edge of vertical force–time waveforms observed during slow human running. Existing theoretical models do provide good rising-edge fits, but require more than a dozen input variables to sum the force contributions of four or more vague components of the body's total mass (*m*_{b}). Here, we hypothesized that the force contributions of two discrete body mass components are sufficient to account for vertical ground reaction force–time waveform patterns in full (stance foot and shank, *m*_{1}=0.08*m*_{b}; remaining mass, *m*_{2}=0.92*m*_{b}). We tested this hypothesis directly by acquiring simultaneous limb motion and ground reaction force data across a broad range of running speeds (3.0–11.1 m s^{−1}) from 42 subjects who differed in body mass (range: 43–105 kg) and foot-strike mechanics. Predicted waveforms were generated from our two-mass model using body mass and three stride-specific measures: contact time, aerial time and lower limb vertical acceleration during impact. Measured waveforms (*N*=500) differed in shape and varied by more than twofold in amplitude and duration. Nonetheless, the overall agreement between the 500 measured waveforms and those generated independently by the model approached unity (*R*^{2}=0.95±0.04, mean±s.d.), with minimal variation across the slow, medium and fast running speeds tested (Δ*R*^{2}≤0.04), and between rear-foot (*R*^{2}=0.94±0.04, *N*=177) versus fore-foot (*R*^{2}=0.95±0.04, *N*=323) strike mechanics. We conclude that the motion of two anatomically discrete components of the body's mass is sufficient to explain the vertical ground reaction force–time waveform patterns observed during human running.

## INTRODUCTION

Running ground reaction forces are of fundamental physical and biological importance. Acutely, they determine the body's state of motion, limb-loading rates and tissue stresses. Over time, they influence the health and functional status of the tissues, limb and runner. However, the mechanical basis of the vertical force–time waveform patterns described most extensively for human runners (Cavanagh, 1987; Munro et al., 1987) continues to be a matter of significant disagreement (Chi and Schmitt, 2005; Clark et al., 2014; Denoth, 1986; Derrick, 2004; Lieberman et al., 2010; Nigg, 2010; Shorten and Mientjes, 2011). The current discordance is at least partially attributable to the mechanical complexities of the limbs and bodies responsible for the forces present at the limb–ground interface. Human and other vertebrate runners are mechanically complex in their body and limb-segment morphology, tissue properties, neural control of muscle forces, and joint and limb stiffnesses. These features allow body mass components to accelerate differently with respect to one another and the ground. Because the total waveform corresponds to the acceleration of the body's entire mass, the summed accelerations of different mass components must somehow determine the instantaneous forces and waveform patterns observed.

At present, the most common approach to explaining the force–time waveform patterns of human runners are lumped element, spring–mass systems that include four or more hypothetical mass components with multiple springs and dashpots (Fig. 1A). Most current versions include 14 or more input variables derived via forward dynamics simulations (Chi and Schmitt, 2005; Liu and Nigg, 2000; Ly et al., 2010; Nigg and Liu, 1999; Nikooyan and Zadpoor, 2011; Zadpoor and Nikooyan, 2010). Per their intended purpose, these models are able to provide close, *post facto* fits to the rising edge of the force–time waveforms that result from rear-foot strike mechanics at jogging speeds under a variety of surface, footwear and other conditions (Ly et al., 2010; Zadpoor and Nikooyan, 2010). However, these models do not attempt to predict the falling edge of the waveform, they do not explain the differently shaped waveforms that typically result from fore-foot strike mechanics, and their ability to fit waveforms from intermediate and fast running speeds is completely unknown.

A scientific justification for including numerous mass components to account for vertical ground reaction force–time waveforms was importantly provided by a direct motion-to-force experiment conducted by Bobbert et al. (1991) a quarter century ago. These investigators demonstrated that the total waveform can indeed be reasonably predicted from the summed accelerations of seven body mass components at modest running speeds. Recently, we have theorized that an alternative approach may allow the number of masses needed for full waveform prediction to be reduced from seven to only two.

Our approach (Clark et al., 2014) divides the body's mass into two, anatomically based, invariant, mass components: the first corresponding to the mass of the lower limb (*m*_{1}, 8% total body mass) and the second corresponding to the remainder of the body's mass (*m*_{2}, 92% total body mass; Fig. 1B). The model theoretically allows the full vertical force–time waveform to be predicted from the force contributions corresponding to the two body mass components. Impulse 1 results from the vertical collision of *m*_{1} with the running surface and impulse 2 results from the vertical accelerations of *m*_{2} throughout the ground contact period (Clark et al., 2014). Our introductory effort (Clark et al., 2014) was able to account for essentially all of the variation present in four vertical force–time waveforms (*R*^{2} range: 0.95–0.98, mean=0.97±0.01) selected for their amplitude, duration and shape heterogeneity. However, the close fits achieved resulted from *post facto* selection of the input parameters to maximize the goodness of each fit.

Here, we undertook a direct, experimental test of the hypothesis that the motion of two discrete body mass components is sufficient to predict running vertical ground reaction force–time waveforms in full. The two-mass model requires only body mass and three stride-specific measures as inputs: contact time, aerial time and lower limb acceleration. Because these three inputs can be readily acquired from a variety of video or other inexpensive motion-sensing technologies, our model potentially offers economical options for generating running ground reaction force waveforms without a force platform. Additionally, the concise running force–motion linkage provided could be applied to footwear, prosthetic and orthotic design, or used to inform gait interventions designed to reduce injuries or enhance running performance.

## MATERIALS AND METHODS

### Model formulation

Although the physical motion of running occurs in three spatial dimensions, the total vertical ground reaction force waveform is determined by the forces due to the instantaneous vertical accelerations of the body mass components. Our computational model utilizes experimental measurements to determine the parameters linking running motion to impulse forces, thus avoiding the limitations of modeling ground reaction forces with single-axis or single-body mass–spring-damper systems (Nikooyan and Zadpoor, 2011). The fundamental premise of the two-mass model (see Appendix for detailed derivation) is that the total vertical ground reaction force waveform is composed of two overlapping bell-shaped impulses due to the vertical collision of the lower limb (*J*_{1}) with the running surface and the concurrent vertical accelerations of the rest of the body (*J*_{2}) during ground contact. The total ground reaction impulse *J*_{T} is the sum of *J*_{1} and *J*_{2} and can be determined by the total stance-averaged vertical ground reaction force *F*_{T,avg} during the ground contact time *t*_{c}:
(1)

The model assumes steady-speed level running where the speed is constant and the net vertical displacement of the center of mass of the body is zero over each step. Thus, the time-averaged vertical ground reaction force must equal the body's weight and *F*_{T,avg} can be determined if contact time *t*_{c} and aerial time *t*_{a} are known:
(2)

where *m*_{b} is body mass and ** g** is gravitational acceleration (

**=9.8 m s**

*g*^{−2}), and the quantity

*t*

_{c}+

*t*

_{a}equals the step time,

*t*

_{step}.

Impulse *J*_{1} corresponds to the vertical deceleration of *m*_{1} during surface impact:
(3)

where *m*_{1} is 8.0% of the body's mass (Plagenhoef et al., 1983; Winter, 1990), Δ*t*_{1} is the time interval between touchdown and the vertical velocity of *m*_{1} slowing to zero, Δ*v*_{1} is the change in vertical velocity of *m*_{1} during Δ*t*_{1}, and *F*_{1,avg} is the average force during the total time interval (2Δ*t*_{1}) of impulse *J*_{1}, here defined as the impact interval. A single ankle marker is used to measure Δ*v*_{1} and Δ*t*_{1}, which determine the vertical acceleration of the lower limb mass *m*_{1} (Fig. 2). The lower limb attains a relatively constant positive velocity after the impact interval (Fig. 2C), resulting in a near-zero acceleration of *m*_{1} (Fig. 2D) and negligible force (Fig. 2E). Thus, *J*_{1} can be represented by a finite impulse during the impact interval. Kinematic data for the ankle marker position during the impact interval for representative rear-foot strike (RFS) and fore-foot strike (FFS) subjects appear in Fig. 3. For highest accuracy of the Δ*t*_{1} measurement, a velocity threshold and projection method is used to eliminate minor fluctuations in the ankle marker velocity profile near zero (see Appendix). To assess the accuracy of this projection method, the measured ankle marker velocity at the projected time to zero was quantified for all footfalls.

Impulse *J*_{2} corresponds to the remainder of the body's mass and is determined from total ground reaction impulse *J*_{T} in Eqn 1 and impulse *J*_{1} in Eqn 3:
(4)

where *F*_{2,avg} is the average force of impulse *J*_{2} during the contact time *t*_{c}.

The raised cosine bell (RCB) curve function was used to generate the *F*(*t*) waveforms of *J*_{1} and *J*_{2} for both foundational and empirical reasons. The RCB function is unique among all bell curve functions in that it can be derived from the first two terms of the Fourier series (see Appendix). Analyses of vertical ground reaction force–time waveforms from jumping in place indicate that using a function consisting of two bell curves provides a superior representation and requires a lower number of Fourier terms than a half-sine curve (Racic and Pavic, 2009). Direct waveform comparisons indicate that the RCB function provides superior descriptions versus both the half-cosine and Gaussian functions (see cos^{2} data in table 2 of Sim et al., 2008) when vertical ground reaction force impulses are generated at frequencies ≥2 Hz.

Thus, each impulse uses the RCB function for the force waveform *F*(*t*) during the interval *B−C*≤*t*≤B+*C*:
(5)

where *A*=2*F*_{avg} is the peak amplitude, *B* is the center time of the peak and *C* is the half-width time interval. For higher accuracy, the waveform function for impulse *J*_{2} includes an offset to account for the force plate threshold and an asymmetry adjustment to account for the center of mass height difference at takeoff and touchdown. A *J*_{2} force peak corresponding to the minimum height of the center of mass was set at 0.47*t*_{c} based on the observations of Cavagna et al. (1977, see their table 4) for human running (see Appendix for additional details).

The total force curve *F*_{T}(*t*) is the sum of the two individual force waveforms representing impulse *J*_{1} and impulse *J*_{2} :
(6)

### Model force–time waveforms

Model-predicted versus measured vertical force waveforms were tested for both individual footfalls and trial-averaged force data. For the individual footfall predictions, the input parameters of body mass, ground contact time, subsequent aerial time, *m*_{1} vertical velocity at touchdown and impact interval time Δ*t*_{1} were used to generate the model-predicted waveform, which was then compared with the measured waveform from that footfall. For the trial-averaged waveform comparisons, the input parameters from each right footfall were averaged for that trial. A minimum of three and a maximum of six right footfalls were included in the trial average, depending on the number of steps completed during the trial. The measured input parameters from each trial average were used to generate a model-predicted waveform, which was then compared with the measured trial-averaged waveform. An example of this data treatment appears in Fig. 4, including a series of original measured waveforms (Fig. 4A), the trial-average model-predicted waveform (Fig. 4B), and the model-predicted versus trial-average measured waveform (Fig. 4C). Both single-footfall and trial-averaged predictions were assessed to evaluate whether the number of footfalls included influences the predictive accuracy of the model.

To validate the anatomical mass fractions in our model, waveforms were alternatively generated with literature-suggested impact-mass minimum and maximum values (Derrick, 2004). A smaller *m*_{1} quantity of 1.5% was used (representing the approximate mass of the foot; Plagenhoef et al., 1983; Winter, 1990; Hamill and Knutzen, 2009) with a corresponding *m*_{2} quantity of 98.5%. A larger *m*_{1} quantity of 16% was used (representing the approximate mass of the entire stance limb; Plagenhoef et al., 1983; Winter, 1990; Hamill and Knutzen, 2009) with a corresponding *m*_{2} quantity of 84.0%. We predicted that across both speed and foot-strike mechanics, values along the rising edge of the waveform would be under-predicted with an *m*_{1} quantity of 1.5% and over-predicted with an *m*_{1} quantity of 16%.

### Statistical analysis

The predictive accuracy of the model was assessed on the 500 individual footfalls acquired using both the *R*^{2} statistic and the root mean square error (RMSE) statistic, global values that quantify goodness of fit in relative and absolute terms, respectively. These statistical assessments are broadly used for a variety of purposes, including quantifying the degree of overlap present in time-series data per prior practice (Cohen, 2013; Clark et al., 2014; Clark and Weyand, 2014; Morin et al., 2005). The footfall sample size was sufficiently large to detect very small differences in model performance across foot-strike type and speed using the *R*^{2} and RMSE statistics.

A two-way ANOVA was used to evaluate whether goodness-of-fit and RMSE values varied as a function of speed (slow, medium, fast) and foot-strike (rear- versus fore-foot categories) classifications (*P*≤0.05, Table 1) on all 500 footfall waveforms acquired.

### Subjects and participation

A total of 42 subjects, 23 men and 19 women, volunteered and provided written informed consent. The consent process and all experimental procedures were approved by, and conformed to, the approval terms granted by the Institutional Review Board of Southern Methodist University. All subjects were between 18 and 37 years of age and engaged in regular physical activity at the time of testing. The mean age and size characteristics of the men and women were as follows: men: age=23.3±5.0 years, range=18–37 years; height=1.79±0.07 m, range=1.69–1.95 m; mass=81.1±8.5 kg, range=71.0–101.5 kg; and women: age=22.5±1.7 years, range=18–36 years; height=1.68±0.06 m, range=1.55–1.78 m; mass=63.3±9.4 kg, range=43.4–82.0 kg. Subjects included recreationally trained individuals, intercollegiate team-sport athletes, and professional track and field athletes, four of whom were Olympic medalists in sprint or hurdle events.

### Data acquisition

Data were collected across a range of speeds (3.0 m s^{−1} to top speed) on a three-axis, custom-built high-speed force treadmill (AMTI, Watertown, MA, USA) capable of speeds of over 20 m s^{−1}. To ensure that the model was being evaluated at speeds that required a normal running gait, only trials above a Froude speed of 1.0 [*v*/√(*g**L*_{0})>1.0] were included in the statistical analysis. For both submaximal and maximal tests, subjects followed testing procedures similar to those described in Weyand et al. (2000, 2010). Thirty-nine of the 42 subjects completed trials up to top speed, and these subjects were habituated to high-speed treadmill running during one or more familiarization sessions before undergoing top speed testing. All subjects wore standardized black compression shirts and shorts, and the same model of running shoes (Nike Waffle Racer version 7.0, Beaverton, OR, USA).

### Kinetic and kinematic data collection and analysis

Ground reaction force data were acquired at 1000 Hz and were post-filtered using a low-pass, fourth-order, zero-phase-shift Butterworth filter with a cutoff frequency of 25 Hz (Winter, 1990). For each footfall, the vertical ground reaction force applied during the contact period was determined from the time during which the vertical force signal exceeded a threshold of 40 N. Additionally, trial-averaged vertical ground reaction force waveforms were generated for individual subjects at different trial speeds by averaging the force from each millisecond of the contact period for the right-foot waveforms that had corresponding kinematic data.

For each trial, 3.46 s of video data were collected using a high-speed video system consisting of three Fastec Imaging HiSpec 2G cameras (Mikrotron GmbH, Unterschleissheim, Germany) operating at 1000 frames s^{−1}. Subjects wore reflective markers on the heel and ball of the foot on the lateral aspect of the right running shoe, as well as on the lateral aspect of the joint axes of rotation of the right ankle, knee and hip to capture these respective positions during the trials. A single-frame video file of each subject was recorded prior to testing to serve as a reference for the kinematic analyses. To assess the predictive accuracy of the model across different foot-strike types, footfalls were classified as RFS if the first surface contact occurred on the posterior half of the foot, and FFS if the first surface contact occurred on the anterior half.

The procedures used in extracting three-dimensional marker coordinates from the high-quality multiple camera videos consisted of correcting image distortions, calibrating the three-dimensional space and digitizing the markers. The calibration and digitization routines extensively used functions from the MATLAB Image Processing Toolbox (MathWorks, Natick, MA, USA). The calibration and digitization MATLAB routines were developed by the Hedrick Lab at the University of North Carolina (Hedrick, 2008). A resolution of 0.7 mm was measured under dynamic conditions using the high-speed video system with the custom MATLAB image correction, calibration and digitization routines. Data acquisition timing for the AMTI DigiAmp amplifier and Fastec Imaging cameras was synchronized through hardware interfaces. The digitized marker data were filtered at 25 Hz using the same filter as for the force data.

## RESULTS

### Overall agreement between model-predicted and measured waveforms

The goodness-of-fit agreement (*R*^{2}) between the 500 ground reaction force waveforms we measured and those predicted by our two-mass model approached unity as hypothesized (Table 1). The corresponding error of prediction expressed in force units standardized to the body's weight was slightly greater than 0.2*W*_{b} and was equal to 11.5% of the mean stance-averaged vertical force from all 500 trials (1.82±0.23*W*_{b}). The *R*^{2} and RMSE statistics for the 108 trial-averaged waveform values were nearly identical to those for the 500 individual waveforms. Specifically, the goodness of fit between model-predicted and measured trial-averaged ensemble waveforms was only 0.01 greater (*R*^{2}=0.96±0.03) than the individual footfall value, while the RMSE of prediction was only 0.02*W*_{b} less (0.19±0.07*W*_{b}) than the corresponding individual footfall values.

### Predictive accuracy across foot-strike types and running speeds

The goodness-of-fit between the model-predicted and measured waveforms (*R*^{2}) was nearly identical (Δ*R*^{2}=0.01) for the 177 RFS versus 323 FFS waveforms (Table 1). Because of the large number of footfalls tested and minimal variability in model predictive accuracy, the 0.01 greater *R*^{2} value for the FFS versus RFS waveforms was significantly different statistically. Corresponding RMSE differences for the RFS versus FFS waveform predictions were small (ΔRMSE=+0.04*W*_{b}) and not statistically different.

The goodness-of-fit between the model-generated and actual waveforms was slightly, but significantly, better for waveforms acquired from slow and medium speeds versus those from fast speeds. The *R*^{2} values, when averaged for the waveforms from the three speed ranges (*R*^{2} total, Table 1) regardless of foot-strike type, were slightly, but significantly lower (Δ*R*^{2}=−0.03) for the faster versus medium and slower speed trials. Similarly, the RMSE of the model fits to the slow and medium speed waveforms were not different from one another; however, both were significantly smaller than the error of prediction for the waveforms from the fastest speeds. Similar across-speed patterns were present for both the *R*^{2} and RMSE statistics when the waveforms were analyzed within either the RFS or FFS mechanics categories (Table 1).

The waveforms generated with the two-mass model accurately predicted the more rapid rising edges of the RFS versus FFS waveforms, including transient rising-edge peaks when present, regardless of the running speed, and total waveform amplitude and duration (Figs 5A,C,E and 6A,C,E). Waveform shape variability across foot-strike types was accurately predicted from the shorter deceleration periods of the *m*_{1} mass segment for RFS versus FFS (Δ*t*_{1}, Table 2) with little difference in the overall *J*_{1} impulse values (Fig. 5B,D,F versus Fig. 6B,D,F; Δ*v*_{1}, Table 2) at similar speeds. From slower to faster speeds, *J*_{1} impulses predicted by the model became greater for both RFS and FFS waveforms as a result of the greater pre-impact ankle velocities (Table 2). The close fits to the middle and later portions of all the waveforms resulted largely or exclusively from the *J*_{2} impulse predicted from the model because the *m*_{1} impact deceleration event concluded during the first quarter to half of the total contact period.

The measured ankle marker velocity at the time our projection technique identified a zero value was, on average, −0.05±0.04 m s^{−1} for the 500 individual trials. All but two of the footfalls had a measured velocity within ±0.20 m s^{−1} of a literal zero value (i.e. 0.00 m s^{−1}).

### Model predictive accuracy with alternative *m*_{1} segment values

Poorer agreement between model-predicted and measured waveforms resulted from using *m*_{1} segment values that were either smaller or larger than the originally assigned anatomical model value of 8.0% of *m*_{b}. As hypothesized, waveforms generated using lesser *m*_{1} values (*m*_{1}=1.5% with *m*_{2}=98.5% of *m*_{b}) consistently under-predicted the force values measured along the rising edge of the waveforms (Figs 7A,D,G and 8A,D,G). Conversely, waveforms generated using greater *m*_{1} values (*m*_{1}=16.0% with *m*_{2}=84.0% of *m*_{b}) consistently over-predicted measured rising-edge force values (Figs 7C,F,I and 8C,F,I). Differences were more pronounced at the faster trial speeds because of the greater *m*_{1} segment decelerations and correspondingly larger *J*_{1} impulses at faster speeds.

The *R*^{2} goodness of fit values for the 500 individual footfall waveforms generated using *m*_{1} values equal to 1.5% and 16.0% of *m*_{b} accounted for 12% and 21% less variance, respectively, versus the original *m*_{1}=8.0% results (*m*_{1}−1.5% *R*^{2}=0.83±0.16 and *m*_{1}−16.0% *R*^{2}=0.74±0.21). Predictive error values for the waveforms generated using the two alternative *m*_{1} segment values were approximately twice as large as those obtained using the original 8.0% value (*m*_{1}−1.5% RMSE=0.36±0.17*W*_{b}; *m*_{1}−16.0% RMSE=0.47±0.24*W*_{b}).

## DISCUSSION

As hypothesized, we found that the force contributions of two discrete body mass components do indeed suffice to predict running vertical ground reaction force–time waveforms in full. Our two-mass, two-impulse model independently predicted nearly all of the variation in measured running ground reaction force waveforms, which differed considerably in their shape, amplitude and duration. Although our prior, best-fit analysis indicated that this outcome might be theoretically possible (Clark et al., 2014), an experimental test incorporating simultaneous motion and force data had not been previously undertaken. Our direct test here indicated that regardless of whether our 42 human subjects jogged, ran at intermediate speeds or sprinted, and whether they first contacted the running surface with the fore or rear portions of their feet, there was near-complete agreement between the model-predicted and measured force–time waveforms across the 500 footfalls we analyzed (Table 1). Thus, we conclude that the vertical ground reaction forces of human runners can be broadly understood from the motion of two discrete portions of the body: (1) the contacting lower limb and (2) the remainder of the body's mass.

### Two-mass model: scientific and technical elements

A primary scientific challenge here was not knowing *a priori* how well the waveform contributions from the 92% of the body's mass could be predicted when modeled as a single mass component. Conceivably, the summed force contributions resulting from the motion of the head, arms, trunk, upper portion of the contacting leg and full mass of the non-contacting leg might defy accurate prediction when treated as a single mass (Bobbert et al., 1991). This large, multi-jointed mass component lacks a fixed, readily measurable center because of the non-uniform motion of the different segments that comprise it. Hence, the complexity of the within-segment and total motion of our model's mass component *m*_{2} and its corresponding force contribution would be difficult to measure and predict from positional data. We ultimately implemented an indirect approach that allowed the force contributions of mass *m*_{2} to be quantified without position and time data. We simply subtracted impulse *J*_{1} from the total ground reaction impulse *J*_{T}, after quantifying the latter from body mass, gravity and the contact proportion of the total step time (Eqns 1–3). The resulting fits supported the efficacy of the approach as the agreement between model-predicted and measured waveforms was consistently excellent over the later portions of the stance period where the total waveform is predominantly determined by impulse *J*_{2} (Figs 5A,C,E and 6A,C,E, waveform trailing edges).

A potential limitation of our impulse subtraction approach was the required assumption that the net vertical displacement of the body's center of mass over the course of one or many strides is zero. However, our analysis indicated that little to no predictive error was introduced by this assumption under our level, treadmill-running test conditions. We found similar levels of predictive accuracy for trial-averaged and individual-footfall waveforms even though non-zero vertical displacements of the center of mass, when considered on a per-step basis, could have been substantially greater over the course of an individual step versus multiple consecutive steps. However, for our across-speed comparisons, it seems likely that our slightly less accurate waveform predictions for the fastest speeds, where step-to-step mechanics are generally less consistent versus medium and slower speeds (Weyand et al., 2010), are probably attributable to marginally greater violations of this assumption for the faster trials (Table 1).

An additional technical challenge involved accurate and consistent quantification of the duration of mass component *m*_{1}'s impact-period deceleration to a zero velocity (i.e. Δ*t*_{1} in Eqn 3; see also Fig. 2) for all footfalls. Following RFS impacts, the ankle marker position–time trajectories consistently provided a discrete vertical position minimum corresponding to the zero velocity needed to quantify half-impact duration Δ*t*_{1} in our model (Fig. 3). However, for some FFS footfalls, the rates of the positional change versus time during the last fraction of the deceleration period were less consistent and more prolonged than the FFS data in Fig. 3. This led us to implement an ankle marker projection technique (see Appendix) to avoid basing impulse *J*_{1} predictions on Δ*t*_{1} values that, in these cases, are not representative of the overall timing of the impact-deceleration event. Upon implementation across all 500 footfalls, we found that the measured ankle marker velocities at the time that the projection technique identified a zero velocity value were very near the true zero value desired (mean *v*_{ankle}=−0.05±0.04 m s^{−1}). Given the overall mean ankle marker vertical velocity, Δ*v*_{1}, of −1.59±0.02 m s^{−1} at impact across all 500 footfalls (Δ*v*_{1}, Table 2), this method, on average, captured 97% of *m*_{1} total post-impact deceleration to zero.

The ability to consistently predict rising-edge waveform peaks that occurred from 15 to 50 ms after initial impact is a noteworthy aspect of our model validation. As the timing of rising-edge peaks on individual waveforms is determined by the overlap of the two impulses in our model, successful prediction required precisely capturing the timing of both the high- (*J*_{1}) and low-frequency (*J*_{2}) components of the waveforms. The timing of impulse *J*_{1} was determined from motion data, and was therefore fully independent of our force data, filtering and processing routines. Conversely, the timing of impulse *J*_{2} was directly dependent on our force data and processing routines, and was therefore fully independent of our kinematic data and processing routines. Had even a minor degree of temporal inaccuracy been present in either our kinematic or kinetic data, the predictive accuracy with which the two-mass model identified rising-edge force peaks would not have been possible.

### Integrating two-mass model and multi-mass model results

Our experimental goal was to identify the most concise mechanical explanation that running ground reaction force waveforms might have. The multi-mass models, in contrast, were formulated to evaluate the potential influence of numerous factors on the waveform rising edge. Many of the features incorporated into the multi-mass models provide reasonable theoretical representations of the numerous, potentially influential musculoskeletal complexities present (Liu and Nigg, 2000; Ly et al., 2010; Nigg and Liu, 1999; Nikooyan and Zadpoor, 2011; Zadpoor and Nikooyan, 2010). These include mass components that vary in stiffness, that are both rigid and wobbling in nature, and that are connected with both serial and parallel elements (Fig. 1A). While the quantitative descriptions derived for these features undoubtedly include uncertainty, their basic influence on the RFS jogging waveforms thus far analyzed are plausible and largely consistent with other experimental and modeling approaches (Gruber et al., 1998; LaFortune et al., 1996; Pain and Challis, 2001; Shorten and Mientjes, 2011). Accordingly, our demonstration that a substantially more concise model can account for running ground reaction force waveforms in full should not be regarded as incompatible with the more theoretical results the multi-mass models have provided. Indeed, the most reasonable conclusion from the different approaches is that the collective influence of the many biological complexities incorporated into the multi-mass models is, in sum, accurately described by the concise impulse–momentum relationships the two-mass model provides.

### The force–motion relationship for human running: general or foot-strike specific?

In contrast to the prevailing view that the impact forces resulting from RFSs and FFSs involve different mass quantities (Lieberman et al., 2010; Nigg, 2010), our results indicate that mass quantities and force–motion relationships do not differ across strike types. This becomes fully apparent when *J*_{1} impact impulses are quantified using measured, foot-strike-specific deceleration durations in conjunction with the invariant, anatomically based mass quantities in our model. With both factors in place, we were able to accurately predict the waveforms in full for both foot-strike types at slow, intermediate and fast running speeds alike.

As is evident from both our impulse *J*_{1} illustrations (Figs 5 and 6) and the Δ*v*_{1} values in Table 2, rising-edge force peaks that are generally visible for RFS but not FFS waveforms can be fully attributed to the different *m*_{1} deceleration durations we measured. Longer FFS deceleration durations reduce and delay the *J*_{1} peak force values; they also result in greater waveform force contributions from impulse *J*_{2} at the time of the *J*_{1} force peak. In combination, these timing-dependent factors typically do not allow the *J*_{1} impulse peak to introduce a localized force peak along the rising edge of the full waveform (Fig. 6B and Fig. 8B,E). We found this to be the case even though the total ground reaction impulses resulting from FFSs in our data set were just as large as those resulting from RFSs (note: *J*_{1} impulses are ∝Δ*v*_{1} in Table 2). One noteworthy exception to these general foot-strike-specific waveform patterns has recently been documented for specialized sprinters (Clark and Weyand, 2014). These fore-foot striking athletes have waveforms that are characterized by prominent rising-edge force peaks that result from brief, large *J*_{1} impulse peaks. The high pre-impact limb velocities (Δ*v*_{1}) and brief impact deceleration periods (Δ*t*_{1}) of specialized sprinters introduce conspicuous rising-edge force peaks, particularly at faster running speeds (Fig. 6E,F and Fig. 8H).

Our waveform predictions using alternative mass fractions for segments *m*_{1} and *m*_{2} also support the validity of describing FFS and RFS waveforms with the same fractional body mass quantities and force–motion relationships. As hypothesized for these alternative mass distributions, substantially reducing the lower limb mass value of *m*_{1} from the anatomical fractional value of 0.08 (Plagenhoef et al., 1983) to the much lower value of 0.015 resulted in predicted rising-edge force values that fell consistently below measured values (Figs 7A,D,G and 8A,D,G). Conversely, substantially increasing the *m*_{1} fractional mass value to 0.16 resulted in predicted rising-edge force values that fell consistently above the measured values (Figs 7C,F,I and 8C,F,I). The magnitudes of the respective predictive errors and trends observed across speed were similar for the two foot-strike types using the aforementioned alternative *m*_{1} fractional mass values. In both cases, the increases in the magnitude of the *J*_{1} impulse across speed introduced greater predictive errors at higher versus lower speeds. More globally, the contrast between the systematic predictive errors present in both sets of alternative-mass generated waveforms versus the near-full agreement achieved with the original values provides strong support for the validity of the anatomical values originally assigned to mass components *m*_{1} and *m*_{2} in our two-mass model.

### General considerations, applications and concluding remarks

The accuracy and conciseness of our mechanical explanation for the variable ground reaction force–time patterns of dozens of human runners performing across their full range of speeds raises a noteworthy possibility. Our two-mass, two-impulse approach may offer a mechanical explanation that generalizes to the ground reaction force–time patterns of other running species. However, for non-human runners, the lesser distal-limb mass segments typically present (Hildebrand, 1960; Rubenson et al., 2011) could alter both the form and applicability of the two-mass approach. Minimally, species-specific mass distributions would need to be incorporated to generate waveforms from the basic stride measures included in the model. In this regard, comparative anatomical variation constitutes both a challenge and an experimental opportunity to evaluate the basic tenets of the model. More broadly, the testable hypotheses the model offers should be tractable using a variety of direct and indirect approaches.

Finally, our empirical validation of a concise model that can fully predict running vertical ground reaction forces offers basic insights with immediate potential for application. In contrast to simplified, single-mass models (Blickhan, 1989; Blum et al., 2009; Farley and Gonzalez, 1996; McMahon and Cheng, 1990; Silder et al., 2015), which are incapable of capturing the high-frequency, impact-phase components of the waveform (Bullimore and Burn, 2007; Clark and Weyand, 2014; Shorten and Mientjes, 2011), and multi-mass models that do not account for the whole waveform and are too theoretical and complex for practical application, the two-mass model requires only body mass and very limited motion data in order to predict the waveform in full. These attributes enable indirect assessment of impact forces, limb loading rates, and other informative, force-based outcomes using video or other inexpensive motion-sensing technologies. Potential model applications include: informing the design of running robots, exoskeletons and prostheses, the customization of running shoes and orthotics from individual gait mechanics, the development of wearable ground force sensors, and the improvement of evidenced-based feedback for gait analyses, intervention and modification.

## Acknowledgements

The authors thank Andrew Udofa for his help with data collection and the artwork in Fig. 1B, and the participants for their time and effort.

## APPENDIX

### Impulse determination

The model assumes steady-speed horizontal running where the net vertical displacement of the center of mass of the body is zero over each step and the speed is constant. Each step is defined by a contact time *t*_{c} with vertical ground reaction force *F*(*t*), and an aerial time *t*_{a} where the force is zero. Under these conditions, a runner supports an average of one body weight during the step time (*t*_{step}*=t*_{c}+*t*_{a}). This can be expressed using the formal mathematical definition of the average value of the force function *F*(*t*) over the interval *t=*0 to *t*_{step}:
(A 1)

where body weight is defined by the product of mass *m*_{b} and gravitational acceleration ** g**=9.8 m s

^{−2}. This equation yields the total average force

*F*

_{T,avg}during the contact time interval: (A 2)

The total ground reaction impulse *J*_{T} can simply be determined from body weight and *t*_{step}:
(A 3)

The ground reaction force is a result of the acceleration *a** _{i}*(

*t*) of body components

*i*with mass

*m*contacting the ground: (A 4)

_{i}The average force due to each body component is: (A 5) (A 6)

The ground reaction force is the sum of two distinct impulse waveforms. Each impulse waveform has an associated *F _{i}*

_{,avg}and time interval. Impulse

*J*

_{1}results from the acceleration of the lower limb during the limb impact interval. Impulse

*J*

_{2}results from the acceleration of the remainder of the body's mass during the entire contact interval. The total ground reaction impulse

*J*

_{T}is the sum of

*J*

_{1}and

*J*

_{2}: (A 7)

Impulse *J*_{1} is quantified from the vertical deceleration of *m*_{1} during surface impact:
(A 8)

where *m*_{1} is 8.0% of the body's mass, Δ*t*_{1} is the time interval between touchdown and the vertical velocity of *m*_{1} slowing to zero, Δ*v*_{1} is the change in vertical velocity of *m*_{1} during Δ*t*_{1}, and *F*_{1,avg} is the average force during the total time interval (2Δ*t*_{1}) of impulse *J*_{1}.

Impulse *J*_{2} corresponds to the remainder of the body's mass and is determined from *J*_{1} in Eqn A8 and total ground reaction impulse *J*_{T} in Eqn A3:
(A 9)

### Force curve function

The bell-shaped force curve *F*(*t*) for each impulse (*J*_{1}, *J*_{2}) can be accurately modeled using the RCB curve (Clark et al., 2014). The raised cosine function can be derived from the first two terms of the Fourier series:
(A 10)

where α_{0} is the mean of the signal, and *f _{n}*, α

_{n}and θ

*are the frequency, amplitude and phase angle of the*

_{n}*n*th harmonic (Clark and Weyand, 2014; Winter, 1990). The first two terms are: (A 11)

The peak of this function can be referenced to *t=*0 by defining phase θ_{1}=π/2:
(A 12)

Term α_{0} is the mean of the function, and term α_{1} is the amplitude of the function. Each term is defined by the total peak amplitude *A*, resulting in α_{0}=α_{1}=*A*/2. The peak is located at center time *B*. Frequency *f*_{1} can be expressed in terms of width parameter *C*, which is defined from the peak at *t*=*B* to the time where *F*(*t*) decays to zero, resulting in *f*_{1}=1/(2*C*). The constants *A*, *B* and *C* are inserted into Eqn A12 to yield the raised cosine periodic function:
(A 13)

The RCB curve is defined over a finite time interval of one period: (A 14)

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 simple properties of this function, peak amplitude *A*=2*F*_{avg}, and the area under the curve is *J*=*AC*.

Impulse *J*_{1} has force waveform *F*_{1}(*t*) during the interval *B*_{1}−*C*_{1}≤*t*≤*B*_{1}+*C*_{1}:
(A 15)

where *A*_{1}=2*F*_{1,avg} using *F*_{1,avg} in Eqn A8, and *B*_{1} and *C*_{1} equal the time Δ*t*_{1} after touchdown for the vertical velocity of *m*_{1} to reach zero.

Impulse *J*_{2} has force waveform *F*_{2}(*t*) during the interval *B*_{2}−*C*_{2}≤*t*≤*B*_{2}+*C*_{2}:
(A 16)

where *A*_{2}=2*F*_{2,avg} using *F*_{2,avg} in Eqn A9, and *B*_{2} and *C*_{2} equal 0.5*t*_{c} for a symmetrical waveform.

The total force curve *F*_{T}(*t*) is the sum of these two individual force waveforms:
(A 17)

### Impulse *J*_{1} ankle marker velocity considerations

Impulse *J*_{1} results from the acceleration of the lower limb during the limb impact interval and is quantified by Eqn A8. Δ*v*_{1} and Δ*t*_{1} are determined using ankle marker kinematics to represent the motion of the lower limb mass *m*_{1}. The lower limb attains a relatively constant positive velocity after the impact interval (Fig. 2C), resulting in a near-zero acceleration of *m*_{1} (Fig. 2D) and negligible force *F*_{1}=*m*_{1}** g** (Fig. 2E) as described in Eqn A6. This force is less than 55 N for a subject with body mass

*m*

_{b}=70 kg and lower limb mass

*m*

_{1}=5.6 kg (

*m*

_{1}=0.08×70 kg). Thus,

*J*

_{1}can be modeled by a finite impulse during the impact interval.

Kinematic data for the ankle marker position during the impact interval for representative RFS and FFS subjects appear in Fig. 3. Δ*t*_{1} is the time interval between touchdown and the vertical velocity of *m*_{1} slowing to zero. For some FFS runners at slower speeds, minor fluctuations in the ankle marker velocity–time profile near the end of the *m*_{1} impact interval can create variability in the Δ*t*_{1} measurements as a result of slow *m*_{1} impact velocities and skin marker motion artifact during the impact interval (Bobbert et al., 1991). Accordingly, Δ*t*_{1} was quantified using a technique that represented the functional end of the *m*_{1} impact time interval. An ankle marker velocity of −0.25 m s^{−1} was used as a threshold point, and the previous 10 ms of data were utilized for a linear projection of the ankle marker velocity to zero. For consistency in analysis, this method was applied to all ankle marker kinematics data, regardless of the speed or foot-strike mechanics of the runner (see Results).

### Impulse *J*_{2} asymmetry considerations

The temporal location of the peak of impulse *J*_{2} is dependent on the relative location of the center of mass at touchdown and takeoff. Idealized spring–mass running has symmetrical center of mass displacement, and thus a symmetrical profile where the location of peak *B*_{2} is 0.50*t*_{c}. However, the stance leg is more extended at takeoff than at touchdown (Blickhan, 1989), and this causes the center of mass height to be lower at touchdown than at takeoff (Cavagna, 2006), which results in an asymmetrical impulse *J*_{2} profile. The model waveform *F*_{2}(*t*) for impulse *J*_{2} can be modified to include width parameters to control the symmetry:
(A 18)

where *A*_{2} is the peak amplitude, *B*_{2} is the center time of the peak, *C*_{2L} is the leading half-width time interval, and *C*_{2T} is the trailing half-width time interval. The location of peak *B*_{2} was set at 0.47*t*_{c} as per the center of mass asymmetry originally reported by Cavagna et al. (1977) (see their table 4). With the symmetry control, *C*_{2L}=*B*_{2} and *C*_{2T}=*t*_{c}*−B*_{2}.

### Impulse *J*_{2} threshold considerations

A baseline noise level is present in all force platforms. To establish accurate and reproducible contact measurements, a threshold value is specified such that any ground reaction force signal below this value is zero. The threshold setting can be incorporated into the impulse *J*_{2} model waveform *F*_{2}(*t*). Eqn A18 can be solved for the width parameters *C*_{2}:
(A 19)

This equation is specifically evaluated for each width parameter at time *t* where the force *F*_{2} is equal to the threshold. A force threshold of 40 N was used for the AMTI high-speed force treadmill system. The leading width parameter *C*_{2}=*C*_{2L} is determined at the first channel (*t*=1 ms) and the trailing width parameter *C*_{2}=*C*_{2T} is determined at the last channel (*t*=*t*_{c}). Leading and trailing width parameters *C*_{2L} and *C*_{2T} are calculated using the same peak location *B*_{2} and peak amplitude *A*_{2}.

The peak amplitude is recalculated after the width parameters are changed in order to preserve the impulse. The impulse after the calculation of width parameters (*J*_{2A}) should approximately equal the impulse before the calculation of width parameters (*J*_{2B}):
(A 20)
(A 21)
(A 22)

As a result of this recalculation, the impulse is approximately equal to the original impulse and the values at the first and last channel are approximately equal to the threshold value.

## FOOTNOTES

**Competing interests**The authors declare competing financial interests. P.G.W., L.J.R. and K.P.C. are the inventors of US Patent no. 8363891 which is owned by Southern Methodist University and contains scientific content related to that presented in the manuscript. The patent is licensed to SoleForce LLC in which the three aforementioned individuals are equity partners.

**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 writing the manuscript.

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

- Received January 25, 2016.
- Accepted October 24, 2016.

- © 2017. Published by The Company of Biologists Ltd