Predicting upper limb motor recovery in subacute stroke patients via fNIRS-measured cerebral functional responses induced by robotic training

Abstract

Background

Neural activation induced by upper extremity robot-assisted training (UE-RAT) helps characterize adaptive changes in the brains of poststroke patients, revealing differences in recovery potential among patients. However, it remains unclear whether these task-related neural activities can effectively predict rehabilitation outcomes. In this study, we utilized functional near-infrared spectroscopy (fNIRS) to measure participants’ neural activity profiles during resting and UE-RAT tasks and developed models via machine learning to verify whether task-related functional brain responses can predict the recovery of upper limb motor function.

Methods

Cortical activation and brain network functional connectivity (FC) in brain regions such as the superior frontal cortex, premotor cortex, and primary motor cortex were measured using fNIRS in 82 subacute stroke patients in the resting state and during UE-RAT. The Fugl-Meyer Upper Extremity Assessment Scale (FMA-UE) was chosen as the index for assessing upper extremity motor function, and clinical information such as demographic and neurophysiological data was also collected. Robust features were screened in 100 randomly divided training sets using the least absolute shrinkage and selection operator (LASSO) method. Based on the selected robust features, machine learning algorithms were used to develop clinical models, fNIRS models, and combined models that integrated both clinical and fNIRS features. Finally, Shapley Additive Explanations (SHAP) was applied to interpret the prediction process and analyze key predictive factors.

Results

Compared to the resting state, task-related FC is a more robust feature for modeling, with screening frequencies above 90%. The combined models built using artificial neural networks (ANNs) and support vector machines (SVMs) significantly outperformed the other algorithms, with an average AUC of 0.861 (± 0.087) for the ANN and an average correlation coefficient (r) of 0.860 (± 0.069) for the SVM. Furthermore, predictive factor analysis of the models revealed that FC measured during tasks is the most important factor for predicting upper limb motor function.

Conclusion

This study confirmed that UE-RAT-induced FC can serve as an important predictor of rehabilitation, especially when combined with clinical information, further enhancing the accuracy of model predictions. These findings provide new insights for the early prediction of patients’ recovery potential, which may contribute to personalized rehabilitation decisions.

Background

More than two-thirds of stroke survivors have upper limb dyskinesia, which seriously affects their quality of life and social interactions [1]. Improving upper limb motor function has therefore become a key goal of rehabilitation therapy [2]. However, the heterogeneity among patients makes it challenging to assess rehabilitation potential in the early stages and may hinder the selection of optimal rehabilitation strategies. Therefore, it is particularly critical to identify predictors of rehabilitation, which are not only central to assessing rehabilitation potential, but also provide the basis for individualized treatment plans.

In recent years, machine learning algorithms such as support vector machines (SVMs) have been used to identify predictors of the potential for rehabilitation, such as time since stroke and initial FMA-UE scores [3,4,5]. These indicators provide a basis for predicting rehabilitation, but they lack sufficient sensitivity to distinguish subtle differences in patients’ functional status. After a stroke, the brain compensates for damaged areas by reorganizing neural networks. Therefore, functional connectivity (FC) between neural networks becomes a powerful indicator to characterize brain damage [6,7,8]. Studies based on functional magnetic resonance imaging (fMRI) and electroencephalography (EEG) have shown that FC at rest correlates with functional impairment of the upper limb and is a valid biomarker for predicting motor recovery [9,10,11]. However, recent studies have shown that FC in the task state (task-FC) more directly reflects differences in functional compensation within the impaired region, whereas these differences are not apparent in the resting state [12, 13]. This may be because the brain allocates more neural resources to task-related processing, reducing noise and nonsmoothness, thereby revealing individual differences more effectively [14, 15].

In this context, effectively utilizing external tasks to elicit functional responses in patients is key to assessing rehabilitation potential. Upper Extremity Robot-Assisted Rehabilitation Training (UE-RAT) simulates real movement tasks by precisely controlling parameters, making the neural responses it evokes more compatible with the demands of daily life. Furthermore, incorporating visual and auditory stimuli in UE-RAT enhances interregional brain interactions and elicits more comprehensive brain functional response patterns [16, 17]. These properties allow UE-RAT-induced neural activity to effectively characterize adaptive changes in the brain following a stroke and reveal subtle differences in rehabilitation potential [18, 19]. However, it is important to note that the application of fMRI and EEG during UE-RAT tasks may be limited due to environmental constraints and sensitivity to motion.

fNIRS is an emerging neuroimaging tool that provides information on regional neural activity by noninvasively monitoring changes in oxygenation and deoxyhemoglobin concentrations [20]. Due to its resistance to motion artifacts and its ability to adapt to freely moving subjects, it has become an ideal method for assessing neural activation during rehabilitation training [21]. Research using fNIRS has demonstrated that patients with varying degrees of upper limb impairment exhibit differences in FC patterns at rest [22]. Further research, such as the studies by Huo [23] and Xie [24], indicates that these differences are more pronounced among patients during UE-RAT tasks. These results suggest that using fNIRS to measure changes in FC in the UE-RAT task using fNIRS can amplify individual differences, which provides an effective tool for accurately assessing rehabilitation potential.

