A total of 79 first-time stroke patients with hemiplegia who were hospitalized in the Rehabilitation Medicine Department of The Affiliated Brain Hospital of Nanjing Medical University from April 2018 to December 2021 were prospectively included in this study. The study protocol was approved by the Ethics Committee of The Affiliated Brain Hospital of Nanjing Medical University, and written informed consent was obtained from all patients or their family members.

The inclusion criteria of the study included: (1) diagnosis of first-ever subcortical ischemic or hemorrhagic stroke by computed tomography or magnetic resonance; (2) patients were within 2–12 months after stroke; (3) patients only had unilateral limb hemiplegia (Modified Brunnstrom classification as grade I–IV); (4) patients did not have severe cognitive impairment; (5) patients did not have other acute diseases or serious complications; (6) right-handedness. The exclusion criteria were: (1) disturbance of consciousness, severe hearing and visual impairment; (2) significant pain in the affected side (Ten-point Visual Analog Scale > 4); (3) severe primary heart, lung, liver, kidney or hematopoietic system diseases; (4) after craniectomy or cranioplasty. (5) MRI contraindications.

The hospital records of the 79 patients were retrospectively reviewed for data collection. The following data on demographics and relevant past medical history were collected: age, sex, (hospitalization time, date of stroke onset and lesion location). Information on motor functional outcomes were also collected, using the following the assessment tools: Brunnstrom recovery stage (BRS) (Brunnstrom, 1966), modified Ashworth scale (MAS) (Bohannon and Smith, 1987), Barthel Index (Mahoney and Barthel, 1965), Fugl-Meyer Assessment (FMA) (Fugl-Meyer et al., 1975), functional ambulation category (FAC) (Holden et al., 1984), and Semans balance scale (Wang, 2018).

The functional classification of stroke patients was based on the BRS scale, which includes three parts: upper limb, lower limb and hand, with each part having six levels (0–6). The higher the level, the stronger the motor ability (Brunnstrom, 1966). The scale has good internal consistency in stroke patients. The MAS scale was used to evaluate spasticity in the biceps brachii and triceps surae muscles, which was graded out of six levels (0, 1, 1+, 2, 3, 4) (Bohannon and Smith, 1987). The higher the level, the more severe the spasm. The Barthel Index was used for the ability of daily living, with a total score of 100. The higher the score, the better the quality of life (Mahoney and Barthel, 1965). FMA was used to evaluate motor function, with a total score of 100, including 66 for upper limbs and 34 for lower limbs. The higher the score, the better the motor ability (Fugl-Meyer et al., 1975). FAC was used for assessing walking ability, and was divided into six levels—the higher the level, the better the walking function (Holden et al., 1984). The Semans balance scale was used to assess balance. It is an observational assessment method mainly used in the balance assessment of pediatric cerebral palsy and hemiplegic patients after stroke (Wang, 2018). It classifies balance into eight levels, with higher levels corresponding to better balance.

Image scanning was performed on the same day as functional assessment, on a 3.0-T MRI scanner (Siemens, Verio, Germany) in the Affiliated Brain Hospital of Nanjing Medical University. All patients lay supine with their head fixed by foam pads with a standard birdcage head coil to minimize head movement. Participants were instructed to remain as still as possible, open their eyes, remain awake, and not think of anything. High-resolution T1-weighted images were acquired by 3D magnetization-prepared rapid gradient-echo (MPRAGE) sequence [repetition time (TR) = 2,300 ms; echo time (TE) = 2.85 ms; flip angle (FA) = 9 degrees; matrix = 256 × 256; field of view (FOV) = 256 × 256 mm²; slice thickness/gap = 1/0.5 mm; 176 slices covered the whole brain] for image registration and functional localization. The imaging took ~260 s. The resting-state of functional MRI (rsfMRI) were subsequently collected in the same slice orientation with a gradient-recalled echo-planar imaging pulse sequence (TR = 2,000 ms; TE = 30 ms; FA = 90 degrees; matrix = 64 × 64, FOV = 240 × 240 mm²; thickness/gap = 4.0/0 mm; voxel size = 3.8 × 3.8 × 4 mm³; slice numbers = 30). A total of 251 volumes were obtained in this acquisition sequence and each functional resting-state session lasted ~500 s.

