The feasibility of predicting ground reaction forces during running from a trunk accelerometry driven mass-spring-damper model

Subjects and protocol

Twenty healthy male athletes (age 22 ± 4 years, height 178 ± 8 cm, mass 76 ± 11 kg) who engaged in running related sports activities on a weekly basis volunteered to participate in this study. The institutional ethics committee at Liverpool John Moores University granted ethical approval for this study (ethics approval number: 09/SPS/010) in accordance with the Declaration of Helsinki, and written consent was obtained from all participants. After a 15 min warm-up (including light jogging, dynamic stretching and individual dynamic tasks) and an individual number of familiarisation trials, the participants were asked to run over a force platform at different running speeds of 2, 3, 4 and 5 m s⁻¹ (±5%) in a randomised condition order. Running speeds were measured with photocell timing gates (Brower Timing System, Utah, USA) placed 2 m apart, with the last gate positioned 2 m before the centre of the force platform as described in Vanrenterghem et al. (2012). The participants completed four trials of each running speed, landing on the force platform with their dominant leg (defined as the self-reported preferred kicking leg (Van Melick et al., 2017)). Trials with unsuccessful foot contacts on the force platform (double foot contact or when the foot was not placed within the force platform) and/or when the desired approach speed was not achieved were repeated until a valid trial was recorded.

Experimental measurements

Resultant ground reaction forces were measured (GRF) with a sampling frequency of 3,000 Hz from a 0.9 ×  0.6 m² Kistler force platform (9287C, Kistler Instruments Ltd., Winterthur, Switzerland). Resultant trunk accelerations (TrunkAcc) were simultaneously collected at 100 Hz using a tri-axial accelerometer (KXP94, Kionex, Inc., Ithaca, NY, USA) with an output range of ±13 g embedded within a commercial GPS device (MinimaxX S4, Catapult Innovations, Scoresby, Australia) with a total weight of 67 grams and 88 × 50 × 19 mm in dimension. The GPS device was positioned on the dorsal part of the upper trunk between the scapulae within a small pocket of a tight fitted elastic vest according to the manufacturer’s recommendations (Boyd, Ball & Aughey, 2011). Different vest sizes were used to ensure the tightest fit for the individual participants. TrunkAcc data (measured in the units g) was pre-processed with the manufacturer’s proprietary filter algorithms (50 Hz low-pass filter, personal communication with the manufacturer), and downloaded as ‘raw accelerometer data’ from the manufacturer’s software (Catapult Sprint, version 5.1.7, Melbourne, Australia) after each test session. Each session also included a static measurement at the beginning and end of the session to detect any calibration drift over time, and none was detected. TrunkAcc and GRF were synchronised using a combination of time and event synchronisation as described in Nedergaard et al. (2017a) and exported to Matlab (version R2016a, The MathWorks, Inc., Natick, MA, USA) where a 4th order recursive Butterworth low-pass filter with a cut-off frequency of 20 Hz was applied to GRF and TrunkAcc. GRF data was collected from a single stance phase per trial, where touch down and take off were defined when the vertical GRF crossed a 20 N threshold.

Mass-spring-damper model

The complex multi-segmental dynamics of the human body during stance phase were modelled as a passive MSD-model (Alexander, Bennett & Ker, 1986; Derrick, Caldwell & Hamill, 2000). This model consists of two masses (Fig. 1); a lower point mass (m₂) on top of a linear spring (k₂) in parallel with a damper (c) representing the support leg; an upper point mass (m₁) representing the dynamics of the rest of the body and another linear spring (k₁) connecting the two masses.

The one-dimensional motion of the MSD-model was described by the acceleration of the two masses Eqs. (5) and (6):