The functional status of stroke patients is critical to rehabilitation prognosis, and the completion of UE-RAT requires activation of the nervous system that controls the patient’s upper limb movements. Therefore, we hypothesize that the neural activity of stroke patients during UE-RAT reflects their functional status. By analyzing these task-related neural activities, it is possible to quantify the patient’s functional status and use it as a quantitative indicator of rehabilitation prognosis. Therefore, the aim of our study is to utilize fNIRS to compute features such as cortical activation and FC in subacute stroke patients both at rest and during UE-RAT training. Subsequently, we will employ the LASSO algorithm to identify robust features and use machine learning algorithms to construct a predictive model for upper limb motor function recovery. Through this approach, we aim to explore the value of the neural activity characteristics induced by UE-RAT in relation to rehabilitation prognosis.

Method

Participant

A retrospective analysis was conducted on eighty-two subacute stroke participants who underwent rehabilitation treatment at the Third Affiliated Hospital of Sun Yat-sen University from October 2021 to May 2023. The inclusion criteria were as follows: (1) onset between 1 week and 6 months, (2) first stroke with a definite diagnosis by cranial CT or MRI, (3) absence of obvious cognitive and language dysfunction (MMSE > 21 points), and (4) age > 18 years. The exclusion criteria were as follows: (1) neurological disorders other than stroke that can cause movement disorders, such as Parkinson’s disease, spinal cord injury, or traumatic brain injury; (2) pregnancy, breastfeeding, or the use of a powered implantable cardiac device to monitor or support cardiac function; and (3) new infarct foci or worsening secondary to bleeding. The trial was registered with the China Clinical Trial Registry (CCTR) under registration number ChiCTR2100054527 on December 19, 2021 and was approved by the Institutional Review Board of the Third Affiliated Hospital of Sun Yat-sen University under number 20,210,233,301.

The following categories of data were obtained from the electronic medical record (EMR) system to form the clinical feature set: (1) demographic data, including age, sex, education level, height, and weight; (2) stroke characteristics, including stroke type, onset time, stroke location, and hemiparesis side; (3) neurophysiological data, including the presence of motor evoked potentials (MEPs), average wave amplitude, and average latency on the healthy/affected side; (4) comorbidities, including history of heart disease, diabetes, hypertension, hyperlipidemia, and hyperuricemia; (5) laboratory data, including the levels of albumin, blood glucose, uric acid, cholesterol, triglycerides, fibrinogen, D-dimer, hemoglobin, white blood cells, and red blood cells; (6) vital signs, including temperature, heart rate, respiratory rate, systolic blood pressure, and diastolic blood pressure; and (7) clinical assessment data, including the Fugl-Meyer Assessment of Upper Extremity (FMA-UE), the Modified Barthel Index (MBI), the National Institutes of Health Stroke Scale (NIHSS), and the Mini-Mental State Examination (MMSE).

Treatment process

Within the first week of admission, clinical measurements, including FMA and manual muscle testing (MMT), among others, were conducted by a licensed physical therapist. Subsequently, the UE-RAT training mode was adjusted based on the MMT results of the participants’ hemiplegic upper limbs. Specifically, muscle strength levels 0–1 corresponded to passive movement (severe group, n = 32), levels 2–3 corresponded to assisted movement (moderate group, n = 35), and level 4 and above corresponded to resistive movement (mild group, n = 15). Throughout the inpatient rehabilitation period, all participants received daily 20-minute sessions of robot-assisted task-oriented upper limb exercise training. In addition to RAT therapy, patients are also required to undergo daily physical therapy, including transfer function training (10 min), balance training (20 min), endurance training (20 min), and neuromuscular electrical stimulation (20 min). Furthermore, they receive occupational therapy consisting of neurodevelopmental therapy (15 min), therapeutic task-oriented activities (15 min), and training in activities of daily living (15 min). These therapies were administered to each patient once per day, five days a week, for two consecutive weeks. The training was administered by professional therapists in a quiet treatment room. Robot-assisted therapy was conducted using the ArmGuider device (ArmGuider, ZD Medtech Co., Ltd., China), a two-degree-of-freedom end-effector robot designed to enhance shoulder and elbow joint flexibility through training in the horizontal plane.

During the preparatory phase of exercise training, participants sat in front of the training platform with their hands and forearms secured to the robotic arm of the device. The participants were then instructed to manipulate the robotic arm along a predefined trajectory to perform motor and cognitive tasks according to the gaming scenario. For example, they controlled a net on the screen using the robot’s moving arm to catch butterflies moving in various directions. During the training, an arrow is displayed on the screen, indicating the direction in which the patient should apply force. After the patient successfully catches the butterfly, a new directional arrow appears, and the process repeats. This game patients will not fail, will eventually catch the butterfly, because it is an assisted training. RAT will automatically detect the patient’s active force, in a single arrival task, when the participant’s active force phase occupies more than 20% of the total trip to successfully catch a butterfly, the screen will show that 10 points will be awarded, and issued a “ding” sound. If the patient’s active participation is less than 20%, and the robot passively completes the butterfly catching, the screen will show that 1 point has been earned, and a “Keep going!” sound will be emitted. Each patient typically completes approximately 200–250 reaching movements per training session, while traditional occupational therapy assists patients in completing around 60–80 upper limb training movements within a 20-minute session.

