Machine Learning-Based Risk Identification and Precision Exercise for Common Running Injuries

Overview

This study aims to identify biomechanical risk factors associated with common running-related injuries in recreational runners using machine learning analysis. It also aims to evaluate whether a precision exercise intervention based on these risk factors can improve injury-related biomechanical and kinematic outcomes. The main questions it aims to answer are: Which biomechanical features identified by machine learning are associated with the occurrence of common running-related injuries, including medial tibial stress syndrome (MTSS), patellofemoral pain (PFP), and chronic Achilles tendinopathy? Whether a precision exercise intervention based on these risk factors can improve injury-related biomechanical and kinematic characteristics？ Participants will: Undergo baseline biomechanical assessment during running, including motion capture, ground reaction force, and surface electromyography Be prospectively followed for 6 months to record the occurrence of running-related injuries Be classified into injury groups based on diagnosis, including MTSS, PFP, and chronic Achilles tendinopathy Following the completion of the follow-up period, participants will be allocated to either a precision multidimensional intervention group, a precision exercise intervention group, or a control group. The precision multidimensional intervention group will receive patient education in addition to an 8-week intervention program, including stretching, strength training, and real-time movement feedback. The precision exercise intervention group will receive the same 8-week intervention program, consisting of stretching, strength training, and real-time movement feedback. The control group will maintain their habitual physical activity patterns without receiving any additional intervention. Participants will be reassessed after the intervention and at the 3-month follow-up using the same biomechanical testing protocol.