(1)${\displaystyle \lambda =\frac{m_1}{m_2}}$ (2)${\displaystyle {\omega }_1^2=\frac{k_1}{m_1}=\frac{\left(1+\lambda \right)k_1}{\lambda M}}$ (3)${\displaystyle {\omega }_2^2=\frac{k_2}{m_2}=\frac{\left(1+\lambda \right)k_2}{M}}$ (4)${\displaystyle \zeta =\frac{c}{2\sqrt{k_2{m}_2}}}$ (5)${\displaystyle a_1=-{\omega }_1^2\left(p_1-p_2\right)+g}$ (6)${\displaystyle a_2=-{\omega }_2^2{p}_2+{\omega }_1^2\lambda \left(p_1-p_2\right)-2\zeta {\omega }_2^2{v}_2+g}$ (7)${\displaystyle GR{F}_{\mathrm{model}}=\frac{M{\omega }_2}{1+\lambda }\left({\omega }_2{p}_2+2\zeta v_2\right)}$ where p₁, p₂ v₁, v₂, a₁, and a₂ are the initial positions, velocities and, accelerations of the two masses (m₁ and m₂), respectively; λ is the mass ratio of the lower mass relative to the total body mass (Eq. 1); ω₁² and ω₂² are the natural frequencies of the springs (Eqs. (2) and (3)) based on the linear spring constants (k₁ and k₂) and the mass of the two masses (m₁ and m₂); ζ is the damping ratio of the damper (Eq. 4); and g is the acceleration from gravitational acceleration (−9.81 m  s⁻¹). The resultant GRF acting on the MSD-model is calculated as in Eq. (7), where M is the sum of the two masses (i.e., total body mass):

Model parameter estimation

To estimate the eight MSD-model parameters (p₁, p₂ v₁, v₂, ω₁² ω₂², ζ, λ), we used gravity corrected TrunkAcc from the stance phase as a surrogate of the MSD-model’s upper mass acceleration (Fig. 2A). For each trial, model parameters (ACC_param) were optimised by fitting the MSD-model’s upper mass acceleration (a₁) to the TrunkAcc signal. A purpose-built gradient descent optimisation routine in Matlab was used, where the two second-order differential equations of the MSD-model’s motion Eqs. (5) and (6) were transformed to four first-order differential equations and solved numerically with a Runge Kutta 4th order method. Root mean squared error (RMSE) between the TrunkAcc and a₁ waveforms were calculated for every iteration to determine the optimal model ACC_param combination that best fitted TrunkAcc for the individual trials. The ACC_param estimated from the TrunkAcc fitting were then used to predict the resultant GRF from Eq. (7). Furthermore, to help understand differences in estimated model parameters and the predicted versus measured resultant GRF, we also estimated the eight model parameters (GRF_param) by fitting the MSD-model to the measured GRF (Fig. 2B), similar to the approach previously described by Derrick, Caldwell & Hamill (2000).

Statistical analysis

Measured and modelled GRF were normalised to the participants’ mass. RMSE between the TrunkAcc and a₁, waveforms, and between the measured GRF and predicted GRF waveforms, were calculated to evaluate the accuracy of the TrunkAcc fitting and the predicted GRF, respectively. RMSE median and interquartile range (25th to 75th percentile) were calculated to determine the variation in the model’s accuracy within and across running speeds. Similarly, the median and interquartile range (25th to 75th percentile) of the ACC_param and GRF_param were calculated to explore the variation within and across running speeds. The median data presented and discussed in the following is the median of all trials within the individual running speeds (N = 80 trials) and the overall median across all running speeds (N = 320 trials).

Model parameter estimation

