Abstract

Background. Detecting fatigue at the early stages of a run could aid training programs in making adjustments, thereby reducing the heightened risk of injuries from overuse. The study aimed to investigate the effects of running fatigue on plantar force distribution in the dominant and nondominant feet of amateur runners. Methods. Thirty amateur runners were recruited for this study. Bilateral time-series plantar forces were employed to facilitate automatic fatigue gait recognition using convolutional neural network (CNN) and CNN-based long short-term memory network (ConvLSTM) models. Plantar force data collection was conducted both before and after a running-induced fatigue protocol using a FootScan force plate. The Keras library in Python 3.8.8 was used to train and tune deep learning models. Results. The results demonstrated that more mid-forefoot and heel force occurs during bilateral plantar and less midfoot fore force occurs in the dominant limb after fatigue (). The time of peak forces was significantly shortened at the midfoot and sum region of the nondominant foot, while it was delayed at the hallux region of the dominant foot (). In addition, the ConvLSTM model showed higher performance (Accuracy = 0.867, Sensitivity = 0.874, and Specificity = 0.859) in detecting fatigue gait than CNN (Accuracy = 0.800, Sensitivity = 0.874, and Specificity = 0.718). Conclusions. The findings of this study could offer empirical data for evaluating risk factors linked to overuse injuries in a single limb, as well as facilitate early detection of fatigued gait.

1. Introduction

The popularity of long-distance running as an easily accessible and promoted sport has increased within the last four decades. However, the incidence of musculoskeletal injuries caused by running has also increased rapidly, especially in the lower extremities [1]. The tendons, muscles, or bones of the lower extremities during long-distance running are repeatedly subjected to chronic submaximal loading over a long period of time [2]. As a consequence, overuse injuries are considered to be the most common running injuries [2]. It should, however, be noted that the etiology of running fatigue-induced injuries is multifactorial and complex [3]. Fatigue from long-distance running can shift foot mechanics, potentially causing structural overload [4, 5]. Investigating the relationship between fatigue and the load distribution pattern during running gait has garnered increased interest [6]. A consensus is that an increase in peak metatarsal head pressure occurs after running fatigue [7, 8]. However, inconsistent results have been reported for the influence of running fatigue on the middle foot and heels [2, 7]. Weist et al. [9] demonstrated a significant increase in midfoot pressure and the impulse in the medial heel after performing a running-induced fatigue protocol. Nevertheless, the study by Bisiaux and Moretto [10] found a reduction in pressure and impulse in the midfoot under similar conditions. Willson and Kernozek [11] observed a significant reduction in peak heel pressure after 30 min of high-intensity running.

In addition, more laterally directed roll-off and inadequate pronated heel strikes have been demonstrated to be the potential triggers for lower limb overuse injuries [12]. Previous studies have reported that midfoot, metatarsal, and medial heel loading increases after running-induced plantar muscle fatigue, while loading on the lateral toes decreases [2]. Similar studies have reported that running fatigue causes a reduction in the medial longitudinal arch, which significantly increases mid-toe pressure [13]. In addition, excessive plantar forces of the forefoot lateral were identified as a potential cause for gait-related Achilles tendinopathy [7]. Most of the studies mentioned above did investigate changes in unilateral plantar pressure distribution after fatigue, especially the dominant side. However, there seems to be a lack of empirical data on the effect of running fatigue on the nondominant plantar load distribution. Due to the altered symmetry after fatigue, excessive loading usually occurs in the unilateral limb, especially in the nondominant limb, which may have a weaker fatigue tolerance. Ignoring asymmetric information about bilateral limbs to explore risk factors for fatigue, although it may simplify the data processing and analysis process, may also produce deceptive results. Assuming that changes in bilateral plantar pressure distribution play a primary role in the course of overuse injury development [14, 15], the role of fatigue on both dominant and nondominant plantar distribution should be examined. Therefore, one of the objectives of this study was to determine the effect of running-induced fatigue on bilateral plantar force distribution.

Previous studies have widely shown that fatigue gait risk is associated with shifts in the distribution of bilateral plantar pressure [1, 4, 6, 11]. Coaches and runners can avoid the occurrence of overuse injuries by monitoring fatigue levels in the context of competitive and recreational sports. In addition, excessive fatigue may affect runners’ performance and cause secondary injuries to the runner [16]. Therefore, human activity recognition (HAR) methods based on wearable sensors and deep learning algorithms have been widely developed in the last decade [1719]. Despite significant strides in gait and biomechanics analysis, research into automated fatigue gait recognition with data-driven models remains insufficient [18, 20]. Typical techniques to detect fatigue are surface electromyogram-based collection of muscle activity signals and optical motion capture-based collection of joint kinematics [21, 22]. However, the limited data collection area and the location of the marker attachments make monitoring limited. On the contrary, force plates or insoles with force sensors are easy to use and save time in the experimental setup for data collection. Therefore, this study intends to use a deep learning algorithm based on bilateral plantar pressure data for the early identification of fatigue gait.