1. Participants This prospective study recruited recreational runners from East China Normal University between November and December 2024 through poster advertisements and WeChat-based recruitment. Eligible participants were male recreational runners aged 18 to 30 years with a weekly running distance greater than 10 km and no history of lower-limb injury within the past three months. Exclusion criteria included: (1) the presence of other musculoskeletal or tissue injuries, such as arthritis, muscle injuries, or meniscal injuries; (2) a history of lower-limb surgery; and (3) regular participation in contact or high-impact sports, in order to minimize potential confounding effects from non-running-related musculoskeletal injuries. A priori sample size estimation was performed using GPower software (GPower V3.1, Heinrich-Heine-Universität Düsseldorf, Germany). With a significance level of α = 0.05, statistical power of 0.80, and an effect size of 1.7, the minimum required sample size was estimated to be 170 participants. Considering an anticipated attrition rate of 20%, a total of 204 participants were ultimately recruited. All participants were subsequently followed prospectively for 6 months and classified into four groups based on injury status: injury-free, chronic Achilles tendinopathy, MTSS, and PFP. Diagnostic criteria for chronic Achilles tendinopathy included: (1) insidious Achilles tendon pain aggravated by weight-bearing activities and worse in the morning and/or at the onset of activity; (2) pain and swelling located 2-6 cm proximal to the Achilles tendon insertion; (3) ultrasonographic evidence of tendon thickening \>7 mm or \>20% compared with the contralateral side; and (4) symptoms persisting for more than 3 months. MTSS was defined as: (1) diffuse medial tibial pain during exercise with a palpation tenderness area of at least 5 cm; (2) symptom duration of at least 3 weeks; and (3) absence of stress fracture or other fractures on radiographic examination. PFP was defined as: (1) pain around or behind the patella during at least two activities such as running, squatting, stair ascent/descent, prolonged sitting, jumping, or resisted knee extension; (2) pain not attributable to direct trauma and persisting for at least 3 months; and (3) a minimum pain intensity of ≥3 on a 10-cm Visual Analog Scale (VAS), where 0 indicates no pain and 10 indicates maximal pain. During the intervention phase, participants were randomly allocated to the Precision Multidimensional Intervention Group, the Precision Exercise Intervention Group, or the control group using simple randomization, with the allocation sequence generated by a computerized random number generator. Allocation concealment was achieved using sealed, opaque envelopes, which were opened sequentially according to enrollment order. The study employed a double-blind design during the intervention period. Unblinding was performed after completion of data analysis. However, in emergency situations-such as serious adverse events where the relationship with the intervention could not be determined and immediate clinical management decisions were required-early unblinding was permitted. 2. Data Collection and Testing Protocol Prior to the prospective follow-up, all participants underwent a comprehensive biomechanical assessment, including condition-specific testing under different temperature environments (e.g., low and high temperature conditions). The assessment protocol was as follows. Participants first performed a 3-minute static stretching session and then changed into standardized laboratory footwear, tight-fitting tops, and shorts. A total of 39 reflective markers (14 mm in diameter) were placed on anatomical landmarks according to the Plug-in-Gait full-body model, including the left and right anterior head, C7, T10, left and right anterior superior iliac spine, left and right posterior superior iliac spine, thigh, lateral knee, shank, ankle, toe, and heel. Static calibration trials were then collected, during which participants stood still in the center of the capture volume with feet shoulder-width apart and arms abducted at the elbows with palms facing downward. For the dynamic trials, participants were instructed to run at a self-selected speed along the laboratory runway, maintaining their habitual running pattern. Three-dimensional kinematic data were collected using a 12-camera Vicon V2.2 motion capture system (Vicon, UK) at a sampling frequency of 200 Hz. Two embedded force plates were installed in the middle of the runway. Participants were instructed to adjust their stride length to ensure that the affected limb accurately contacted the force plate during running, allowing for the acquisition of kinetic data at 1000 Hz. Surface electromyography (sEMG) data were collected using a 32-channel Noraxon system (Noraxon, USA) at 1000 Hz to record muscle activity of the rectus femoris, biceps femoris, semimembranosus, semitendinosus, tibialis anterior, and lateral gastrocnemius, which were used for validation of the OpenSim musculoskeletal model. A total of five successful trials were collected for each participant, with a 1-minute rest interval between trials. 3. Data Processing Data preprocessing was performed using Vicon Nexus software (Vicon, Oxford, UK), with a primary focus on biomechanical analysis during the stance phase of running. Initial contact was defined when ground reaction force (GRF) exceeded 20 N, and toe-off was identified when GRF fell below 20 N, thereby enabling precise determination of the stance phase. Marker trajectories (200 Hz) and ground reaction force (GRF) data (1000 Hz) were filtered using a fourth-order, zero-lag Butterworth low-pass filter at a cutoff frequency of 20 Hz. Subsequently, transformation matrices were constructed, joint centers were estimated, and segmental local coordinate systems were defined to calculate three-dimensional kinematics of the trunk, hip, knee, and ankle joints. Joint moments were derived using inverse dynamics. Surface electromyography (sEMG) data were processed in Noraxon software (Noraxon MR3, USA) using the root mean square (RMS) method. Muscle activation levels during the stance phase were normalized to maximal voluntary contraction (MVC), yielding standardized activation values ranging from 0 (no activation) to 1 (maximal activation). All derived datasets were used for validation of the OpenSim musculoskeletal model. 4. OpenSim Musculoskeletal Modeling Kinematic data were exported as .trc files, and kinetic data were converted into .mot format files prior to import into the OpenSim platform for musculoskeletal modeling. The modeling pipeline included model scaling, inverse kinematics, residual reduction analysis, and muscle control estimation. Finally, model validity was assessed by comparing OpenSim-derived outputs with experimentally measured data. 4.1 Model Scaling A generic musculoskeletal model was first constructed using participant-specific anthropometric data, including height, body mass, and individualized muscle characteristics. The model was subsequently scaled based on experimentally acquired marker data to adjust segment lengths and masses. 4.2 Inverse Kinematics Inverse kinematics was employed to identify the optimal correspondence between the OpenSim model and experimental laboratory data. A weighted least-squares optimization approach was used to minimize discrepancies between experimentally measured three-dimensional marker coordinates, coordinate systems, and model-derived outputs, thereby reducing overall error to a minimum. 4.3 Residual Force Calculation Due to discrepancies between experimentally measured ground reaction forces and those computed using Newtonian mechanics, residual forces (F\_residual) arise, defined as: F\_residual=F\_exp-m×a In OpenSim, such residuals are reduced through trajectory optimization and segment mass adjustment. Residual forces are generally considered acceptable within a range of 0-10 N, while residual moments are considered acceptable within 0-50 N·m. 4.4 Muscle Control Estimation Lower-limb muscle force was estimated using a muscle control framework. Subsequently, an objective function was defined, typically comprising dynamic residuals, marker tracking errors, and muscle activation costs. Optimization algorithms (e.g., gradient descent and interior-point methods) were then employed to iteratively update muscle activation levels in order to minimize the objective function. The optimized muscle force parameters were ultimately obtained. 4.5 Model Validation To evaluate simulation accuracy, the validation approach proposed by Błażkiewicz et al. was adopted, in which OpenSim-derived kinematics, kinetics, and muscle activation outputs were compared with experimentally measured data. The model was considered acceptable if the peak sagittal plane angles of the hip, knee, and ankle differed by less than 5°, pelvic residual forces and moments were below 20 N and 75 Nm, respectively, and the muscle activation profiles showed satisfactory agreement between simulation and experimental data. Following validation, 28 lower-limb muscles were included in the machine learning model, comprising: ipsilateral internal oblique muscle force, ipsilateral external oblique muscle force, gluteus maximus, gluteus medius, gluteus minimus, rectus femoris, vastus medialis, vastus lateralis, vastus intermedius, tensor fasciae latae, adductor longus, adductor magnus, adductor brevis, gracilis, semimembranosus, semitendinosus, biceps femoris, tibialis anterior, peroneus longus, peroneus brevis, medial gastrocnemius, lateral gastrocnemius, tibialis posterior, soleus, flexor hallucis longus, flexor digitorum longus, extensor digitorum longus, and extensor hallucis longus. 5. Machine Learning 5.1 Data Selection and Preprocessing Missing data were first handled by excluding features with more than 30% missing values. For variables with less than 30% missingness, imputation was performed using the k-nearest neighbor (KNN) algorithm. The dataset was then randomly split into training and testing sets using stratified sampling at a ratio of 8:2. Multicollinearity among selected features was subsequently assessed by calculating the variance inflation factor (VIF) and correlation coefficients (R values). Features with VIF \> 10 or \|R\| \> 0.8 were considered to exhibit multicollinearity and were removed from further analysis. Given the potential class imbalance in the dataset, the synthetic minority over-sampling technique (SMOTE) was applied to the training set to augment minority class samples. SMOTE generates synthetic samples for the minority class based on k-nearest neighbor interpolation 18, thereby ensuring balanced class distributions in the training dataset. All preprocessing procedures were implemented using the "imblearn" package in Python. 5.2 Machine Learning Model Development Three machine learning algorithms-XGBoost, random forest, and support vector machine (SVM)-were employed to analyze and predict biomechanical risk factors associated with common running-related injuries, including chronic Achilles tendinopathy, MTSS, and PFP. Both XGBoost and random forest are ensemble learning methods based on decision trees and are widely used for classification and regression tasks, with the capability to handle high-dimensional datasets containing millions of samples and features. In contrast, SVM can effectively address nonlinear separable problems through kernel functions such as linear, polynomial, and radial basis function (RBF) kernels. All models were implemented using the scikit-learn package in Python 3.8. Model performance was evaluated using cross-validation, in which the training dataset was partitioned into multiple subsets; the model was trained on one subset and validated on another. This process was repeated iteratively, and the results were averaged to assess model generalizability on unseen data, while also evaluating potential underfitting and overfitting. Hyperparameter tuning was performed for XGBoost, random forest, and SVM to identify optimal parameter configurations. Specifically, the hyperparameter space was discretized into a grid, and grid search was conducted to identify the combination that maximized cross-validation performance. The optimal hyperparameters were then used to retrain and optimize all models. 5.3 Model Validation and Evaluation The predictive performance of each machine learning model was evaluated using the following metrics: accuracy, recall, precision, and F1-score. In addition, receiver operating characteristic (ROC) curves were plotted for each model, with the true positive rate on the y-axis and the false positive rate on the x-axis, and the area under the curve (AUC) was calculated to quantify overall discriminative ability. The optimal machine learning model was selected based on the above evaluation metrics and subsequently used to identify biomechanical risk factors associated with common running-related injuries. The formulas for accuracy, recall, precision, and F1-score are defined as follows: Accuracy = (TP + TN) / (TP + TN + FP + FN) Recall = TP / (TP + FN) Precision = TP / (TP + FP) F1-score = (2 × precision × recall) / (precision + recall) where TP represents true positives, FP false positives, TN true negatives, and FN false negatives. 5.4 Visualization of Model Outputs In this study, SHapley Additive exPlanations (SHAP) values were computed for each feature within the four-class machine learning model. SHAP, originally proposed by Lundberg and Lee 19, is a model-agnostic interpretability framework whose robustness and explanatory power have been extensively validated 20. SHAP enables both local (instance-level) and global (population-level) interpretation of model predictions. Compared with other model interpretation methods, SHAP is grounded in solid theoretical foundations derived from cooperative game theory. Accordingly, SHAP analysis was employed to interpret the outputs of the machine learning models, enabling the identification of biomechanical risk factors associated with common running-related injuries, as well as quantification of the relative importance of these risk factors. 6. Construction of the Precision Exercise Intervention Model 6.1 Precision Multidimensional Intervention Group Based on the above machine learning results, a targeted intervention protocol was developed. Participants in the Precision Multidimensional Intervention Group received an 8-week precision intervention program designed according to biomechanical risk factors identified by machine learning. The intervention comprised stretching exercises, strength training, and real-time movement feedback training (see Tables 1 and 2 for the strength training components). It should be noted that Tables 1 and 2 are provided as illustrative examples only, and the specific exercise prescription was individually tailored based on the risk factors identified through machine learning analysis. The real-time feedback component was implemented using reflective markers in combination with a high-speed infrared motion capture system. This setup enabled continuous extraction of kinematic parameters, which were provided to participants in real time. This feedback loop allowed ongoing correction of movement patterns, promoting a gradual transition toward a running pattern associated with reduced injury risk. In addition, participants in the Precision Multidimensional Intervention Group received three individualized, one-on-one educational sessions over the 8-week intervention period, each lasting approximately 30-45 minutes. The first session focused on characterizing the participant's pain profile, including pain location, the contexts in which pain was most severe during daily activities, their understanding of pain etiology, and expectations regarding treatment outcomes. Based on this information, the subsequent two sessions were tailored by the therapist to individual characteristics. The educational content included: (1) an explanation of the mechanisms underlying the specific injury type, including the influence of ankle and proximal joint kinematics and muscle activation patterns on tibial loading; (2) pain management strategies during running, such as avoiding excessive ankle plantarflexion and hip external rotation while maintaining appropriate lower-limb alignment and shock attenuation strategies; (3) training load regulation strategies based on cadence control and load management principles, for example reducing activity or load when pain exceeds 5/10; and (4) addressing patient-specific questions. Education was delivered using a combination of verbal explanation and visual illustrations to enhance comprehension. 6.2 Precision Exercise Intervention Group Participants in the Precision Exercise Intervention Group received an 8-week exercise-based intervention. The stretching and strength training components were identical to those implemented in the Precision Multidimensional Intervention Group; however, the structured patient education module was not included. Accordingly, this group received a standardized exercise prescription based solely on the identified risk factors derived from machine learning analysis. 6.3 Control group The control group maintained their habitual running activities and daily physical activity levels throughout the study period and did not receive any additional intervention. 6.4 Follow-up Assessment A follow-up biomechanical assessment was conducted 3 months after the completion of the intervention. The follow-up testing protocol was identical to that used at baseline and immediately post-intervention, including all biomechanical measurements of kinematics, kinetics, muscle activation, and OpenSim-derived musculoskeletal parameters. All participants underwent testing under the same experimental conditions to ensure consistency and comparability across time points. 7. Statistical analysis Statistical analyses were performed using SPSS 22.0. For normally distributed data, a two-way repeated-measures analysis of variance (ANOVA) was used to examine differences across groups (intervention vs control) and time points (pre-intervention, post-intervention, and follow-up) for chronic Achilles tendinopathy, MTSS, and PFP outcomes. Data are presented as mean ± standard deviation (M ± SD). In the presence of significant interaction effects, between-group differences were further examined using independent-samples t-tests, while within-group pre-post comparisons were conducted using paired-samples t-tests. Main effects were adjusted using the least significant difference (LSD) method. All statistical tests were two-sided, and a significance level of P \< 0.05 was considered statistically significant.