The Fugl-Meyer Assessment of Upper Extremity (FMA-UE) provides a standardized method for comprehensive evaluation of motor function in the shoulder, elbow, forearm, wrist, and fingers, as well as coordination of the hand and fingers, enabling detailed analysis of upper limb motor function following stroke. Hence, the FMA-UE score at discharge was selected as the primary outcome measure for improvement in upper limb motor function. The difference in scores before and after treatment (∆FMA-UE) was calculated, with improvements exceeding 9 points serving as the threshold for the minimal clinically important difference (MCID) binary classification for subacute stroke patients [25, 26], as the MCID is considered a meaningful clinical improvement that is beneficial for patients’ daily activities postintervention. Participants with ∆FMA-UE ≥ 9 were classified as responders.

fNIRS measurements and preprocessing

During the two-week training period, fNIRS measurements were performed primarily on the first day of UE-RAT training. To ensure the accuracy of the measurements, all participants were asked to sit quietly for 5–10 min prior to the measurements to eliminate hemodynamic responses induced by previous activity. Participants then completed a 10-minute resting state measurement in the seated state followed by a 10-minute fNIRS measurement in the RAT training state. Multichannel tissue oxygenation monitors with continuous-wave technology (NirSmart, Danyang Huichuang Medical Equipment Co., Ltd., China) operating at wavelengths of 740 and 850 nm were utilized for fNIRS measurements. All differential path length factors were initially set to 7.0 with a sampling rate of 10 Hz. The calibration function of the instrument and the corresponding template were used to precisely position the channels to correspond with the 10/10 electrode positions based on different head sizes. A total of 38 measurement channels, comprising 18 light source probes and 16 detector probes, were symmetrically placed over the regions of the ipsilesional and contralesional prefrontal cortex (IPFC/CPFC), dorsolateral prefrontal cortex (IDLPFC/CDLPFC), superior frontal cortex (ISFC/CSFC), premotor cortex (IPMC/CPMC), primary motor cortex (IM1/CM1), primary somatosensory cortex (IS1/CS1), and occipital cortex (IOC/COC), as depicted in Fig. 1.

After acquiring the signals for HbO2 and HHb, data preprocessing was performed according to methods described in our previous studies [24, 27]. fNIRS data preprocessing was performed with the following procedures using customized routines in MATLAB (The MathWorks, Inc.). First, the absorbance signals recorded by fNIRS were subjected to bandpass filtering at 0.0095–2 Hz (zero-phase, sixth-order Butterworth filter) to reduce unrelated noise components and low-frequency baseline drift. Next, principal component analysis (PCA) and independent component analysis (ICA) were conducted on the HbO2 and HHb signals of each channel to identify components potentially related to noise and artifacts, such as cardiac pulsations, respiratory signals, and blood pressure changes. Components exhibiting significant spectra in the range of 0.01–0.08 Hz, corresponding to relevant temporal processes, were visually identified and retained, indicating functional hemodynamic responses in the brain. Finally, a sliding average filter with a time window of 3 s was utilized to remove obvious outliers in the signals, and pseudoartifacts were eliminated through cubic spline interpolation.

Fig. 1

Experimental setup. (A) Schematic representation of near-infrared spectroscopy setup, comprising 18 source probes, 16 detector probes, and 38 measurement channels. (B) Upper limb movement training guided by robot assistance. PFC: prefrontal cortex; DLPFC: dorsolateral prefrontal cortex; SFC: superior frontal cortex; PMC: premotor cortex; M1: primary motor cortex; S1: primary somatosensory cortex; OC: occipital cortex

fNIRS feature extraction

Continuous wavelet transform (CWT) is a wavelet analysis method used for analyzing near-infrared brain oxygenation signals. Here, the Morlet wavelet was utilized to identify oscillatory signals in the range of 0.01–0.08 Hz as hemodynamic responses to neural activity [28]. By scaling the wavelet scale, the spectral information of the frequency can be obtained, and by translating in time, the time information of the desired frequency component can be obtained. At a particular frequency ƒ and time point $\:{\text{T}}_{\text{n}}$, the WT wavelet coefficient is defined as follows:

$$\:{\text{w}}_{\text{k}}\left({\text{T}}_{\text{n}}\right)={\text{W}}_{\text{k}}\left({\text{f},\text{T}}_{\text{n}}\right){\cdot\:\text{e}}^{{\text{i}{\varnothing}}_{\text{k}}\left({\text{f},\text{T}}_{\text{n}}\right)}={\text{a}}_{\text{k}}\left({\text{f},\text{T}}_{\text{n}}\right)+{\text{i}\text{b}}_{\text{k}}\left({\text{f},\text{T}}_{\text{n}}\right)$$

Time-domain averaging of the CWT results produced the wavelet amplitude (WA) of the HbO2 and HHb signals for each channel at each time and frequency. These values reflect the fluctuation amplitude of local cerebral blood flow at a given frequency caused by task-triggered cortical activity and are used as a feature of cortical activation.