Since larger spatial dependencies exist in the pressure data of each plantar region throughout the gait cycle [23, 24], the CNN model has been reported to be better at extracting local spatial features [15]. Similarly, time series data-based plantar pressure data are considered to possess time dependence [25]. However, the plantar pressure distribution data based on time series features may be regarded as static spatial data by the CNN model, and the time-dependent information within the series is lost. Previous studies have shown that long short-term memory network (LSTM) models perform better for the prediction of long-time dependence and nonlinear dynamic changes in a time series [26]. However, LSTM models are less effective in handling spatial relationships of data. The spatial characteristics of the pressure distribution in different plantar regions and the dynamic time characteristics of the variation with time should be considered in the model selection for this study.

The ConvLSTM model will be used in this study on the ground that it transforms the structure of recurrent neural networks into a convolutional structure, thereby preserving the spatial and temporal information of plantar pressure [27]. To verify the performance of the ConvLSTM model for fatigue gait recognition, we used a CNN model to compare the performances. Two hypotheses were proposed: (1) The metatarsal, midfoot, and heel pressures increased in the dominant and nondominant feet after the fatigue intervention, with more significant changes in the nondominant foot; (2) The ConvLSTM model has better performance than the CNN model for automatic recognition of fatigue gait.

2. Materials and Methods

2.1. Participants

Thirty healthy amateur runners (males) were enlisted from universities and local running clubs for this study. The anthropometric information of the participants is presented in Table 1. The inclusion criteria for this study were that the dominant extremity side was the right extremity side (preferred leg when kicking a ball), the absence of any lower extremity or pelvic musculoskeletal pain in the last 6 months, and running at least 2–3 times per week and for <45 min or <10 km at per running event. The Ethics Committee at Ningbo University approved the study (code: RAGH20210827), and all subjects signed the informed consent.

2.2. Data Collection

Subjects were guided by the experimental operator to familiarize the experimental environment (includes ground running tests with barefoot before and after the running-induced fatigue protocol) and process and participated in a 10-min jogging warm-up on a treadmill (Satun h/p/cosmos, Nussdorf–Traunstein, Germany) in advance. A previously identified and validated protocol was employed for building a running-induced fatigue model [28]. With reference to our previously built approach [29], a heart rate sensor band (Polar RS100, United States) and a Borg [30] RPE scale (6–20 scales) were utilized for monitoring fatigue during running. Every participant commenced walking on a treadmill at a velocity of 6 km/h. The pace of gait was augmented by 1 km/h every 2 min until an exertion level of 13 on the Borg scale was attained. Participants sustained the running pace at the established equilibrium until achieving a Borg rating of 17% or 90% of their maximum heart rate (maximum heart rate = 220-age), at which juncture they persisted in running for an extra 2 min. New neutral running shoes were given to every participant for the protocol involving running-induced fatigue.

Pedobarographic data collection was done before and after the running-induced fatigue protocol. Dynamic plantar force data were measured during running using a FootScan pressure plate (size: 2 × 0.4 m, frequency: 480 Hz, RsScan International, Olen, Belgium) embedded in the middle of a 20-m runway. The pressure plate is calibrated using the individual’s body weight prior to measurement to avoid errors. Two sets of infrared photocells were placed on either side of the data collection area to monitor the running speed. All participants were required to run barefoot over the data collection area at a speed of 3.3 m/s ±5% [30]. Participants were instructed to use the nondominant foot as the first step on the force plate and to ensure that two consecutive steps were recorded for each trial. Attempts to change the operating mode to strike the pressure plate were ruled out until three valid trial data points were measured before and after the running-induced fatigue protocol.

2.3. Data Processing

