Consumer-priced wearable sensors combined with deep learning can be used to accurately predict ground reaction forces during various treadmill running conditions

Data collection

Fifty distance runners (25 male and 25 female) were recruited to participate in this study if they typically completed at least one 30+ min run per week and were currently free of injury. The intention during participant recruitment was to include a representative pool of runners with a wide variety of anthropometric characteristics and performance capabilities (Table 1). All participants provided written informed consent and all experimental procedures were approved by the University of Bath’s Research Ethics Approval Committee for Health (EP20/21-069).

Before the start of data collection a pair of Nurvv Run insoles of the appropriate size was placed into the participant’s typical running shoe, placed underneath the shoes’ normal insole if possible. The Nurvv Run system (consumer price: $250 USD) is made up of a pair of thin pressure insoles that each contain 16 force sensitive resistors (FSRs) spread evenly across the foot and an attached device that clips on to the lateral aspect of the shoe. The device connected to the left insole contains an inertial measurement unit (IMU). During data collection the Nurvv system collected the FSR output value from all 32 sensors at the frequency of 1,000 Hz. Due to storage limitations on the Nurvv system these values were only logged at 50 Hz (20 ms moving average windows) between the identified initial contact and toe off times of each ground contact, values outside of these events were deleted. The IMU built into the device clipped to the lateral aspect of the left shoe collected and logged data continuously at 1,138 Hz (range of $\pm$4,000 degrees per second for gyroscope and $\pm$30 g for accelerometer). The footstrike classifications used within this article were exported directly from the Nurvv Run system, which uses a custom algorithm that compares the timing of pressure increases under each zone of the foot upon initial contact, informed by the results of Giandolini et al. (2014). The participants completed a running protocol on an instrumented treadmill that collected three-dimensional GRFs at 1,000 Hz (Bertec, OH, USA). The protocol was made up of a flat (+1%, chosen to reflect energetic cost of outdoor running (Jones & Doust, 1996)), uphill (+6%), and downhill (−4%) stage and the running speeds were almost all based around each participant’s typical easy training pace. During the flat condition participants completed 5 min at their easy training speed, 2 min 10% slower and 2 min at 10% faster than this speed. Participants also completed 3 min at the fixed speed of 12 km/h during the flat condition regardless of their chosen speed (three participants were not able to complete this speed). During the uphill condition, 2 min at their chosen speed and 2 min at 10% below their chosen speed was completed. During the downhill condition 2 min was completed at their chosen speed, as well as at 10% above and 10% below this speed.

Data processing

A fourth order low pass Butterworth filter (15 Hz cut off) was applied to the treadmill force data before a 40 N threshold was used to identify the timing of all foot touchdown and toe off events. Before any further processing of the pressure data was completed the logged 50 Hz FSR data was up-sampled back to the recorded frequency of 1,000 Hz using a quadratic interpolator in the Python package Scipy (v1.7.3). The pressure data (FSR signals) were then segmented into contact periods using an internal Nurvv algorithm, which uses the participant’s mass to calculate a total FSR output threshold that is then used to identify the same foot contact events. The same event timings were identified in the IMU data using an accelerometer-based step detection algorithm previously validated for treadmill running (Benson et al., 2019). As all three systems were running on separate clocks, the segmented foot contact data had to be synchronised to ensure the same physical foot contact was being evaluated by each independent system. This was done by correlating the stride times from two systems against each other. The alignment of these two vectors were then adjusted to find the strongest correlation between the sets of stride times calculated from each system, indicative of the correct alignment (Fig. 1). The alignment was accepted if the Pearson correlation coefficient was over 0.85. However, the correct alignment was often very clear ( $r>0.99$) in comparison to the incorrect alignment $(r<0.50)$ occurring one contact before and after. This process was performed once with stride times from the force data and the pressure data, before then being repeated with the force data and IMU data, ensuring foot contacts from all three systems were synchronised.

Each isolated foot contact had a varying number of data points due to the duration of the ground contact, across all participants the longest contact was 361 ms. Therefore, we normalised all signals to the length of 400 data points, to ensure no resolution was lost from the time series signals. The treadmill force data and the foot IMU data were filtered with fourth order low pass Butterworth filters with a 15 Hz and a 10 Hz cut off, respectively. As the pressure data was logged at 50 Hz (using 20 ms window averaging) before it was then interpolated back to 1,000 Hz, following visual inspection it was deemed that no filtering of the FSR time series data was needed. Based on the pressure data from each of the 16 FSR channels under the left foot, a two-dimensional centre of pressure coordinate ( $Co{P}_x$, $Co{P}_y$) throughout each foot contact was calculated using a weighted average approach and the known geometry of the sensors within the insole.