Wavelet phase coherence (WPCO) is used to calculate functional connectivity, which describes the statistical interdependence between two hemoglobin oscillatory components by examining how the phase difference between two signals aligns over a specific frequency range [29]. With the CWT, two time series, $\:{\text{x}}_{1}\left({\text{t}}_{\text{n}}\right)$ and $\:{\text{x}}_{2}\left({\text{t}}_{\text{n}}\right)$, in the instantaneous phase under a specific frequency f can be obtained for $\:{{\varnothing}}_{1}\left(\text{f},{\text{t}}_{\text{n}}\right)$ and $\:{{\varnothing}}_{2}\left(\text{f},{\text{t}}_{\text{n}}\right)$, respectively. Thus, the instantaneous phase difference between the two relative oxygenated hemoglobin concentration signals can be expressed as follows:

$$\:{\Delta\:}{\varnothing}\left(\text{f},{\text{t}}_{\text{n}}\right)={{\varnothing}}_{1}\left(\text{f},{\text{t}}_{\text{n}}\right)-{{\varnothing}}_{2}\left(\text{f},{\text{t}}_{\text{n}}\right)$$

Then, based on $\:\text{c}\text{o}\text{s}{\Delta\:}{\varnothing}\left(\text{f},{\text{t}}_{\text{n}}\right)$ and $\:\text{s}\text{i}\text{n}{\Delta\:}{\varnothing}\left(\text{f},{\text{t}}_{\text{n}}\right)$ in the time domain average, $\left\langle \text{c}\text{o}\text{s}{\Delta\:}{\varnothing}\left(\text{f}\right) \right\rangle$ and $\left\langle \text{s}\text{i}\text{n}{\Delta\:}{\varnothing}\left(\text{f}\right) \right\rangle$ can be obtained, resulting in the WPCO defined as follows:

$$\:\text{W}\text{P}\text{C}\text{O}\left(\text{f}\right)=\sqrt{{\left\langle \text{c}\text{o}\text{s}{\Delta\:}\phi\left(\text{f}\right) \right\rangle }^{2}+{\left\langle \text{s}\text{i}\text{n}{\Delta\:}\phi\left(\text{f}\right) \right\rangle }^{2}}$$

The WPCO lies between 0 and 1 and quantifies the degree of agreement between the instantaneous phases of two signals over the continuum of a time series to determine possible connectivity. Higher WPCO values indicate coherence between the two cortical regions, and lower values indicate a weaker relationship between the two cortical signals. More detailed fNIRS features are shown in Table S1.

Selected features

We obtained 38 clinical characteristics from the clinical data. Focusing on fNIRS measurements in the resting state and the first motor training state, we calculated 315 fNIRS features, including activation of brain areas (n = 28), functional brain connectivity (n = 182), and fNIRS difference features derived from subtracting the resting state from the task state (n = 105). The details are presented in Table S1. To obtain robust features, we performed 100 random divisions of the dataset at a 3:1 ratio, with 75% of the data classified as the training set and 25% as the test set. In the training sets, samples were balanced for the classification task using the synthetic minority oversampling technique (SMOTE) [30]. Finally, features were normalized using z scores to mitigate the dominant effect of features at larger numerical scales [31].

To address the challenge of irrelevant or redundant predictors within the model training process, LASSO was applied. LASSO is a regression technique that improves model accuracy by adding a penalty to the regression coefficients. This penalty, controlled by a tuning parameter λ, shrinks some coefficients to zero, effectively performing feature selection and reducing model complexity [32, 33]. We repeatedly conducted feature selection across 100 different training sets, utilizing tenfold cross-validation to select the tuning parameters each time. The final results were statistically analyzed, and features with the highest selection frequency were considered robust. This method provided a way to assess feature stability, ensuring that key features were identified across varying data distributions.

Model construction and evaluation

The feature sets were designated the clinical feature set, fNIRS feature set, and combined feature set, which incorporated both clinical and fNIRS features. Robust features and target variables were selected from these feature sets to construct the model at an empirical ratio of 10:1 for the number of datasets: features, which helps to reduce the risk of model overfitting. The dataset was randomly resegmented 100 times, and in each segmented training set, six machine learning algorithms were used to construct the predictive model: support vector machine (SVM), random forest (RF), logistic regression (LR), k-nearest neighbor (KNN), artificial neural network (ANN) and elastic network (EN). Tenfold cross-validation combined with a grid search was used to select the optimal hyperparameters for the different algorithms. The specific parameter settings are listed in Table S2. The experimental flowchart is shown in Fig. 2. All modeling was carried out using Python version 3.7.15.

For the binary classification model of responders and nonresponders, we calculated the area under the receiver operating characteristic curve (AUC) and accuracy to assess the predictive performance of the different models. For the FMA-UE scores, we calculated the Spearman correlation coefficient (r) and root mean square error (RMSE) between the predicted and true values to evaluate the model’s performance, and the results are presented as the mean ± standard deviation. Finally, we employed Shapley Additive Explanations (SHAP) to elucidate the prediction process of the best-performing model on the test set. The SHAP interpretation elucidates the output of a machine learning model by employing the Shapley value principle from game theory. This method quantifies the influence of each feature on the model’s output, providing insight into how the model utilizes these features to make predictions [34].