For each trial, 10 plantar anatomical regions were identified by the FootScan application. To avoid recognition errors, the pixels of each area were manually calibrated by an operator. These areas were defined as the hallux (H), other toe (OT), metatarsal 1–5 (M1–M5), midfoot (MF), medial heel (HM), and lateral heel (HL). Time-series attributes of force information for each region and the sum area were interpolated to 101 frames using linear interpolation for statistical comparison. To reduce the effect of individual weight and gait speed differences on the data, all data in this study were annotated using Zavg (total force over the entire support period divided by the total number of frames) [25]. As shown in Figure 1, to preserve asymmetric information of bilateral limbs before and after fatigue, the plantar force data of the nondominant and dominant foot were stitched longitudinally to obtain the bipedal force distribution information of one gait cycle for machine learning training [23].

2.4. CNN Model Building

This study uses the Keras Application Programming Interface (API) in Python 3.8.8 for CNN and ConvLSTM model building. CNN models have good performance for feature extraction of input data through convolutional operations of different topological structure kernels. The convolution layer in the model preserves the spatial relationships of the data by using the same convolution operation for each position of the original data. Each type of feature that is extracted generates a feature matrix Z. Therefore, after k times convolution calculations, the corresponding output matrix Zk can be represented by Equation (1). In addition, the convolution operation for one-dimensional time series data is also a nonlinear transformation of the original series. Applying a convolution kernel of length l to a univariate time series X of length T, Equation (2) is obtained.where Wk and k are the convolution kernels (size: k1 × k2) and the number of convolution kernels, respectively. b is biased, and the convolution operator is defined as . f is the activation function that performs a nonlinear transformation in the convolution layers.

As shown in Figure 2, The optimal CNN model for the recognition of fatigue gait is obtained through repeated debugging parameters. We used a total of eight convolutional layers, three maximum pooling layers, one average pooling layer, one dropout layer, and one dense layer to build the convolutional neural network model. The number of convolution kernels is set to (128, 128, 128, 128, 64, 64, 32, 32). The time step settings are (10, 10, 10, 10, 10, 10, 4, 4). In addition, “RELU” and “Softmax” are set as the activation functions for the convolutional and dense layers, respectively.

2.5. ConvLSTM Model Building

The convolutional layer extracted the temporal characteristics from the pressure data, while the LSTM layer handled the spatial characteristics (Figure 3(b)). In our ConvLSTM model, operations are depicted by Equations (3)–(8), where symbolizes the convolution process, and denotes the Hadamard product.

where it, ft, and ot are the input gate, oblivion gate, and output gate, respectively, in the proposed model; xt represents the data input at the current moment, while ht−1 refers to the output from the hidden layer at the preceding moment. ct denotes the cell state.

Figure 3(a) shows the framework for building the ConvLSTM model in this research. In this study, we try to make l choose a variety of different division lengths, such as 51, 101, 151, 202, and so on, for modeling, and finally find that the model’s classification performance is optimal when the subsequence length is l = 101. The optimal ConvLSTM for the recognition of fatigue gait is obtained through repeated debugging parameters. We sequentially set up a ConvLSTM layer (number of convolution kernels = 64, kernel size = (1,5)), a dropout layer (random discard ratio = 0.5), a Flatten layer and two dense layers (first: units = 50, activation = RELU”; second: units = 2, activation = “Softmax”) in the final model.

To ensure fast convergence during the training of the binary classification model, the cross-entropy loss function was chosen as the loss function for this study, as shown in Equation (9).where N is the number of samples and and are defined as the true and predicted values.

2.6. Statistical Analysis

In this study, a total of 90 cases were sampled, and 80% of the samples were set as the training set and 20% as the test set, where 20% of the training samples were set as the validation set for cross-validation. Therefore, the training set, validation set, and test set samples in this study are 72, 14, and 18. To avoid the occurrence of model underfitting, the number of model iterations was set to 300. This study uses accuracy, sensitivity, and specificity as quantitative metrics for the performance of two classification models. We used fatigue gait as positive samples and normal gait as negative samples. Therefore, Equation (10) was used to assess the overall classification capability of the models. Equations (11) and (12) were used to evaluate the classification capability of negative samples and positive samples of the models, respectively.where TP and TN are the number of samples correctly identified as fatigue gait and normal gait, respectively, and FP and FN are the number of samples incorrectly identified as fatigue gait and normal gait, respectively. To avoid accidental error, each model is run five times on the test set, and the corresponding classification results are collected.

The Shapiro–Wilk test was performed to check the normality of the data distribution. The paired sample T-test of open-source statistical parameter mapping 1d (SPM1d) was used to check the differences between pre- and post-fatigue time-series forces at dominant and nondominant foot. The discrete values of the percentage of time of peak force were checked using a paired sample T-test in Python 3.8.8 with the SciPy library. The significance levels were set at 0.05.

3. Results

3.1. Force Development in Toe and Metatarsal Areas