Following the preparation of the time series data at each foot contact, a data matrix had to be created that would serve as the input to the deep learning model for the estimation of the GRF traces. These input features should include as much information about the foot contact as possible whilst remaining feasible to collect continuously during an individual’s running session. The segmented and normalised time series data from the 16 FSR pressure channels, the two calculated centre of pressure values ( $Co{P}_x$, $Co{P}_y$), and the three axes of the IMU accelerometer were combined with the following discrete characteristics; duration of ground contact, participant mass, current running speed, current gradient/slope, insole length. All of these input features were chosen due to their feasibility to be collected outside the lab. The discrete values were duplicated for the length of the time sequence $[1\times 400]$. Each of these features were then concatenated together into a single input matrix of shape $[N\times 400\times 26]$ where N is the number of foot contacts included. Before this input matrix was passed to the model, the data was first standardised to z-scores. This ensures that features with greater magnitudes did not inherently have greater influence on the model’s output, as well as improving numerical stability and convergence speed during model training. The training set was first transformed to z-score values based on each features mean and standard deviation values, then the same transformation was performed to the validation set using the descriptive statistics calculated from the training set. This avoids the ‘leaking’ of any information from the validation set into the training set, replicating the pre-processing conditions upon model implementation.

Machine learning models

Independent deep learning models were created to estimate $GR{F}_v$ and anteroposterior ground reaction forces ( $GR{F}_{ap}$) from the input matrix in the Python package Pytorch (v1.10) (Paszke et al., 2019), but the same architecture was used to create each model. The matrix of input data was firstly passed to a bi-directional LSTM layer (Hochreiter & Schmidhuber, 1997), which output two hidden state vectors (one for each direction) at each time point in the sequence. These hidden state vectors are the same size as the hidden size of the LSTM network and contain temporal information about the input sequence from the forward and backward pass of the LSTM. The hidden state value at each time point was then passed through three linear layers, transforming the hidden state vector at that time point into an estimated force value (Fig. 2). This is repeated for the LSTM output at every timepoint in the sequence until a full time series of estimated force values are produced. This estimated force trace that the deep learning model produces was then compared to the measured force trace from the treadmill, the PyTorch autograd engine can then leverage its gradient tracking functionality to back propagate through the LSTM model, identifying how the model weights can be appropriately adjusted to decrease the difference between the measured and estimated force trace. The direction and strength of these model weight updates were informed by an Adaptive Moment Estimation (ADAM) optimiser (Kingma & Ba, 2014), aiming to minimise the RMSE between the measured forces and estimated forces across the training batch.

To increase the objectivity of the model creation process we chose to use hyperparameter optimisation to identify which hyperparameter choices, used during initial model definition, led to the highest GRF estimation accuracy. A search window for each hyperparameter was chosen and the Python package Optuna (v3.2, Tree-structured Parzen Estimator) was used to find the combination of hyperparameter options that maximised model accuracy. Running a hyperparameter optimisation on the full dataset can lead to the overfitting of hyperparameters to the current dataset, serving to inflate the model’s performance on the dataset within this study but reducing its ability to perform well on unseen data. To evaluate for the presence of this effect, a representative 25 of the 50 participants were separated out based on their anthropometric and running characteristics and split into five groups of five. The Optuna optimiser aimed to select the hyperparameters that minimised the average validation loss across the five folds, set to 250 iterations as a compromise between computational time and exploration of solution space. Within the presentation of the results the 25 participants included in the hyperparameter optimisation will be separated from those who were not, so the presence of hyperparameter overfitting can be evaluated for. After the final hyperparameters had been identified by the optimiser for the $GR{F}_v$ and $GR{F}_{ap}$ model they were used to construct the final model architecture. To ensure an accurate evaluation of the model’s performance on unseen data all data belonging to an individual remained contained in either the training or validation set at any given time. Whilst more computationally intense, Leave One Subject Out (LOSO) validation was performed. This evaluation technique keeps the data for all but one participants within the training set before evaluating the trained model’s performance on the individual’s data left out of training, repeated as many times as there are participants. In comparison to cross fold evaluation, this method provides greater insight into the current model limitations, setting expectations for future use of the trained model. To evaluate how accurately the model was able to estimate force traces at a range of magnitudes a Bland Altman analysis was performed.

To gain a greater understanding of which variables being input into the model were having the largest influence on the accuracy of the model output, Permutation Feature Importance (PFI) was also calculated. The PFI value is calculated by replacing the time series data of a given input variable with a sequence of random numbers, before passing the input matrix to the model (Molnar, 2020). The differing estimated force traces can then be compared, with the PFI value being the ratio between the RMSE after permuted input and the original RMSE value. This process was repeated 100 times for each variable in the input set to minimise the influence of random input values. Larger PFI values indicate a variable having a greater positive influence on the model’s force estimation accuracy, with values under 1.0 indicating the variable having a negative influence on the model’s accuracy.