Fig. 2

Flowchart of data analysis. Ninety-four participants were included. Cortical activation and connectivity features were measured using fNIRS during resting and UE-RAT tasks at the participants’ first rehabilitation session, and baseline clinical characteristics were collected for model construction, with the predictor target being the FMA-UE score at discharge. After excluding 12 participants with missing FMA-UE scores, 82 participants met this criterion. The dataset of these patients was divided 100 times, and the training set after each division was subjected to feature selection, and the features with the highest screening frequency were used to build clinical, fNIRS, and combinatorial models. FMA-UE: Fugl-Meyer Assessment for Upper Extremity; MEP: Motor Evoked Potential; S1: Session 1

Statistical analysis

All the statistical analyses were performed using SPSS version 22.0 (IBM, Armonk, NY). The normality of the data was assessed using the Kolmogorov‒Smirnov test. Depending on the distribution (normal or nonnormal), continuous data are presented as either the mean ± SD or the median and interquartile range (IQR). Differences in continuous variables were compared using either the Mann‒Whitney U test or the independent samples t test. Fisher’s exact test was used to compare differences in categorical variables. Categorical variables are described by the number of participants (percentage), and two-tailed tests with a P value less than 0.05 were considered to indicate statistical significance.

Results

Of the 94 participants who underwent upper limb robotic rehabilitation treatment, 12 participants were excluded due to accidental discharge, lack of assessment, or incomplete medical records. Ultimately, a total of 82 subacute stroke participants were included. Among the participants with FMA-UE scores less than 21, 29 were in the passive group, 17 were in the assisted group, and 4 were in the resistance group. All participants showed improvements in FMA-UE scores after training. In the responder group (n = 20), the average FMA-UE score was 18.60 ± 17.13 before training and 36.90 ± 19.23 after training. In the nonresponder group (n = 62), the average FMA-UE score was 22.19 ± 17.79 before training and 24.79 ± 17.64 after training. A summary of the demographic characteristics, stroke features, and clinical test results of the participants is presented in Table 1. Significant differences (p < 0.05) were found in the time since onset and baseline MBI between responders and nonresponders. In the responder group, the average time of onset was 48.50 ± 63.00 days, and the MBI score was 45.45 ± 22.58 points.

Table 1 Baseline characteristics of participants

Modeling feature selection

Feature selection was performed for classification models predicting participants with clinically significant improvement and regression models predicting continuous FMA-UE values, as shown in Fig. 3. For the classification model, task-state FC features (such as IPFC-CDLPFC, CM1-IOC, CPFC-OC, and CSFC-IOC) were most frequently selected, with frequencies exceeding 90%. In contrast, the selection frequency of resting-state and difference features was less than 88%, indicating that task-state features are more robust. For more detailed results see Fig. S1. The time since stroke onset was the most significant clinical feature, with a selection frequency of 100%; the baseline MBI score and diabetes status were also found to be highly reliable predictive factors. In the regression model, key fNIRS predictors included resting-state FC for CSFC-IOC and task-state FC between the CDLPFC-COC and IPFC-IPMC.

These robust features were used to construct the models. For classification, both the fNIRS model and the clinical model utilized the 9 features with the highest selection frequencies. The combined model included IPFC-CDLPFC-T, CM1-IOC-T, CPFC-OC-T, CSFC-OC-T, COC-T, CSFC-PMC-T, time since stroke onset, MBI score, and diabetes status. In the regression model, the combined model used CDLPFC-COC-T, IPFC-IPMC-T, CSFC-IOC-R, baseline FMA-UE score, time poststroke, and baseline NIHSS score.

Fig. 3

Feature selection results. (A) An example of feature selection using the LASSO model, illustrating the tuning of the parameter λ. The red dot represents the average mean squared error (MSE) computed from a 10-fold cross-validation. Gray vertical lines indicate the standard deviation of the MSE, while the dashed vertical line marks the optimal value of λ corresponding to the minimum average MSE, a critical point for feature selection. (B) Coefficient paths for feature selection via LASSO, with each colored line representing the coefficient trajectory of an individual feature as a function of λ. Features with nonzero coefficients, determined by the optimal λ marked by the vertical dashed line, show a trend toward zero beyond this point, indicating stabilization in feature selection. (C) Frequency of feature selection targeting MCID classification by FMA-UE. (D) Frequency of feature selection targeting continuous values of FMA-UE

Model outcomes

In the classification models, we initially constructed task-state, resting-state, and difference models based on the fNIRS feature selection results from these three states, as shown in Table 2. We found that the model built solely using task-state data consistently outperformed those based on resting-state or difference features in prediction accuracy, highlighting the critical role of task-state measurements in evaluating stroke recovery potential. However, we observed that the robustness of certain features was not ideal. For instance, key features such as CS1-IOC-T and IDLPFC-IS1-T were selected less than 80% of the time, suggesting greater variability in their performance across data partitions. Therefore, we selected the highest frequency features from task-state, resting-state, and difference features to enhance the robustness of the model.

