We enrolled 65 adults with chronic stroke and 20 age-matched healthy controls. Stroke participants were allocated to a model-development cohort (n = 47) and an independent test cohort (n = 18), with the test cohort held out throughout model development and tuning. All participants were recruited under identical eligibility criteria. Inclusion criteria were: first-ever unilateral stroke confirmed by neuroimaging with persistent motor impairment of the paretic upper limb; chronic stage (≥ 6 months post-stroke); sufficient voluntary movement to complete the hand-to-mouth task, operationalized as active elbow flexion ≥ 30°; ability to sit unsupported and follow verbal commands; and adequate cognition (MMSE ≥ 24). Exclusion criteria were: recurrent or bilateral stroke; marked upper-limb spasticity (Modified Ashworth Scale > 3); musculoskeletal deformity or pain limiting upper-limb motion; severe neglect, aphasia, or apraxia compromising task performance; and cardiovascular or orthopedic conditions judged unsafe or likely to interfere with standardized testing.

A standardized hand-to-mouth (HTM) task (Cai et al., 2021), known for good test–retest reliability, was used as the primary functional assessment to measure upper limb motor ability during drinking, eating, and face-approach activities. The task required participants to reach from a starting position on a table to touch a marker, then lift their hand to touch their lips with the thumb, return to touch the starting marker, and finally return to the initial position. The marker was positioned 30 cm from the table edge at the body midline, calibrated to approximately 80% of individual arm length. Participants sat on an adjustable-height chair with their back against the chair back, but sitting posture was not rigidly constrained to allow necessary compensatory movements. The initial posture required the upper arm hanging naturally, elbow flexed approximately 90°, and the tested hand resting palm-down on the table. All movements were performed at a self-selected pace to better reflect natural ADL performance and avoid inducing abnormal compensatory patterns that may arise from speed-focused instructions. Only the paretic upper limb was tested in the stroke group; the dominant hand was tested in the control group. Each participant completed five task trials; the average of the middle three trials was used for final analysis. Adequate rest was provided between trials to avoid muscle fatigue.

All movements were completed during a single laboratory visit, with synchronous acquisition of surface electromyography (sEMG) and inertial measurement unit (IMU) data. sEMG signals were recorded at 2000 Hz using a wearable system (Noraxon USA, Inc., Scottsdale, AZ, United States) from seven upper limb muscles: Brachioradialis (BR), Biceps Brachii (BB), Triceps Brachii Lateral Head (TRL), Anterior Deltoid (AD), Middle Deltoid (MD), Posterior Deltoid (PD), and Pectoralis Major (PM). Bipolar Ag-AgCl circular surface electrodes, spaced 2 cm apart, were placed longitudinally over the muscle belly center according to recommendations from the International Society of Electrophysiology and Kinesiology (ISEK), following shaving and alcohol cleansing. IMU modules recorded data at 400 Hz and were fixed to the sternum, upper arm, and wrist. Data were wirelessly transmitted via Bluetooth to a host computer for synchronized storage and processing. The placement of sensors and electrodes is illustrated in Figure 1.

The analysis pipeline commenced with preprocessing of the raw sensor data. Joint kinematics were derived from the raw IMU data (angular velocity and acceleration) using established sensor fusion and calibration procedures (Nazarahari and Rouhani, 2021). Concurrently, raw sEMG signals were band-pass filtered, full-wave rectified, and low-pass filtered to create linear envelopes; these envelopes were then normalized to the peak value recorded during the task cycles (Lloyd and Besier, 2003). These processed kinematic and sEMG data served as the primary inputs for the subsequent musculoskeletal modeling pipeline.