The eight model parameters are fundamental to calculating the resultant GRF acting on the MSD-model, and though fitting TrunkAcc was successful, the ACC_param estimated from this approach resulted in poor GRF predictions. Based on the equation of the upper mass acceleration (Eq. 5) and the ACC_param estimated from TrunkAcc, it seems that the MSD-model was able to fit the TrunkAcc by keeping the initial position of the upper mass (p₁) and lower mass (p₂) low, and by keeping the spring stiffness of the upper spring (ω₁²) low. Whereas p₁ has minor influence on the predicted GRF, the velocity of the upper mass at initial contact (v₁) is indirectly influenced by changes in the initial upper mass position ($v_1=\dot{p_1}$). Derrick, Caldwell & Hamill (2000) found that decreased v₁ has a large impact on the duration of the stance phase and therefore could have contributed to the overestimation of foot-ground contact (Fig. 3B). Similarly, the MSD-model decreased the spring stiffness of the upper spring (ω₁²) to better fit the two acceleration peaks typically observed in the TrunkAcc data, which has previously been shown to increase the duration of the stance phase (Derrick, Caldwell & Hamill, 2000). Furthermore, the MSD-model lowered the initial position of the lower mass (p₂), which previously has been shown to both increase the GRF at touch down and decrease the magnitude of the impact peak (Derrick, Caldwell & Hamill, 2000). We therefore believe that the high GRF values observed in our GRF predictions at touch down (Fig. 3B) were primarily related to the lower initial position of the lower mass (p₂) required to fit the upper mass acceleration to the TrunkAcc. Finally, the MSD-model also kept the damping ratio (ζ) low to better fit the magnitude of the two acceleration peaks in the TrunkAcc. Decreasing the damping ratio, has however previously been shown to increase the oscillation in the model’s GRF (Alexander, Bennett & Ker, 1986; Derrick, Caldwell & Hamill, 2000), and may therefore explain why our GRF predictions to a large extent include oscillating characteristics (Fig. 3B).

The comparison between the ACC_param estimated from the TrunkAcc and the GRF_param estimated from measured GRF, clearly demonstrates that the model is unsuitable for predicting GRF from TrunkAcc. A closer look at the GRF_param, showed that the median position and velocity of the lower mass (p₂ and v₂) was constant across running speeds and only varied marginally within running speeds (Fig. 4C). In addition, only small differences in median damping ratios (ζ) were observed between running speeds in this study (ζ between 0.31 and 0.39 au). It was in fact kept constant (ζ = 0.35 au) in the study by Derrick, Caldwell & Hamill (2000). Based on these observations we explored the effect of keeping p₂, v ₂, and ζ ACC_param constant for all trials (using the median GRF_param across running speeds), and for the remaining five MSD-model parameters use the trial specific ACC_param to re-calculate the predicted GRF (Fig. S1). Whilst this decreased the variability of the GRF prediction (RMSE interquartile range) both within and across running speeds, only minor improvements were observed in the GRF prediction. This indicated that keeping selected ACC_param constant would not substantially improve the GRF prediction in future studies. Furthermore, when selected ACC_param were kept constant, their original interaction was broken.

MSD-model hypothesis

If the trunk accelerometry data accurately represents the model’s upper mass acceleration one would at least expect that the ACC_param related to the motion and stiffness of the upper mass and spring (p₁, v₁, ω₁²) would be close to the GRF_param estimated when fitting measured GRF. This was however not the case, and therefore naturally raises the questions as to whether the upper mass acceleration is equivalent to the acceleration measured from trunk accelerometry during running. The trunk accelerometry driven MSD-model approach introduced in this study is based on the hypothesis that the model’s upper mass primarily represents the mass and motion of the trunk segment (Alexander, Bennett & Ker, 1986; Derrick, Caldwell & Hamill, 2000). Our results suggest however that this is not the case, and that independent accelerations of other body segments (e.g., the swing leg and arms) significantly contribute to the MSD-model’s upper mass accelerations. We therefore conclude that the primary model hypothesis for this study was false, and that trunk-mounted accelerometry alone is inappropriate as input for the MSD-model to predict meaningful GRF waveforms.

A high initial peak related to the attenuation of the shock impact from the foot’s collision with the ground (Hamill, Derrick & Holt, 1995; Derrick, 2004) dominated the TrunkAcc signals across running speeds. In contrast, a higher second peak related to the COM displacement during the stance phase dominated the upper mass acceleration when the MSD-model was fitted to measured GRF. This raised the technical question as to whether the poor GRF predictions observed from the measured accelerometer signal were partly a consequence of an artificially high frequency of that initial peak and whether the application of lower filter cut-off frequencies (cut-off frequencies of 20 Hz in the present study) would improve GRF predictions. To explore this, trunk accelerometry data of 10 representative participants was low-pass filtered with cut-off frequencies of 15, 10 and 5 Hz (Fig. S2). Whilst low cut-off frequencies (especially 10 and 5 Hz) to a large extent successfully removed the initial high-frequency peak in the accelerometry signal, and the RMSE between TrunkAcc and upper mass acceleration decreased, it only had a minor influence on the RMSE of the predicted GRF across running speeds (Fig. S2). Therefore, accelerometry post-processing did not improve the GRF predictions from TrunkAcc. This suggests that the trunk accelerometry signal in itself was not the main reason for the poor GRF predictions, but rather an incorrect hypothesis that the MSD-model’s upper mass acceleration primarily represents the acceleration of the trunk segment.