As shown in Figure 4, starting from initial nondominant foot contact, the force progression in forefoot regions differed between pre- and post-fatigue states. Specifically, the force in M3 shows a significant increase at 12%–79% of contact duration after fatigue. However, the force of M5 has decreased at 50%–69% of contact duration after fatigue (). For the dominant foot, there was a significant increase of force at OT (83%–95% (), 96%–100% () of contact duration), M2(17%–97%, ), and M3(14%–97%, ) after running-induced fatigue. However, the force of M4 decreased at 0%–3% of contact duration after fatigue ().

3.2. Force Development in the Middle Foot, Heel, and Sum Areas

As shown in Figure 5, there was no difference in MF at the nondominant foot. Interestingly, there was a significant decrease in plantar force at the dominant foot (0%–65%, ). The heel regions were directly affected by running fatigue. For nondominant and dominant plantar forces, they were significantly increased at HM (nondominant: 30%–36%, ; dominant: 11%–49%, ) and HL regions (nondominant: 11%–19%, ; dominant: 3%–51%, ). However, the sum of forces from all ten regions at nondominant (33%–46%, ) and dominant (34%–47%, ) significantly decreased after running-induced fatigue.

3.3. Relative Time of Peak Force

As shown in Table 2, the relative time of peak force was significantly shortened at MF () and SUM () regions at nondominant feet in a fatigued state. Similarly, there was a significant shortening in the relative time of peak force at M5 of the dominant foot after fatigue. Interestingly, for H regions, the relative time of peak force was significantly delayed.

3.4. Representations of Deep Learning Models

The classification results of total plantar pressure at CNN and ConvLSTM model are shown in Figure 6, and the confusion matrix and ROC of each model are shown in Figure 7. Table 3 presents the average accuracy, specificity, and sensitivity derived from the five test sets. The ConvLSTM model outperformed the CNN with an accuracy of 86.7% versus 80%. Likewise, ConvLSTM’s specificity was superior at 85.9%, compared to CNN’s 71.8%. Nonetheless, both models matched with a sensitivity rate of 87.4%.

4. Discussion

The aim of this investigation was to investigate the effect of running fatigue on the bilateral plantar force distribution of the foot and the effectiveness of CNN and ConvLSTM models for fatigue gait recognition. The results of this study showed that running fatigue changed the distribution pattern of load on the plantar of the dominant and nondominant limbs. These changes are similar to previous studies [6, 8, 11]. The force distribution of the dominant plantar of runners shown major differences reflected in reduced force under the midfoot at the expense of increased force under the H, M2, and M3. This increased loading of the medial forefoot region is in agreement with previously demonstrated higher pressures under the forefoot and lower peak pressures under the midfoot, which were reported by Bisiaux and Moretto [10] after fatigue induced by an intensive 30-min run. These results may indicate that the load was transferred from the midfoot to the toes and metatarsals [10]. The increased loading in M2–3 may be related to reduced activity of the toe flexors and posterior tibial muscles after running fatigue [31]. In addition, Arndt et al. [32] have also reported that higher strain rates and deformation of metatarsal bones can also occur after muscle fatigue caused by running. These findings could be a risk factor for a metatarsal stress fracture [10]. Especially the M2 and M3 are vulnerable because of the difference between the applied plantar pressure and bone strength [31].

Previous studies have demonstrated that dominant feet play a propulsive role, while nondominant feet are more likely to function as a stable gait [15]. Excessive force at the H region after fatigue appears in the dominant limb may be a compensatory effect of the functionally driven winch mechanism [5]. Willson and Kernozek [11] reported that running fatigue could cause changes in plantar surface loading characteristics and running technique. This study showed that the force of M5 at the nondominant foot has decreased at the metaphase (50%–69%) of contact duration, while the force of M3 has increased significantly at most of the contact duration (12%–79%) after fatigue, suggesting that the fatigue transferred foot loading from the lateral region toward the inside of the foot, especially in the nondominant foot [13]. This finding may be a weakening of the function of the nondominant limb to stabilize gait after muscle fatigue [21]. In addition, the relative time of peak force of MF at the nondominant was significantly shortened after fatigue, suggesting that more impulse was concentrated in the MF region. This finding may be due to the damage to the active control mechanism of the MF during the contact stage, leading to a reduction in the cushioning function of the nondominant plantar, which was the potential factors for plantar fasciitis [33]. Interestingly, the relative time of the peak of H-region force at the dominant foot was significantly delayed after fatigue. This finding may be a compensatory mechanism to maintain the propulsive function of the dominant limb, making the gravitational torque between the heel and toe region more even in the later stages of push-off [6].