Diffusion-weighted imaging was performed with the SE-EPI sequence scanning. Using iPAT technology, axial scanning, TR = 11,000 ms, TE = 90 ms, slice number = 67, slice thickness = 2 mm, slice spacing 0 mm, FOV = 256 mm × 256 mm, matrix: 128 × 128, voxel size = 2 × 2 × 2 mm³, diffusion weighted scan b = 1,000 s/mm² in 30 gradient directions, and another non-diffusion weighted image with b = 0. Excitation times (NEX) = 2. Scan time = 11 min and 57 s.

Diffusion weighted imaging (DWI) was processed using the Infinitome software (Omniscient Neurotechnology, 2020). Infinitome relies on standard processing steps using Python (Garyfallidis et al., 2014). The diffusion image was resliced to obtain isotropic voxels. A rigid body alignment was used for motion correction, and slices with excess motion, defined as DVARS >2 sigma from the mean slice were eliminated. The T1 image was then skull stripped using a convolution neural net, described elsewhere (Isensee et al., 2019). This was then inverted and aligned to the DWI image using a rigid alignment, which was then used as a mask to skull strip the DWI image. A diffeomorphic warping method was then used to correct for gradient distortion to similarize the DWI and T1 images, and eddy current correction was performed. Finally, the fiber response function was estimated and constrained spherical deconvolution was used calculate the diffusion tensors. Deterministic tractography was then performed, with random seeding, manifesting in around 300,000 streamlines per brain (Doyen et al., 2022).

A parcellation scheme is necessary for measuring FC changes at an anatomical level. Most available atlasing methods to parcellate the brain into functional regions however are derived from healthy subjects. These atlases also are based on a group average of these healthy cohorts, and are therefore unable to account for individual morphological differences such as gyral variation. This is especially problematic when applying these atlases to brains altered by pathology, such as in stroke. In order to account for any morphological deficits caused by stroke, we utilized a machine learning based method to individually parcellate each brain based on the Human Connectome Project Multimodal Parcellation (HCP-MMP1) atlas (Glasser et al., 2016), creating a subject-specific atlas. While we have described this method in detail elsewhere (Doyen et al., 2022), briefly, a machine learning model was trained using DWI data from 178 healthy controls from the Schiz Connect database, which were pre-processed as above. The model consequently learnt the structural connectivity pattern between each voxel. The unalted HCP-MMP1 was then warped to each brain in the current study's cohort, and the machine learning model was applied to each brain to re-appoint voxels to their most likely parcellation based on structural connectivity features. This method manifested in “reparcellation” of voxels, creating a subject-specific version of the HCP-MMP1 atlas with 180 cortical regions, 9 subcortical regions, and the brainstem.

Each brain region was then automatically mapped to a known large-scale brain network by the Infinitome software, which itself relies upon previous meta-analyses mapping each large-scale brain network. The network template described by Yeo et al. (2011) was utilized, including the core networks: the Central Executive Network (CEN), Default Mode Network (DMN), Dorsal Attention Network (DAN), Limbic Network (LN), Salience Network (SN), Sensorimotor Network (SMN), and the Visual Network (VN), along with several networks which are either part of the extended versions of the core networks, or additional networks described in the literature, including the Accessory Language and Language Networks (part of the extended DMN), Auditory System (part of the SMN), Multiple demand network, and Ventral Attention Network (VAN).

rsfMRI images were pre-processed prior to analysis according to standard pre-processing steps. A rigid body alignment was used to perform motion correction on the T1 and blood-oxygen-level-dependent (BOLD) images. Slices with excess movement (defined as DVARS > 2 sigma from the mean slice) were eliminated. A CNN was used to skull strip the T1 image (Isensee et al., 2019). In order to then skull strip the rsfMRI image, the T1 image was inverted and aligned to the BOLD image using a rigid alignment, and then used as a mask. Slice timing correction, global intensity normalization, and gradient distortion correction were performed, the latter using a diffeomorphic warping method to register the rsfMRI and T1 images. The CompCor methodwas used to calculate high variance confounds (Behzadi et al., 2007), which were regressed out of the rsfMRI image along with the motion confounds. The linear and quadratic signals were detrended (note this method does not perform global signal regression). A 4 mm FWHM Gaussian kernel was used to perform spatial smoothing. The subject-specific brain atlas created in the previous steps were then registered to the T1 image, and the regions were aligned with the regions in each subject's scan. A visual check was then performed by two neuroanatomists, independently, to ensure the methodology accounted for the morphological changes caused by the lesion. The atlas was therefore ideally positioned to extract a BOLD time series, averaged over all the voxels within a region, from 379 brain regions (180 cortical regions from two hemispheres, along with 19 subcortical structures). A Pearson correlation coefficient was then calculated between the BOLD signals of each unique pair of regions, yielding. Thus, it is ideally positioned for extracting a BOLD time series, averaged over all voxels within a region, from all 379 regions (180 parcels from two hemispheres, plus 19 subcortical structures). The Pearson correlation coefficient is calculated between the BOLD signals of each unique area pair (self to self-inclusive), which yields 71,631 unique correlations.