The fNIRS, clinical, and combined models for classification and regression predictions were constructed using robust features selected by LASSO. The results of the constructed models are shown in Table 3; Fig. 4. Among these algorithms, the fNIRS models based on ANN and SVM exhibited superior performance, achieving moderate predictive capabilities: the ANN model had an AUC of 0.773 (0.111), and the SVM model had a correlation coefficient (r) of 0.461 (0.153), confirming the effectiveness of fNIRS features in rehabilitation prediction (Table 4; Fig. 5). Furthermore, integrating fNIRS features with clinical characteristics significantly enhanced the performance of both the ANN and SVM models: the AUC for the ANN model increased to 0.861 (0.087), and the correlation coefficient for the SVM model increased to 0.862 (0.069). These results underscore the importance of combining fNIRS and clinical features, which significantly boosts the accuracy of rehabilitation predictions. For the combined model, we also compare the modeling results of the original data and the balanced processing of the SMOTE data, as detailed in Table S3. The comparison reveals that the SMOTE sample balancing treatment improves the model prediction results in most cases, especially in the sensitivity, which achieves a significant improvement, and the key metrics, such as accuracy and AUC, are also improved. This indicates that after applying sample balancing to the original data, the model’s ability to identify the minority class has significantly improved, which contributes to enhancing the model’s predictive performance.

Table 2 Classification outcomes of MCID patients were predicted using machine learning models constructed from resting-state, task-state, and the difference between the two fNIRS datasets

Table 3 Classification results of predicting postintervention MCID participants using six machine learning methods with three types of datasets

Fig. 4

AUC results of classification analysis (∆FMA-UE < 9 vs. ≥ 9). Violin plots display the median, first and third quartiles, and minimum and maximum values of the AUC distributions calculated by six machine learning methods. (A) fNIRS feature dataset, (B) clinical dataset, (C) fNIRS + clinical dataset

Table 4 Correlation coefficient (r) and root mean square error (RMSE) of model predictions for postintervention FMA-UE score

Fig. 5

Regression analysis results for root mean square error (RMSE) predicting postintervention FMA-UE score. Violin plots display the median, first and third quartiles, and minimum and maximum values of the RMSE distributions calculated by six machine learning methods. (A) fNIRS feature dataset, (B) clinical dataset, (C) fNIRS + clinical dataset

Predictive factor analysis

We recorded the Shapley values of the features in each prediction model and then calculated the mean to obtain the average feature importance, as illustrated in Fig. 6A and B. IPFC-CDLPFC-T, CM1-IOC-T, CPFC-COC-T, and CSFC-IOC-T were the features with the greatest decision contributions in the classification model. Figure 6C shows that participants with higher baseline values of IPFC-CDLPFC-T and CM1-IOC-T are more likely to achieve significant functional improvement from training, whereas lower values of CPFC-COC-T and CSFC-IOC-T suggest limited recovery. The FMA-UE and NIHSS scores at baseline, as well as CDLPF-COC-T, contributed the most to the regression model. Plus and minus signs indicate the direction of a feature’s influence, with a positive value indicating an increase in model output and a negative value indicating a decrease.

SHAP analysis also allows visualization of the prediction process for individual participants to understand how these features play a role in prediction. Examples of individual predictions from the classification and regression models are shown in Fig. 6E and F. The horizontal axis represents the Shapley values, with 0.517 and 26.91 denoting the “base values” of the contributions of the classification and regression models, respectively. These values represent the mean or expected output (predicted value) of the models in the absence of any feature influence. Figure 6E shows that the overall contribution of the Shapley values for this patient’s features (0.572) is greater than the expected value (0.517), classifying it as positive, with CSFC-CPMC-T and CSFC-IOC-T being the main contributing features. Figure 6F shows that a higher baseline FMA-UE score positively affects the predictions, whereas a longer time since stroke diminishes the model’s predictive output.

Fig. 6

Interpretability analysis of model predictions. (A) Average contribution of input features to the predictions of the best classification model (ANN). (B) Average contribution of input features to the predictions of the best regression model (SVM). (C) Swarm plot showing the average Shapley values and feature distributions for the classification model identifying responders (∆FMA-UE ≥ 9), with the color of the points indicating the positive or negative contribution of the feature values to the prediction. (D) Average Shapley values and feature distributions for the regression model predicting patient FMA-UE scores. Example of model output for a single patient, showing how the Shapley value calculated from the patient’s eigenvalues determines the categorical prediction result (E) for responders and the regression prediction result (F) for the FMA-UE value (true value of 26). E[f(X)] denotes the base value of the model, f(X) denotes the output value of the model

Discussion

This study used fNIRS to assess cortical activation and FC changes in subacute stroke patients in the resting state and during UE-RAT. These fNIRS features, combined with demographic and clinical data, were used to develop models through machine learning methods to predict upper limb motor function recovery after two weeks of intensive training. The main findings of the study are as follows: (1) Compared to resting-state measurements, the FC induced by UE-RAT demonstrates superior performance in predicting short-term upper limb motor function recovery. (2) combining task-related FC with clinical features further enhances the accuracy of the predictive models. Additionally, through model interpretation analysis, this study revealed that initial functional brain responses during UE-RAT were associated with short-term recovery of upper limb motor function. For example, patients with greater connectivity between regions such as the IPFC-CDLPFC, CMC-IOC, and CSFC-CPMC have greater short-term rehabilitation potential. This information is crucial for enhancing the accuracy of rehabilitation predictions.