Vertical ground reaction force estimation

Across all 50 participants the average RMSE was 0.15 bodyweight (BW) (Fig. 3) and the average correlation between measured and estimated $GR{F}_v$ traces was 0.991. The accuracy with which the model could estimate the $GR{F}_v$ trace varied considerably between participants. For the participant on which the model performed best, the average $GR{F}_v$ estimation error across all foot contacts was 0.05 BW, with an average correlation of 0.999. In comparison, for the worst performing participant the average RMSE was 0.35 BW with an average correlation of 0.979. The inclusion of a participant within hyperparameter optimisation did not appear to considerably influence model performance, with just a 0.01 BW difference between the average RMSE of the two groups.

The mean relative RMSE (rRMSE) (RMSE normalised to the range of the gold standard waveform) across all participants was 3.2%, the mean percentage difference in peak $GR{F}_v$ value was 6.3%, and the mean percentage difference in $GR{F}_v$ impulse was 5.5% (Fig. 4). The majority of participants that the $GR{F}_v$ estimation model achieved low rRMSE values for also had low error in the estimation of peak vertical force values and area under the force trace.

Figure 5 shows examples of individual foot contacts across the three gradient stages of the treadmill protocol and across three representative participants. The individual contacts were chosen based on the RMSE of that contact being closest to the mean RMSE for that participant and protocol section. Most notably, it is clear that with the participant that the model performs worst on, the model consistently under estimates the peak $GR{F}_v$ across all conditions.

Anteroposterior ground reaction force estimation

The best combination of hyperparmaters for the $GR{F}_{ap}$ model was used to create the final architecture for the $GR{F}_{ap}$ estimation model and LOSO on all 50 participants was performed. Across all participants the average RMSE was 0.04 bodyweight (BW) (Fig. 6) and the average correlation between measured and estimated $GR{F}_{ap}$ traces was 0.982. The accuracy with which the trained model could estimate $GR{F}_{ap}$ appeared to be more consistent than the $GR{F}_v$ estimation results. The same participant, who typically landed very far forward on the forefoot (P24), who had the lowest model performance for both $GR{F}_v$ estimation also had one of the highest levels of error for $GR{F}_{ap}$ estimation. However, it was not consistently the case that those participants that received poorer estimation accuracy from the $GR{F}_v$ model also performed poorly during $GR{F}_{ap}$ estimation.

The mean relative RMSE (rRMSE) across all participants was 3.1%, the mean absolute difference in the peak braking $GR{F}_{ap}$ value was 22.1 N (0.03 BW), and the mean absolute difference in the peak propulsion $GR{F}_{ap}$ value was 19.7 N (0.03 BW) (Fig. 7). It was not the case that as rRMSE increased the error in the estimated peak braking and propulsive forces also increased, which suggests that temporal alignment and estimated curve shape varied between participants. Figure 8 gives examples of how these individual contacts can be constructed. As would be expected the estimated $GR{F}_{ap}$ traces for the participant which the model performed best on showed excellent temporal alignment and shape replication, matched with very close peak braking and propulsive forces during the foot contact. In contrast, the worst performing participant had estimated $GR{F}_{ap}$ traces with poor temporal alignment and atypical signal shape just prior to midstance.

Model sensitivity

Due to the similarity of GRF traces during running, a machine learning model could feasibly tend to output an ‘average’ force trace without much sensitivity to the input values, and still achieve relatively high levels of accuracy. It was seen that during $GR{F}_v$ estimation there was a relatively consistent trend whereby the magnitude of force was overestimated at lower force values and underestimated for larger force values (positive correlation of the line of best fit) (Fig. 9). This finding does indicate that the output of the $GR{F}_v$ estimation models do have a tendency to ‘regress towards the mean’ at the extremities of the training data. This behaviour could be expected within deep learning models due to the loss minimisation approach that is used within their training. However, there are many different techniques that are employed to try and minimise this behaviour. Whilst the identification of this behaviour is important, understanding of the magnitude of it has much greater practical implications. The gradient of the line of best fit across the points presented in Fig. 9 can help quantify this regression magnitude. For example, for the presented model, it could be could be concluded that for every 1 BW increase in peak $GR{F}_v$ above 2.5 BW there would be a 0.3 BW underestimation on average.

The presence of ‘regression towards the mean’ behaviour was less apparent in the estimation of peak $GR{F}_{ap}$ braking force (Fig. 10A), with the error in the estimation of peak braking remaining relatively consistent across the range of magnitudes. However, for those contacts with higher propulsion forces their was a tendency for the model to underestimate their magnitude (Fig. 10B).