Several reports have investigated the influence of the range of motion in the coronal plane of the foot on shock attenuation at heel strikes [34, 35]. In addition, the ability of the musculoskeletal system to attenuate the shock magnitude generated during heel strikes also decays with fatigue [36]. In our study, the plantar forces recorded under the HM and HL regions of both dominant and nondominant revealed the changes which running fatigue during the loading stage. Without other direct measurements, we can only speculate that excessive heel loading after fatigue may be linked to weaker muscle strength, which controls the movement of the ankle joint in the coronal plane after fatigue [10]. These observations were consistent with several previous studies [2, 7]. Interestingly, the sum of forces from all ten regions at nondominant of 33%–46% of contact duration significantly decreased, and the relative time of peak force was significantly shortened after running-induced fatigue, suggesting that dorsiflexor fatigue led to more vertical loading rate on the plantar. A significant interaction between loading rate and running-related calf, foot, and ankle injuries was demonstrated in a study by Gerlach et al. [37]. The results of this study could provide necessary enlightenment about the condition of different running-related injuries among runners with limbs on different sides.

In addition, by applying the feature set of the time series bilateral plantar force data to specific deep learning predictive models for running fatigue gait, the results showed that both CNN and ConvLSTM models have good performance in predicting fatigue gait automatically. As expected, the ConvLSTM model (85.9%, 88.9%, 83.3%, 85.9%, and 89.2%, respectively) has better accuracy compared to the CNN model (85.6%, 73.8%, 82.7%, 80.3%, and 75.4%) in all five tests, suggesting that ConvLSTM performs better for multi-feature data with simultaneous spatiotemporal dependence [27]. Traditional time series biomechanical datasets are all characterized by high dimensionality, high variability, time dependence, and nonlinearity [38]. Therefore, with the promising findings from this study as a foundation, future research suggests applying the Convlstm model to other analyses, such as marker trajectories, ground reaction forces, myoelectric signals, and other prediction and classification needs. In addition, as shown in Table 3, the specificity of the ConvLSTM model was also higher than that of the CNN, indicating that it could detect fatigue gait better, while the performance of sensitivity was consistent in both models, indicating that both models were equally effective in predicting non-fatigue gait.

There are four limitations to this study. Possible differences in plantar pressure distribution patterns between overground conditions and treadmill conditions were reported by García-Pérez et al. [39]. In this study, the running-induced fatigue protocol was carried out on a treadmill. Thus, further investigation is needed to support our findings in overground conditions. The runners were evaluated under barefoot conditions, potentially overlooking the impact of footwear on post-fatigue running posture [40]. Additionally, we selected only two deep learning models (CNN and ConvLSTM) for data training based on data features, and more comparisons of classifiers (such as deep neural network) for plantar pressure feature discovery should be developed in future studies. At the end, only the pedobarographic data of amateur male runners was included in this study; whether the model developed in this study applies to female or elite runners should be verified in future studies.

5. Conclusions

A running-induced fatigue protocol caused different changes in the distribution of plantar force on the dominant and nondominant limbs. These changes may be part of the underlying mechanism of unilateral limb overuse injuries. Future discussions of lower limb lesions or running-related injuries should take this into account. Furthermore, the ConvLSTM model showed high performance (acc = 0.867) in detecting fatigue gait, and it outperformed the CNN model (0.800). This will broaden the possibilities for future research on running-related gait biomechanical feature recognition and enhance the development of fatigue monitoring tools.

Data Availability

Data are available upon request due to privacy restrictions. The data presented in this study may be available upon request from the corresponding author and with the authorization of funding origination.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to acknowledge the contributions of all participants and staff involved in this study. This study was supported by the Research Academy of Medicine Combining Sports, Ningbo (No. 2023001), the Project of NINGBO Leading Medical & Health Discipline (No. 2022-F15, No. 2022-F22), Ningbo Natural Science Foundation (20221JCGY010532, 20221JCGY010607), Public Welfare Science & Technology Project of Ningbo, China (2021S134), Zhejiang Provincial Natural Science Foundation of China for Distinguished Young Scholars (LR22A020002), Zhejiang Provincial Natural Science Foundation (LTGY23H040003), Zhejiang Provincial Key Research and Development Program of China (2021C03130), the János Bolyai Research Scholarship of the Hungarian Academy of Sciences (BO/00047/21/6), and K. C. Wong Magna Fund in Ningbo University.