Replicating GRF from measured GRF

Although TrunkAcc was unsuccessful in predicting GRF during running with a simple MSD-model, the MSD-model could successfully replicate measured GRF during slow to moderate running speeds. In fact, the inclusion of all eight GRF_param in our optimisation routine, compared to only optimising the spring constants of the upper and lower springs (k₁ and k₂) and the position of the lower mass (p₂) (Derrick, Caldwell & Hamill, 2000) allowed us to replicate the measured GRF with higher accuracy. These findings illustrate that despite the MSD-model’s simplicity it has the ability to replicate and potentially predict GRF for a range of running speeds. Since the MSD-model parameters associated with the lower mass and spring are crucial to predict GRF (Eq. 7), this may open opportunities to use segmental kinematics and/or accelerometry from lower extremities to estimate MSD-model parameters. This does however require that the lower limb accelerations measured from e.g., a tibia-mounted accelerometer are similar to the MSD-model’s lower mass acceleration required to accurately predict GRF, something which is not a given. Recent studies have for example shown promising results in predicting GRF during sprinting, in high level sprinters, when contact and flight time, in combination with kinematics from the ankle were used as input for a two-mass model (Udofa, Ryan & Weyand, 2016; Clark, Ryan & Weyand, 2017). Future studies are however still need to explore the use of body-worn micro sensor technology to drive simple human body models and predict GRF waveforms for a range of running speeds.

Model limitations

A limitation with the MSD-model and the associated model parameters is that multiple parameter combinations exist when fitting the MSD-model to measured TrunkAcc or GRF waveforms. Whilst it could be of interest to further explore the physical meaning of the individual model parameters (ACC_param or GRF_param) and their interactions, or within and between subject parameters variations, this was not possible due to the existence of multiple model parameter solutions. Trunk-mounted accelerometry has a major benefit that it is already in use in many field contexts, but a limitation is that it may not very well represent the acceleration of the trunk segment. We have in previous work (Nedergaard et al., 2017b) shown that vertical trunk accelerations, measured from a high-end lab-based motion capture system, improved the upper mass acceleration fitting (median RMSE: 0.03 g across all running speeds) and lowered the average median RMSE of the GRF predictions to 5.18 N kg⁻¹ (vertical GRF) across all running speeds, compared to 8.99 N kg⁻¹ in the current study. Importantly, the accuracy and reliability of the GRF predictions are considered poor in both cases, suggesting that our hypothesis that the MSD-model’s upper mass acceleration primarily represents the trunk acceleration is most likely the weakest link. Secondly, the MSD-model is a one-dimensional model, and therefore only allows the magnitude of the resultant GRF to be estimated. We decided to predict the magnitude of the resultant GRF in our study, considering that we wanted to estimate the overall external biomechanical loading on the body, however we accept that others may prefer to predict the magnitude of the vertical GRF only. Ultimately, we believe that it is important to recognise that the MSD-model approach omits any direction specific load variations across running speeds, and that these may well be relevant in how the musculoskeletal tissues are exposed to stresses. Finally, the MSD-model is a passive elastic model and therefore does not account for additional energy generated by the body’s “active” structures (muscles). Whilst a more complex model could account for this (Zadpoor & Nikooyan, 2010; Nikooyan & Zadpoor, 2011), it is questionable if this would allow for better GRF predictions from TrunkAcc. The complexity of such model would probably also defeat the overall purpose of using a simple model that is still applicable in field settings.