Previous research relying on traditional statistical models achieved AUCs of 0.72 and 0.58 for predicting patients with clinically significant improvement (ΔFMA-UE ≥ 9) based on hand motor scales and time since onset, respectively [25]. Additionally, time since onset was not a good predictor of the FMA-UE score at discharge, with a percentage of variance explained (R2) < 35% [35]. In another study [36], kinematic data measured by robots were used to predict FMA-UE scores, with correlation coefficients ranging from 0.65 to 0.82. We used combined models constructed with ANNs and SVMs that obtained an AUC of 0.861 (0.087) and a correlation of 0.862 (0.069), respectively. Compared to previous studies, we achieved better accuracy, which may be attributed to the modeling methods and feature parameters used. Rehabilitation is influenced by multiple factors, and these influences are not simply linear changes but rather a complex nonlinear process. Machine learning effectively handles and analyzes complex multidimensional data, particularly in situations with large datasets and high feature dimensions, as it can automatically capture hidden patterns and nonlinear relationships within the data. This advantage allows machine learning methods to perform better outcomes for rehabilitation prediction [37, 38], and our study supports this finding. However, it is essential to emphasize that although machine learning outperforms traditional methods in many tasks, it does not imply it is the best choice in all cases. For instance, in the classification of responders versus non-responders shown in Table 4, the predictive performance of LR outperforms the RF, KNN, and EN models, indicating that traditional statistical models may have advantages in dealing with linear pattern data and smaller sample sizes [39]. Thus, when selecting models, it is crucial to consider the specific characteristics of the task and the distribution of the data to ensure the most appropriate model and method are chosen.

Compared to bedside assessments and demographic clinical information, evaluating neuron activity or functional abnormalities induced by stroke may more directly reflect the patient’s neurological status. Using discriminant function analysis, Hu [40] found that compared to clinical variables, resting-state FC based on fMRI had greater discriminative power in predicting the prognosis of patients with Moyamoya disease (MMD). Chen [41] identified biomarkers for upper limb rehabilitation from clinical variables and EEG data and determined that brain network connectivity features were the best predictors. Bian [42] and Mohanty [43] used machine learning to construct classification models based on standardized tests, and their results emphasized the importance of resting-state brain functional connectivity in predicting rehabilitation outcomes, providing more accurate predictions. However, motor recovery is a complex process influenced not only by spontaneous recovery but also by intervention-induced neuroplastic changes. Resting-state measurements cannot directly provide neuroimaging information about how the brain responds to treatment. Several fMRI studies have suggested that task-related brain functional network reorganization provides valuable insights and that identifying intervention-specific biomarkers may more effectively predict upper limb motor recovery [12, 44, 45]. Cole [46] used FC-based machine learning models to predict task activation, and the results showed that task-related FC changes significantly improved the accuracy of predicting cognitive task activation. Another study used FC data from both resting-state and multiple task conditions to predict fluid intelligence scores, and the results indicated that task-based models enhanced the prediction of individual characteristics [47]. Sutoko [48] used brain activation and connectivity features measured by fNIRS to classify children with attention deficit/hyperactivity disorder (ADHD). The results showed that task-based classification outperformed baseline measurements (AUC: 0.799–0.862 vs. 0.845–0.853). Additionally, connectivity features performed better than activation-based biomarkers. This result is consistent with our findings, where FC measured during tasks is highly important and plays a crucial role in enhancing the accuracy of rehabilitation predictions. Overall, FC provides valuable predictors at the level of brain function by assessing neuronal activity and functional brain networks. Task state measurements capture the complex changes in the brain during task performance, and these patterns of change are closer to the actual functional needs of the patient, improving the accuracy of rehabilitation predictions.

In recent years, treatments involving UE-RAT that incorporate visual or auditory stimuli or VR technology have shown better clinical outcomes in improving upper limb function [49, 50]. This finding may be because motor tasks accompanied by multimodal stimuli during training induce more robust cortical responses and strengthen connections within brain interaction networks [16, 51]. Jiang’s study [52] confirmed that the fNIRS responses induced by UE-RAT demonstrate good temporal and spatial reliability (R2 > 0.6), making them effective biomarkers for characterizing UE-RAT-induced neural responses. We observed a potential association between the extensive brain network connectivity induced by this training modality and, in particular, the initial level of functional connectivity of these brain networks and the recovery of upper limb function in stroke patients: higher connectivity of the IPFC-CDLPFC, CMC-IOC, and CSFC-CPMC during task states may indicate greater potential for functional improvement. The occipital cortex plays a critical role in integrating multisensory information, including spatial attention, decision-making, sensorimotor integration, and motor planning during visually guided movements [53]. On the other hand, the primary motor cortex is primarily responsible for executing motor commands and coordinating muscle activities [54]. The enhanced interaction between the CMC and IOC may reflect the brain’s effort to compensate for sensory loss by strengthening connections between the ipsilateral sensory cortex and contralateral motor control regions. Moreover, greater connectivity suggests that the brain reinforces processing of visual information and motor execution, showing that robot-assisted movement rehabilitation also involves other modalities, such as cognitive or visual processing [55].