Machine learning was used to model the test performance of each participant based on the pairwise functional correlation or structural connectivity between the 379 regions of each individual's brain atlas. XGBoost Classifier (Chen and Guestrin, 2016), a boosted trees approach was used to fit each model. All models included age and sex as nuisance predictors. A total of 16 models were trained—one utilizing structural connectivity (SC), and another functional connectivity (FC)—to predict performance on the BRS on upper extremity, FMA on upper and lower extremities, MAS on upper and lower limbs, Barthel Index, FAC, and Semans balance scale. The scores for each of these tests were binarized for the models to classify subjects into either category, as detailed in Table 1.

Five-fold cross-validation was used for each model, with different training and test data splits to ensure the hyperparameters were not being optimized just for the given training and test data splits. The models were evaluated with the mean area under the receiver operating characteristic curve (AUC-ROC) ± standard deviation. Given the small sample size, several strategies were utilized to minimize the impact of class imbalances on the performance of each model. Five-fold cross-validation ensured that each fold had an equal ratio of both classes; a stopping criterion was applied to stop training if the model performance did not improve over consecutive iterations of hyperparameter tuning—an attempt to prevent overfitting; and, feature importance analysis only relied upon the correctly predicted cases, therefore being independent of whether one class was predicted slightly better than the other. These strategies aimed to approximate a balanced class AUC even in cases of class imbalance. Our predominant aim in each model was also feature importance analysis, with a focus on identifying what resulted in the classifications. Consequently, while AUC is inherently affected by class imbalances, the guardrails applied, along with the significantly higher than chance performance of each model counteract any AUC bias, and make it largely irrelevant to the aim of our analysis.

The feature importance analysis for each model was visualized using a SHAP plot of the top 20 features, either functional or structural connectivity between two regions of the HCP-MMP1 atlas, contributing to the model. The SHAP method, based on cooperative game theory, calculates Shapley values for each feature, which estimate the marginal contribution of each feature to the outcome of models relying on every possible permutation of features (Lundberg and Lee, 2017). Each SHAP plot shows a list of features in descending order of importance, quantifying the impact on the model along the x-axis, with the color of each point representing a single observation, indicating whether a high (red), or low (blue) value of that feature is associated with the model. A network-based summary plot is also produced, aggregating the contribution of each parcel by their network affiliation. This gives an indication of which brain networks are most influential in each model's classification.

Subject demographics are detailed in Table 2. The median age (±IQR) of our sample was 62 ± 18. There were 22 females (27.8%) and 57 males (72.2%). The median time (±IQR) from initial stroke to scanning was 35 days (33). There was a small difference in the stroke side, with 55.7% of patients having a left-sided stroke; however, speech deficits were only seen in 21.5% of patients, and most patients had a subcortical stroke. Motor deficits were present in 96.2% of patients.

Machine learning was applied to classify subjects into binarized categories as outlined in Table 1. For each test, two models were constructed to separately evaluate functional and structural connectivity. Both models included age and sex as covariates.