Prior to kinetic analysis, the open-source OpenSim (v4.3) upper limb musculoskeletal model (Seth et al., 2019) was scaled to match each participant’s anthropometry. This critical step utilized participant-specific measurements (e.g., height, weight, and segment lengths derived from IMU positions) to adjust the model’s body segment parameters and muscle-tendon paths, thereby mitigating potential anthropometric bias. Once scaled, kinematic data from IMUs and preprocessed sEMG signals were input into the model. The analysis was performed using the Calibrated EMG-Informed Neuromuscular Modeling toolbox (CEINMS) (CEINMS: a toolbox to investigate the influence of different neural control solutions on the prediction of muscle excitation and joint moments during dynamic motor tasks Pizzolato et al., 2015). The pipeline commenced with the OpenSim Scaling tool, followed by Inverse Kinematics (IK) to obtain joint angles, and then Inverse Dynamics (ID) to calculate net joint torques. Crucially, to enhance the subject-specific physiological accuracy, the CEINMS calibration procedure was performed. This standard step iteratively optimized key musculoskeletal parameters (e.g., tendon slack length) to minimize the error between the EMG-driven model-predicted joint torques and the ID-calculated net joint torques. Subsequently, this fully calibrated, subject-specific model was implemented in an EMG-informed mode: recorded EMG signals directly informed the activations of the corresponding measured muscles, while activations for unmeasured muscles were estimated through optimization (using the standard static optimization algorithm in the CEINMS toolbox, which aims to minimize the sum of squared muscle activations). This EMG-driven approach, which incorporates muscle activation-force generation dynamics models, was employed to estimate force output time-series for major muscles during each task. The entire kinetic computation pipeline followed the standard, validated CEINMS workflow (CEINMS: a toolbox to investigate the influence of different neural control solutions on the prediction of muscle excitation and joint moments during dynamic motor tasks Pizzolato et al., 2015; Seth et al., 2019; Tahmid et al., 2022), ensuring scientific rigor and reproducibility.

The HTM task was segmented into four consecutive subphases: Reach to object phase (Phase 1), Transfer to mouth phase (Phase 2), Transfer to object phase (Phase 3) and Return to start phase (Phase 4). Phase segmentation was based on 3D spatial trajectory and velocity characteristics recorded by the wrist IMU. Phase segmentation was performed by identifying key kinematic landmarks within the 3D velocity profile recorded by the wrist IMU (e.g., movement onset, velocity peaks, and movement offset), following established protocols for HTM task analysis (Schwarz et al., 2022; Huang et al., 2024). This algorithmic segmentation was then visually verified by trained researchers to ensure precise temporal delineation of the movement sequence, providing a clear structure for subsequent phase-specific feature extraction.

Based on the IMU- and sEMG-driven musculoskeletal model, key kinetic features were systematically calculated and extracted during each phase of the movement, providing structured input data for subsequent ML-based multi-phase sequence regression analysis and FMA-UL scores prediction. Joint kinematic features were extracted based on the upper limb degrees of freedom defined in the OpenSim model, where elb denotes elbow flexion/extension, elv denotes shoulder elevation/depression, rot denotes shoulder internal/external rotation, and fle denotes shoulder flexion/extension. For each HTM task phase, five categories of parameters were calculated:

Mechanical work (W) (Seth et al., 2019; Wright et al., 2020) performed by joints and muscles was calculated by integrating muscle power time-series, quantifying an individual’s force output contribution during task execution.

The average value of the middle three trials for each participant was used for statistical analysis. Feature data were first tested for normality. For normally distributed features, independent samples t-tests were used to compare differences between the healthy and stroke groups. For non-normally distributed features, the Mann–Whitney U test was used. Descriptive statistics: normally distributed features are reported as mean ± standard deviation; non-normally distributed features are reported as median (interquartile range, IQR). All feature labels and abbreviations (e.g., TD_2, S_elb_1) used in this study are defined in the note of Table 1. To account for multiple comparisons, all p values were adjusted using the Benjamini-Hochberg False Discovery Rate (FDR) procedure. Statistical significance was set at an FDR-adjusted p < 0.05. To quantify the magnitude of group differences, effect sizes were calculated using Cohen’s d. Statistical analyses were performed using SPSS 24.0 (IBM, United States).

To systematically evaluate the predictive capability of the models and ensure their generalizability, a two-stage validation strategy was employed: internal validation (based on the model derivation cohort, n = 47) and validation on an independent test cohort (n = 18).

Within the derivation cohort, eight typical supervised regression algorithms were included for internal comparison: Linear Regression, Ridge Regression, Lasso Regression, Elastic Net Regression, Decision Tree, Extra Trees Regressor, Random Forest Regressor, and Gradient Boosting Regressor. To identify the optimal model and prevent data leakage, all training, tuning, and evaluation procedures were conducted within a rigorous 5-fold cross-validation (CV) pipeline. Crucially, within each of the 5 folds, all preprocessing steps—specifically data standardization (using StandardScaler to achieve zero mean and unit variance) and any feature selection (e.g., RFE)—were fitted only on the training partition (4 folds). The fitted pipeline was then applied to transform both the training partition and the held-out validation partition. Hyperparameter tuning for each algorithm (e.g., the regularization parameter α for Lasso) was also performed within this nested CV structure, using the coefficient of determination (R²) as the optimization metric. The average performance metrics [R², Mean Squared Error (MSE), Mean Absolute Error (MAE)] derived from the held-out validation folds are reported as the model’s internal validation performance.