Previous research has shown that robot-assisted motor training can improve brain regions associated with motor learning and attention, such as the superior frontal cortex (SFC) [24]. The SFC is important for cognitive functions, including attentional control, decision-making, and motor planning [56]. The connectivity of the CSFC with the CPMC reflects the role of the contralateral brain region in compensating for impaired ipsilateral function after stroke [57]. We found that patients with greater connectivity between the CSFC and CPMC were more likely to benefit from training, possibly because UE-RAT, combined with visual and auditory stimuli, provides rich sensory inputs that enable more effective integration of cognitive and motor planning to guide and optimize motor execution [58]. A study validated the efficacy of this training paradigm in improving upper limb function and cognition in acute and subacute stroke patients [59]. The increased connectivity between brain networks suggested a modification in the pattern of brain region recruitment during multisensory stimulation. This modification enhances the interaction among sensory-motor, visuospatial, and cognitively associated motor and nonmotor areas, helping patients acquire more clinically relevant motor skills [60, 61]. However, patients performing more complex motor tasks may require greater attention and sensorimotor processing to integrate visual, proprioceptive, and tactile feedback information with motor output [62]. The study by Ward [63] revealed that severely impaired poststroke patients exhibited an enhanced activation trend in the contralateral occipital lobe during task execution, which was significantly negatively correlated with early motor performance. Premature performance on complex tasks may lead to overloading of the contralateral hemisphere, particularly an overactive state in the contralateral occipital lobe, which could have an unfavorable effect on the patient’s future functional recovery. Several studies have also reported increased activation in the occipital lobe of stroke patients, which is linked to various symptoms and consequences, including visual field defects, cognitive function recovery, and visual impairments [64, 65].

The results of this study also highlight key clinical features, such as time since stroke, as critical in predicting patient functional improvement. Sale [66] and Lee [25] et al., among others, have suggested that upper limb robotic-assisted therapy may result in greater functional improvement in stroke patients with a shorter time to onset of stroke and that the potential for recovery may diminish as time since stroke progresses [67]. At the same time, baseline clinical tests (e.g., FMA-UE, NIHSS, MBI) were important in at least one of the prediction models, and for the regression task in particular, the baseline score was identified as the most important predictor, consistent with the findings of previous studies [3, 68]. This result suggests that the initial level of impairment may need to be taken into account in the prediction of stroke rehabilitation. We also found that diabetes was an important factor contributing to poor recovery in subacute stroke patients. A recent meta-analysis [69] has shown that diabetes is commonly associated with poorer functional outcomes after stroke. This association may be because diabetes hinders the formation of new blood vessels in the brain, thus failing to meet the brain’s need for proper blood flow and nutrient delivery [70]. Our study suggests using both regression and classification models together, as the predicted values generated by the regression model can serve as a valuable supplement to the classification results. For instance, even if certain patients are classified as nonresponders by the classification model, the regression model can still indicate minor progress in their rehabilitation therapy, aiding in the formulation of subsequent treatment plans or adjustments to rehabilitation strategies. The combined use of these two models can offer a more comprehensive and nuanced assessment of therapeutic efficacy.

Study limitations

This study provides valuable insights into the use of fNIRS measures and machine learning methods to predict the response to upper limb robotic rehabilitation therapy in stroke patients. However, several limitations remain. First, the relatively limited number of subjects in the study may have restricted the statistical power of some analyses. To ensure the reliability and applicability of the model, expanding the sample size and conducting external validation are necessary steps for further research. Second, the study adopted a retrospective design, and although patients received similar treatment plans, additional conventional physiotherapy may have been administered during the rehabilitation process, potentially impacting the results. Another limitation is the lack of motor data related to the quality of upper limb robotic training completed by patients, which is important for interpreting patient participation and rehabilitation outcomes. Additionally, we only evaluated the neural responses during the patients’ first training session. Dynamic monitoring at multiple time points could potentially reveal the recovery process and trends, providing a more accurate prognosis. Finally, we would like to point out that this study only included subacute stroke patients for short-term recovery predictions. The spontaneous recovery during the subacute phase may obscure more complex rehabilitation patterns, so our results may not be applicable to chronic stroke patients. Future research should consider predicting long-term recovery in subacute stroke patients and conduct related studies in chronic stroke patients to validate the effectiveness of UE-RAT task-related brain functional responses in rehabilitation predictions.

Conclusion

In this study, we used fNIRS to explore the potential associations between brain functional responses induced by UE-RAT and upper limb motor function recovery. By integrating machine learning algorithms, we constructed a model that significantly improved the prediction accuracy. Our findings emphasize that brain functional response features measured during task states, particularly FC, are more valuable than resting-state measurements in predicting short-term rehabilitation outcomes. These results suggest that measuring baseline brain responses induced by UE-RAT with fNIRS and employing machine learning for recovery prediction could be effective methods for assessing patients’ rehabilitation potential. Although these initial results are encouraging, further comprehensive and extensive research is needed to validate and expand these findings.