Although our models pointed to several large-scale networks which are known to be important in functional recovery after stroke, such as the DMN, CEN, and SMN (Tuladhar et al., 2013; Larivière et al., 2018; Wu et al., 2020; Olafson et al., 2021; Vicentini et al., 2021), the attention networks, the DAN and VAN were among the top contributors in several of our models. Classically, The DAN is considered to be crucial for sustaining attention on a given task (Ptak and Schnider, 2010), while the VAN is thought to act as a circuit breaker by redirecting attention in response to salient events (Corbetta and Shulman, 2002; Corbetta et al., 2008). Generally, impairment in the DAN following stroke has been associated with neglect (Corbetta and Shulman, 2011; Barrett et al., 2019) and the severity of motor impairment (Siegel et al., 2016), however, the exact nature of changes within the attention networks in stroke patients remains unknown. Some studies have reported increased connectivity (Rehme et al., 2012; Wang et al., 2014; Siegel et al., 2016; Lee et al., 2018), contextualized by the greater need for patients to direct attention to movements following stroke. While others have demonstrated lower FC within the VAN and DAN in both task-free and task-based settings (He et al., 2007; Baldassarre et al., 2016; Adhikari et al., 2017; Barrett et al., 2019). This may be pointing to increased inter-network connectivity and decreased intra-network connectivity, representing the loss of network segregation and specialization following stroke (Caliandro et al., 2017; Marebwa et al., 2017; Guo et al., 2019). Interestingly, congruent to our findings, Romeo et al. recently demonstrated a positive correlation between decreased connectivity of the DAN and the motor component of the Functional Independence Measure (FIM) in the sub-acute to chronic period following stroke (Romeo et al., 2021). Furthermore, several studies have shown that the outcome of motor rehabilitation may be associated with preserved FC within the ipsilesional DAN (Cheng et al., 2021; D'Imperio et al., 2021). Therefore, patients with preserved connectivity within the DAN and VAN may be performing better across the tests in our models. We however cannot validate this given the task-free nature of our study. Nonetheless, it is evident that the attention networks play a role in mediating functional performance, though further studies on large, homogeneous samples are necessary to elucidate the mechanism behind this, which may provide further therapeutic avenues.

Interestingly, the language and accessory language networks were highlighted in our models classifying the BRS, FMA, Barthel Index, Semans balance scale and FAC. This may be due to impairments in the language system typically accompanying strokes with worse motor and functional impairment, and indeed those with expressive aphasia tend to have an 80% chance of also suffering hemiplegia (Romeo et al., 2021). Due to our small sample size, we were unable to perform separate analyses on patients reporting a brain lesion in the left compared to the right hemisphere. However, given most of the strokes within our cohort were subcortical, only nine out of 17 strokes affecting speech were left-sided cortical strokes, and most patients with motor deficits did not have speech deficits, we speculate that there may be alternative reasons for why the language network contributed to our classifications. This finding may instead be pointing to the functional symbiosis between language and motor function (recently reviewed by Anderlini et al., 2019). Recently, several studies have demonstrated that not only do deficits in these domains co-occur, but they also interact in their treatment. For example, Romeo et al. also demonstrated an association between the FIM motor score and the language network, specifically a higher integration within the region of the precentral gyrus, pointing to a possible role of the language network in motor sequence planning (Romeo et al., 2021). Indeed, Maitra et al. have reported an improvement in motor function through self-vocalization (Maitra et al., 2006), while Arya et al. demonstrated inadvertent improvement in language following upper limb therapy (Arya and Pandian, 2014). A similar effect was also found following transcranial direct current stimulation over the motor cortex, which led to improved language outcomes (Meinzer et al., 2016). More recently, Hybbinette et al. also demonstrated that similar mechanism of recovery may be involved in language and motor recovery post-stroke, further contributing to the so-termed “shared recovery hypothesis” (Hybbinette et al., 2021). While well-controlled studies examining this phenomenon are lacking, our findings may suggest either: functional compensation by the language system associated with better performance, or impaired connectivity within the language system manifesting in poor sensorimotor outcomes. Since our models are unable to point to which is the case, further studies are necessary to elucidate the interaction of the language and motor networks following stroke. Nonetheless, both our findings, and the literature highlight the possible need to simultaneously target the language and motor systems for effective neurorehabilitation.

Our models classifying the BRS Upper Limb and FMA Upper Limb both suggested that the limbic/paralimbic system may be playing an important role in motor function following stroke. A similar finding was reported recently by Li et al., who found that the connectivity between the limbic network and DAN was predictive of the FMA (Nishimura et al., 2011). It is therefore possible that the limbic system may be contributing to motor learning after stroke. The limbic system in the Infinitome atlas comprises the orbitofrontal cortex, amygdala, hippocampus, along with other medial temporal regions. Previous studies have established that motor sequence learning relies on the interaction between the prefrontal cortex and hippocampal and striatal networks, along with cortico-cerebellar networks (Albouy et al., 2008; Fernández-Seara et al., 2009; Rose et al., 2011; Burman, 2019; Schapiro et al., 2019; Gann et al., 2021). The activity within the fronto-hippocampal networks decrease as learning progresses (Albouy et al., 2013). As motor learning and re-learning is a key part of neurorehabilitation in the post-stroke period, the limbic system may be key to this process. Alternatively, as the limbic system is also associated with reward and motivation, this finding may be implicating motivation as a key factor in performance. Although one study has examined the role of the motivation in motor recovery following spinal-cord injury in macaques, demonstrating increased functional connectivity between the primary motor cortex and the ventral striatum, orbitofrontal cortex, and anterior cingulate (Nishimura et al., 2011), the role of the limbic network in motor recovery in humans has not garnered much attention. However, it is widely assumed that motivation is a significant component of neurorehabilitation (Widmer et al., 2017; Oyake et al., 2020). Consequently, the limbic system may potentially be a target to improve motor recovery in the post-stroke period, though this remains conjectural and further studies are required to substantiate it.

The prediction performance of our models was comparable to those reported in previous studies examining acute/subacute stroke patients. Wang et al. achieved an AUC of 0.899 using random forest models in acute stroke patients (Wang et al., 2019), while Thakkar et al. had an AUC of 0.77 (Thakkar et al., 2020). It should be noted that direct comparisons of performance are not possible, as most other studies have not combined both DWI based structural and rsfMRI-based FC metrics in machine learning models. Riahi et al. utilized EEG-based FC to estimate the FMA scores in chronic stroke, achieving an R² of 0.97 in their regression model (Riahi et al., 2020). Tozlu et al. compared the performance of several types of models in predicting FMA scores based on a combination of demographic, clinical, neurophysiological, and structural connectivity based metrics, achieving an R² between 0.70 and 0.91 (Tozlu et al., 2020). In order to compare and generalize the performance of our models, future research should incorporate functional and structural connectivity measures into holistic brain connectivity models.

Interestingly, our models utilizing FC tended to have better performance than those utilizing structural connectivity. While this needs to be further investigated and validated in independent cohorts, we speculate that this may be reflective of the rapidity of functional compensation, compared to the structural changes which may occur over a longer time span. This is in line with several studies suggesting that functional reorganization is observed within the first 4–5 weeks following stroke (Golestani et al., 2013; Nijboer et al., 2017; Xia et al., 2021), while structural changes occur from 3 to 12 months (Lin L. Y. et al., 2018); though these timings are contentious. For example, a recent study found no changes in the connectivity of motor areas over a year, despite motor improvement (Branscheidt et al., 2022). This may suggest that functional connectivity changes occur distributively to promote recovery, or that some reported connectivity changes are lesion-related, and not a form of compensation. These discrepancies must be addressed through much larger studies utilizing longitudinal data looking at whole-brain functional connectivity. Though there are a number of outstanding questions, our findings make clear that structural and functional connectome changes are related to functional outcomes, and that machine learning models may be necessary to elucidate the complex patterns of neural change occurring in response to stroke.

Overall, our analysis highlighted several networks which may be playing key roles in motor recovery post-stroke, including the DAN, VAN, the language network, and the limbic system. It is possible that enhancing these networks through non-invasive stimulation techniques, such as repetitive transcranial magnetic stimulation, may improve functional outcomes post-stroke. This is however entirely speculative, and clinical trials are necessary to test these targets. rTMS has already demonstrated benefit for motor recovery the post-acute period, however, treatment often targets the primary motor cortex (Dionísio et al., 2018; Fisicaro et al., 2019). It is therefore necessary to investigate whether a more individualized network-based approach to target selection would improve outcomes. While some of our group has demonstrated that this may be a promising and efficacious technique in a case study (Yeung et al., 2021), larger, controlled trials are needed.

Our study also reveals the need to adopt a global-approach when examining deficits in stroke, and therefore looking beyond classically defined motor-regions when investigating motor deficits. This may consequently allow for the development of deficit-specific brain network—networks of functionally connected regions which may be associated with upper limb or lower limb function, or even a spasticity network. These may then enable for further identification of specific targets for treatment.

Our study is limited by a small, heterogeneous sample, and our findings require larger scale prospective studies for validation. Owing to this, we cannot assume that our results are not biased by overfitting to our sample and class imbalances, though our methodology employs several approaches to mitigate this risk, including 5-fold cross-validation and early stopping criteria. Furthermore, future studies should attempt to replicate functional and structural connectivity-based models in a longitudinal cohort study to investigate whether connectivity changes may reveal pathological and compensatory mechanisms in the acute and chronic stages. Finally, models may further benefit from multimodal inputs, considering baseline information beyond connectomic data as we have done. This may enable the establishment of more individualized diagnostic trajectories for patients.