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 (mb). 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, m1=0.08mb; remaining mass, m2=0.92mb). 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 (R2=0.95±0.04, mean±s.d.), with minimal variation across the slow, medium and fast running speeds tested (ΔR2≤0.04), and between rear-foot (R2=0.94±0.04, N=177) versus fore-foot (R2=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.

Multi-mass model and two-mass model. (A) Representative multi-element spring–mass-damper model (diagram adapted from Nigg, 2010). This type of model relates variations in ground reaction forces to mechanical characteristics of specific elements in the model. (B) The two-mass model is a computational model incorporating the known mass distributions of the human body as illustrated. Mass 1 (m1) represents the lower limb (8% total body mass, mb) and mass 2 (m2) represents the remainder of the body's mass (92% total body mass). Prior to touchdown, m1 typically has a greater vertical velocity than m2 as a result of the extension of the leg prior to impact.

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 (m1, 8% total body mass) and the second corresponding to the remainder of the body's mass (m2, 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 m1 with the running surface and impulse 2 results from the vertical accelerations of m2 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 (R2 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 (J1) with the running surface and the concurrent vertical accelerations of the rest of the body (J2) during ground contact. The total ground reaction impulse JT is the sum of J1 and J2 and can be determined by the total stance-averaged vertical ground reaction force FT,avg during the ground contact time tc:
(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 FT,avg can be determined if contact time tc and aerial time ta are known:
(2)

where mb is body mass and g is gravitational acceleration (g=9.8 m s−2), and the quantity tc+ta equals the step time, tstep.

Impulse J1 corresponds to the vertical deceleration of m1 during surface impact:
(3)

where m1 is 8.0% of the body's mass (Plagenhoef et al., 1983; Winter, 1990), Δt1 is the time interval between touchdown and the vertical velocity of m1 slowing to zero, Δv1 is the change in vertical velocity of m1 during Δt1, and F1,avg is the average force during the total time interval (2Δt1) of impulse J1, here defined as the impact interval. A single ankle marker is used to measure Δv1 and Δt1, which determine the vertical acceleration of the lower limb mass m1 (Fig. 2). The lower limb attains a relatively constant positive velocity after the impact interval (Fig. 2C), resulting in a near-zero acceleration of m1 (Fig. 2D) and negligible force (Fig. 2E). Thus, J1 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 Δt1 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.

Lower limb motion and force during ground contact. (A) A stick figure illustration of mass segment m1 motion (a–d) during the foot–ground portion of a running step. The red circle represents the axis of rotation of the ankle joint. (B–E) Corresponding schematic graphs for vertical position (B), velocity (C), acceleration (D) and force (E) of lower limb mass m1 versus time during the ground contact phase. After the impact interval, m1 reaches a relatively constant positive velocity, resulting in near-zero acceleration of m1 and a negligible force contribution from the lower limb mass for the remainder of the ground contact phase. A simplifying assumption of the two-mass model is that the force resulting from the acceleration of m1 can be accurately modeled by a finite impulse during the impact interval.

Ankle marker vertical position versus time data just prior to touchdown and during the initial ground contact phase for a representative rear-foot strike (RFS) and fore-foot strike (FFS) runner. After touchdown, the ankle marker decreases in vertical position until it reaches its lowest position above the running surface (Min.); the time interval for this deceleration is bracketed by Δt1 for the RFS and FFS, respectively.

Impulse J2 corresponds to the remainder of the body's mass and is determined from total ground reaction impulse JT in Eqn 1 and impulse J1 in Eqn 3:
(4)

where F2,avg is the average force of impulse J2 during the contact time tc.

The raised cosine bell (RCB) curve function was used to generate the F(t) waveforms of J1 and J2 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 cos2 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=2Favg 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 J2 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 J2 force peak corresponding to the minimum height of the center of mass was set at 0.47tc 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 FT(t) is the sum of the two individual force waveforms representing impulse J1 and impulse J2 :
(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, m1 vertical velocity at touchdown and impact interval time Δt1 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.

Measured and predicted force–time waveforms. (A) Five consecutive vertical ground reaction force–time waveforms from a RFS subject jogging at 4 m s−1. Force units were standardized to body weight (Wb). The first, third and fifth step illustrated were footfalls from the right leg. (B) The average values from the three right footfalls were used as input parameters to create a trial-averaged impulse 1 and trial-averaged impulse 2, which combined to form the trial-averaged model force waveform. (C) The model-predicted force waveform and measured average force waveform for the three right footfalls.

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 m1 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 m2 quantity of 98.5%. A larger m1 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 m2 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 m1 quantity of 1.5% and over-predicted with an m1 quantity of 16%.

Statistical analysis

The predictive accuracy of the model was assessed on the 500 individual footfalls acquired using both the R2 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 R2 and RMSE statistics.

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/√(gL0)>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 (R2) 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.2Wb and was equal to 11.5% of the mean stance-averaged vertical force from all 500 trials (1.82±0.23Wb). The R2 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 (R2=0.96±0.03) than the individual footfall value, while the RMSE of prediction was only 0.02Wb less (0.19±0.07Wb) 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 (R2) was nearly identical (ΔR2=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 R2 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.04Wb) 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 R2 values, when averaged for the waveforms from the three speed ranges (R2 total, Table 1) regardless of foot-strike type, were slightly, but significantly lower (ΔR2=−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 R2 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 m1 mass segment for RFS versus FFS (Δt1, Table 2) with little difference in the overall J1 impulse values (Fig. 5B,D,F versus Fig. 6B,D,F; Δv1, Table 2) at similar speeds. From slower to faster speeds, J1 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 J2 impulse predicted from the model because the m1 impact deceleration event concluded during the first quarter to half of the total contact period.

Vertical ground reaction force–time waveforms for RFS trials. Trials for an individual subject took place at 3.0 m s−1 (A,B), 6.0 m s−1 (C,D) and 8.0 m s−1 (E,F). A, C and E illustrate the model-predicted waveform (solid blue line) from the average kinematics measured during the trial compared with a measured average of the vertical ground reaction force waveforms (solid black line). B, D and F illustrate the sum of the first impulse (dotted red line) and second impulse (dashed green line) to form the total model-predicted waveform (solid blue line).

Vertical ground reaction force–time waveforms for FFS trials. Trials for an individual subject took place at 4.0 m s−1 (A,B), 6.0 m s−1 (C,D) and 11.1 m s−1 (E,F). A, C and E illustrate the model-predicted waveform (solid blue line) from the average kinematics measured during the trial compared with a measured average of the vertical ground reaction force waveforms (solid black line). B, D and F illustrate the sum of the first impulse (dotted red line) and second impulse (dashed green line) to form the total model-predicted waveform (solid blue line).

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 m1 segment values

Poorer agreement between model-predicted and measured waveforms resulted from using m1 segment values that were either smaller or larger than the originally assigned anatomical model value of 8.0% of mb. As hypothesized, waveforms generated using lesser m1 values (m1=1.5% with m2=98.5% of mb) 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 m1 values (m1=16.0% with m2=84.0% of mb) 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 m1 segment decelerations and correspondingly larger J1 impulses at faster speeds.

The vertical ground reaction force–time waveforms generated by the model for RFS trials with varied m1 values. For all panels, the solid black line represents the average of the vertical ground reaction force data measured during the trial, the solid blue line represents the model-predicted waveform from the average kinematics measured during the trial and the input m1 and m2 values, and the dotted red line represents the impulse resulting from the impact of m1 with the running surface. (A–C) Measured and predicted at 3.0 m s−1; (D–F) measured and predicted at 5.0 m s−1; and (G–I) measured and predicted at 7.2 m s−1. A, D and G illustrate waveforms predicted using m1 values of 1.5% total body mass and m2 values of 98.5% total body mass; B, E and H illustrate waveforms predicted using m1 values of 8.0% total body mass and m2 values of 92.0% total body mass; and C, F and I illustrate waveforms predicted using m1 values of 16.0% total body mass and m2 values of 84.0% total body mass.

The vertical ground reaction force–time waveforms generated by the model for FFS trials with varied m1 values. For all panels, the solid black line represents the average of the vertical ground reaction force data measured during the trial, the solid blue line represents the model-predicted waveform from the average kinematics measured during the trial and the input m1 and m2 values, and the dotted red line represents the impulse resulting from the impact of m1 with the running surface. (A–C) Measured and predicted at 3.0 m s−1; (D–F) measured and predicted at 5.0 m s−1; and (G–I) measured and predicted at 10.0 m s−1. A, D and G illustrate waveforms predicted using m1 values of 1.5% total body mass and m2 values of 98.5% total body mass; B, E and H illustrate waveforms predicted using m1 values of 8.0% total body mass and m2 values of 92.0% total body mass; and C, F and I illustrate waveforms predicted using m1 values of 16.0% total body mass and m2 values of 84.0% total body mass.

The R2 goodness of fit values for the 500 individual footfall waveforms generated using m1 values equal to 1.5% and 16.0% of mb accounted for 12% and 21% less variance, respectively, versus the original m1=8.0% results (m1−1.5% R2=0.83±0.16 and m1−16.0% R2=0.74±0.21). Predictive error values for the waveforms generated using the two alternative m1 segment values were approximately twice as large as those obtained using the original 8.0% value (m1−1.5% RMSE=0.36±0.17Wb; m1−16.0% RMSE=0.47±0.24Wb).

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 m2 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 m2 to be quantified without position and time data. We simply subtracted impulse J1 from the total ground reaction impulse JT, 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 J2 (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 m1's impact-period deceleration to a zero velocity (i.e. Δt1 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 Δt1 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 J1 predictions on Δt1 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 vankle=−0.05±0.04 m s−1). Given the overall mean ankle marker vertical velocity, Δv1, of −1.59±0.02 m s−1 at impact across all 500 footfalls (Δv1, Table 2), this method, on average, captured 97% of m1 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- (J1) and low-frequency (J2) components of the waveforms. The timing of impulse J1 was determined from motion data, and was therefore fully independent of our force data, filtering and processing routines. Conversely, the timing of impulse J2 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 J1 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 J1 illustrations (Figs 5 and 6) and the Δv1 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 m1 deceleration durations we measured. Longer FFS deceleration durations reduce and delay the J1 peak force values; they also result in greater waveform force contributions from impulse J2 at the time of the J1 force peak. In combination, these timing-dependent factors typically do not allow the J1 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: J1 impulses are ∝Δv1 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 J1 impulse peaks. The high pre-impact limb velocities (Δv1) and brief impact deceleration periods (Δt1) 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 m1 and m2 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 m1 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 m1 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 m1 fractional mass values. In both cases, the increases in the magnitude of the J1 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 m1 and m2 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 tc with vertical ground reaction force F(t), and an aerial time ta where the force is zero. Under these conditions, a runner supports an average of one body weight during the step time (tstep=tc+ta). 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 tstep:
(A 1)

where body weight is defined by the product of mass mb and gravitational acceleration g=9.8 m s−2. This equation yields the total average force FT,avg during the contact time interval:
(A 2)

The total ground reaction impulse JT can simply be determined from body weight and tstep:
(A 3)

The ground reaction force is a result of the acceleration ai(t) of body components i with mass mi contacting the ground:
(A 4)

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 Fi,avg and time interval. Impulse J1 results from the acceleration of the lower limb during the limb impact interval. Impulse J2 results from the acceleration of the remainder of the body's mass during the entire contact interval. The total ground reaction impulse JT is the sum of J1 and J2:
(A 7)

Impulse J1 is quantified from the vertical deceleration of m1 during surface impact:
(A 8)

where m1 is 8.0% of the body's mass, Δt1 is the time interval between touchdown and the vertical velocity of m1 slowing to zero, Δv1 is the change in vertical velocity of m1 during Δt1, and F1,avg is the average force during the total time interval (2Δt1) of impulse J1.

Impulse J2 corresponds to the remainder of the body's mass and is determined from J1 in Eqn A8 and total ground reaction impulse JT in Eqn A3:
(A 9)

Force curve function

The bell-shaped force curve F(t) for each impulse (J1, J2) 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 fn, αn and θn are the frequency, amplitude and phase angle of the nth 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 f1 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 f1=1/(2C). 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=2Favg, and the area under the curve is J=AC.

where A2=2F2,avg using F2,avg in Eqn A9, and B2 and C2 equal 0.5tc for a symmetrical waveform.

The total force curve FT(t) is the sum of these two individual force waveforms:
(A 17)

Impulse J1 ankle marker velocity considerations

Impulse J1 results from the acceleration of the lower limb during the limb impact interval and is quantified by Eqn A8. Δv1 and Δt1 are determined using ankle marker kinematics to represent the motion of the lower limb mass m1. The lower limb attains a relatively constant positive velocity after the impact interval (Fig. 2C), resulting in a near-zero acceleration of m1 (Fig. 2D) and negligible force F1=m1g (Fig. 2E) as described in Eqn A6. This force is less than 55 N for a subject with body mass mb=70 kg and lower limb mass m1=5.6 kg (m1=0.08×70 kg). Thus, J1 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. Δt1 is the time interval between touchdown and the vertical velocity of m1 slowing to zero. For some FFS runners at slower speeds, minor fluctuations in the ankle marker velocity–time profile near the end of the m1 impact interval can create variability in the Δt1 measurements as a result of slow m1 impact velocities and skin marker motion artifact during the impact interval (Bobbert et al., 1991). Accordingly, Δt1 was quantified using a technique that represented the functional end of the m1 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 J2 asymmetry considerations

The temporal location of the peak of impulse J2 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 B2 is 0.50tc. 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 J2 profile. The model waveform F2(t) for impulse J2 can be modified to include width parameters to control the symmetry:
(A 18)

where A2 is the peak amplitude, B2 is the center time of the peak, C2L is the leading half-width time interval, and C2T is the trailing half-width time interval. The location of peak B2 was set at 0.47tc as per the center of mass asymmetry originally reported by Cavagna et al. (1977) (see their table 4). With the symmetry control, C2L=B2 and C2T=tc−B2.

Impulse J2 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 J2 model waveform F2(t). Eqn A18 can be solved for the width parameters C2:
(A 19)

This equation is specifically evaluated for each width parameter at time t where the force F2 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 C2=C2L is determined at the first channel (t=1 ms) and the trailing width parameter C2=C2T is determined at the last channel (t=tc). Leading and trailing width parameters C2L and C2T are calculated using the same peak location B2 and peak amplitude A2.

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 (J2A) should approximately equal the impulse before the calculation of width parameters (J2B):
(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.

Similar articles

Subject collections

Other journals from The Company of Biologists

Neuropeptide evolution and function

Neuropeptides are a diverse assemblage of signalling molecules that have key roles in the regulation of behaviour. Understanding the evolutionary relationships and functions of the plethora of neuropeptides has presented a considerable challenge to biologists. Based on presentations and discussions at a Royal Society meeting in 2017, three companion Review articles by Elphick et al., Jékely et al. and DeLaney et al. discuss advances in our knowledge of neuropeptide evolution and function and the techniques that have facilitated progress in this field of research.

The exquisite bright colours of Pachyrhynchus weevils were thought by Alfred Russel Wallace to warn off potential predators, but whether this warning related to their hard exteriors, their spiky legs or some irritating taste had never been tested. Now, a century and half later, a team from Taiwan revisits this question and suggest that hardness itself acts as an effective secondary defence.

In our latest early-career researcher interview, Brooke Flammang, Assistant Professor at the New Jersey Institute of Technology, tells us about her research journey (including writing her Master's thesis in an ambulance while working as a paramedic), the importance of collaboration in integrative biology, and her approach to teaching.

"The paper provided the first quantitative field evidence of the way that animals might gain protection from predation by seeking cover in a group of other similar animals. This protection is known as the dilution effect."

William Foster discusses ‘Evidence for the dilution effect in the selfish herd from fish predation on a marine insect’, the 1981 classic he published in Nature with John Treherne, former JEB Editor-in-Chief.