Based on the internal validation results, the best-performing algorithm (Lasso Regression), which also offered the best balance of predictive accuracy and model interpretability, was identified. Subsequently, a final model of this algorithm was refit using the entire derivation cohort (n = 47), with its hyperparameters fixed to the optimal values identified during the cross-validation. Finally, this trained model was evaluated on the completely independent test cohort (n = 18). The R², MSE, and MAE calculated on this dataset are reported as the model’s independent test performance, providing an unbiased assessment of its generalizability to new patient data.

To enhance model interpretability, feature importance was extracted and ranked based on models with intrinsic feature weights (e.g., linear models and tree-based models). Furthermore, SHAP (SHapley Additive exPlanations) was applied to the representative model to generate both model-level (SHAP summary plots) and instance-level (SHAP force plots) explanations, thereby revealing the direction and magnitude of contributions from key features to the predictions.

A total of 65 stroke patients and 20 healthy controls were included in this study. The stroke patients were divided into a model development cohort (n = 47) and an independent test cohort (n = 18). As summarized in Table 2, there were no significant differences in demographic or key clinical characteristics, including age, sex, assessed side, type of brain injury, time since injury, and baseline FMA-UL scores, between the two cohorts (all p > 0.05).

To identify phase-specific kinetic impairments that constitute the basis for subsequent machine learning prediction, we first conducted a comprehensive group comparison. The analysis revealed a constellation of phase-specific neuromuscular alterations (Table 1), which delineate the pathophysiology of upper limb dysfunction and provide the foundational feature pool for FMA-UL score prediction.

During the Reach to object phase (Phase 1), dominated by shoulder elevation and elbow extension, patients exhibited impaired motor output stability. This was characterized by significantly decreased smoothness (all p ≤ 0.036) of joint moments (S_elb_1, S_elv_1) and multiple muscles (S_PD_1, S_MD_1, S_PM_1). These differences showed moderate to large effect sizes (e.g., S_BRA_1, d = −0.817; S_PD_1, d = −0.511). Elevated co-contraction of shoulder synergistic pairs (CCI_MD_PM_1; p = 0.008, d = 0.908) further indicated abnormal stiffness.

In the Transfer to mouth phase (Phase 2), which involves elbow flexion and shoulder elevation, patients demonstrated diminished propulsive power and disrupted coordination. This was evidenced by a 34% reduction in elbow flexion concentric work (W_elb_2, p = 0.017, d = −0.931) and a 38% decrease in biceps mechanical work (W_BB_2, p = 0.017, d = −0.805). Coordination between elbow and shoulder motion was impaired in patients (IC_elb_elv_2, p = 0.015, d = −0.502). This deficit was accompanied by markedly reduced smoothness across the elbow joint (S_elb_2) and multiple muscles (S_fle_2, S_rot_2, S_BRA_2, S_TRI_2, S_AD_2, S_PD_2; all p ≤ 0.034), with several showing large effects (e.g., S_TRI_2, d = −1.330; S_BRA_2, d = −0.942). Patients also showed increased compensatory trunk displacement (TD_2, p = 0.036, d = 0.680) and elevated co-contraction across all measured muscle pairs (all CCI pairs, p ≤ 0.012), all with large effect sizes (d ≥ 0.854).

The Transfer to object phase (Phase 3), which requires controlled elbow extension, revealed deficits in eccentric control. Patients showed significantly altered elbow eccentric work (W_elb_3, p = 0.017, d = 0.866) and altered biceps and triceps work (W_BB_3, p = 0.023; W_TRI_3, p = 0.036). Elbow-shoulder coordination was significantly degraded (IC_elb_elv_3, p = 0.008, d = −0.625). Smoothness of the elbow moment (S_elb_3) and brachioradialis force (S_BRA_3) was also compromised (both p ≤ 0.037). A pronounced increase in antagonist co-contraction across all muscle pairs (all p ≤ 0.004) and increased trunk displacement (TD_3, p = 0.012) were observed, with large effect sizes (d ≥ 0.682).

Finally, in the Return to start phase (Phase 4), dominated by shoulder abduction, the results highlighted profound compensatory mechanisms. The most salient finding was a large increase in trunk displacement (TD_4, p = 0.015, d = 0.612). This was associated with altered elbow-shoulder coordination (IC_elb_elv_4, p = 0.033) and decreased smoothness across the shoulder joint (S_elv_4) and multiple muscles (S_fle_4, S_rot_4, S_BRA_4, S_TRI_4, S_MD_4, S_PD_4; all p ≤ 0.046). Concurrently, all muscle pairs showed elevated co-contraction (all p ≤ 0.008), indicating inefficient control and stiffness, with all differences showing large effect sizes (d ≥ 0.911).

To establish the optimal FMA-UL score prediction model, we systematically compared eight supervised learning regression algorithms. The performance of all models was evaluated on the development cohort using 5-fold cross-validation (internal validation) and on the independent hold-out test cohort (n = 18), with detailed results presented in Table 3. The evaluation results indicated that regularized linear models (Lasso, Ridge, ElasticNet) generally outperformed tree-based models (e.g., Random Forest, ExtraTrees) and standard Linear Regression for this task. Among all tested algorithms, the Lasso Regression model demonstrated the best overall performance in both internal validation (R² = 0.932) and independent test.

On the more challenging independent test set, the Lasso model achieved the highest coefficient of determination (R² = 0.881) and the lowest Mean Squared Error (MSE = 1.802) and Mean Absolute Error (MAE = 0.954). This indicates the model possesses strong generalization capability. Figure 2 illustrates the relationship between the predicted values from the Lasso model and the actual FMA-UL scores on the independent test set. Data points are tightly clustered around the ideal fit line (y = x), visually confirming the model’s high predictive accuracy. Furthermore, the color of the scatter points represents the absolute error, with most points skewed toward blue (low error), further demonstrating that the model’s predictions are reliable at the individual level.

To understand how the Lasso model makes predictions, we first analyzed its internal structure. Figure 3 displays the 10 features with the largest absolute Lasso regression coefficients (weights) in the final model. As the plot indicates, TD_2 (trunk displacement in phase 2) has the largest coefficient, underscoring its prominent role in the model’s structural composition, followed by IC_elb_elv_3 (elbow–shoulder coordination in phase 3) and TD_3 (trunk displacement in phase 3).

However, while these coefficients define the model’s structure, their rank order can be sensitive to collinearity among correlated features, which is common in biomechanical data. Therefore, to move from reporting the model’s structure (Figure 3) to a more robust analysis of each feature’s practical impact on the prediction, we employed SHAP (SHapley Additive Explanations).

Figure 4 is a SHAP summary plot, which ranks features based on their mean absolute SHAP value and illustrates the impact of each feature’s value (High = Red, Low = Blue) on the model prediction (positive SHAP values push the prediction higher; negative SHAP values push it lower). From the SHAP summary plot, it can be observed that: TD_2 was re-confirmed as the most important feature. Its low values (blue dots) are primarily distributed in the positive SHAP value region, indicating that smaller trunk displacement in phase 2 tends to increase the model’s prediction for the FMA-UL score. IC_elb_elv_3 (phase 3 coordination) and TD_4 (phase 4 trunk displacement) also showed significant impacts. For instance, high values (red dots) for IC_elb_elv_3 are mostly on the positive side, while low values (blue dots) are more skewed to the negative side. TD_4 also exhibited a trend where low values skewed positive. Additionally, high values (red dots) for W_PD_4 (phase 4 posterior deltoid work) and S_TRI_2 (phase 2 triceps smoothness) were also mainly distributed in the positive SHAP value region, suggesting that increases in these kinetic and EMG-related metrics positively contribute to the model’s prediction of functional recovery.

In addition to global importance, we also used SHAP force plots to explain how the model makes predictions for specific individual patients. Figure 5 presents two representative individual samples.

In these plots, the f ( x ) value is the model’s final predicted FMA-UL score for that patient. The base value (approx. 24.2) represents the model’s average predicted score across the entire training dataset. Red arrows represent features that “push” the prediction higher than the base value (i.e., factors increasing the FMA-UL score). Blue arrows represent features that “pull” the prediction lower than the base value (i.e., factors decreasing the FMA-UL score). The length of each arrow represents the magnitude of that feature’s impact (i.e., the magnitude of its SHAP value).

Figure 5A shows a sample with a high predicted score [ f ( x ) = 27.26]. For this patient, multiple features (red) acted in concert to push the predicted score above the base value, with the most impactful factors including CCI_AD_PD_4, CCI_AD_PD_2, and S_elb_2. The only significant negative factor (blue) was CCI_AD_PD_3. Figure 5B, in contrast, shows a sample with a low predicted score [ f ( x ) = 20.52]. This prediction is significantly lower than the base value. The main driving factors (blue) for this lower prediction include features such as IC_elb_elv_2, W_TRI_1, and S_TRI_2. Meanwhile, features like W_elb_3 and W_TRI_4 (red) prevented the score from being even lower to some extent. These individualized explanations are crucial for clinically understanding the specific functional impairment patterns of a given patient.