Prediction of motor outcome after acute stroke using rich-club organizational analysis

Subjects and motor assessment

Stroke patients (n = 16; 11 male; mean age 65.8 ± 11.0; infarct side: 50% left) with first-time acute ischemic stroke were recruited between September 2015 and July 2018 with informed consent. MRI and assessment of motor functions were performed for at less than 1 week (n = 12), and 1 (n = 16), 3 (n = 13) and 6 (n = 9) months after acute stroke. Motor functions were examined by the Upper-Extremity Fugl-Meyer assessment (UE-FM) scale and Barthel index (BI). The UE-FM was developed to quantitatively assess the severity of motor impairment of upper extremity due to hemiplegic stroke, and is based on the well-defined stages of motor recovery³¹. BI was developed to quantitatively assess disability and functional outcome³². All patients received rehabilitation at the Acute Stroke Unit of Queen Mary Hospital after admission due to acute stroke. After an average of 5 days after admission, patients were transferred to the Stroke Rehabilitation Ward of Tung Wah Hospital for more intensive rehabilitation. All patient underwent conventional occupational rehabilitation therapy sessions, including activity of daily living training, upper limb functional training, cognitive perceptual training and functional task training. All procedures were carried out following operational guidelines of Human Research Ethics Committee, and all protocols were approved by the Institutional Review Board of the University of Hong Kong/Hospital Authority Hong Kong West Cluster.

Image acquisition

All MRI scans were performed using a 3.0 T MRI scanner (Achieva TX, Philips Healthcare, Best, The Netherlands) with body coil for excitation and 8-channel head coil for reception. Diffusion tensor imaging (DTI) data was performed using single-shot spin-echo echo planar imaging, consisting of non-diffusion-weighted image (b0) and diffusion-weighted images (DWIs) with b-values = 1000 s/mm² along 32 gradient directions, were acquired with the following parameters: TR/TE = 4000/81 ms, field of view = 230 × 230 mm2, reconstructed resolution = 3 × 3 mm2, 33 contiguous slices with thickness of 3 mm, SENSE factor = 2, number of averaging = 2, total scan time ≈ 5 minutes. TI-weighted images data were obtained by a 3D Magnetization Prepared-Rapid Gradient Echo (MPRAGE) with the following parameters: TR/TE/TI = 7/3.17/800 ms, field of view = 240 × 240 mm², reconstruction resolution = 1 × 1 × 1 mm3, 160 contiguous slices, scan time = 6 min 1 s.

Image pre-processing

All of the pre-processing procedures were performed by SPM12 (https://www.fil.ion.ucl.ac.uk/spm/). Head motion correction was performed by registering DWI to b0 images³³. MPRAGE images were first reoriented taking anterior commissure³⁴ as the origin. The reoriented MPRAGE images were normalized to MNI152 template to obtain a transformation matrix M. MNI152 template was then inverse normalized to MPRAGE images by the inverse of the transformation matrix M⁻¹. The inverse normalized MNI152 template was non-linearly registered to diffusion-weighted imaging and a transformation matrix T was obtained. In this way, an inverse warping transformation from the standard space to the native DTI space can be obtained.

Network construction

Tractography: Diffusion tensor and diffusion metrics were obtained using the Diffusion Toolkit³². White matter tractography was obtained using the TrackVis (http://trackvis.org) with the algorithm of Fibre Assignment by Continuous Tracking (FACT)^36,37 with angle threshold of 45° and random seed of 32. To smooth the corners of the fibres, spline filter was applied. The results of streamlines were used to construct the structural brain network.

Network nodes definition: Automated Anatomical Labelling 2 (AAL2) atlas³⁸ in the standard space was inversely wrapped to the individual DTI native space according to M⁻¹ and T obtained in image pre-processing. 94 Cortical regions (47 for each hemisphere) were obtained and each region was regarded as a node, of these 8 regions were excluded, namely left and right inferior occipital gyrus, left and right superior temporal pole, left and right middle temporal pole and left and right inferior temporal pole, due to variations in brain coverage.

Network link definition: UCLA Multimodal Connectivity package was used to calculate the number and the length of fibres between every two regions. Two regions were defined as connected regions if a fibre occurs between them. The link weight was defined as the fibre count between two regions.

Network topology metrics

MATLAB (2018b) and brain Connectivity Toolbox (https://sites.google.com/site/bctnet/) was used to calculate brain network topology metrics. The metrics included: (1) node degree, defined as the number of connected links of a node; (2) node strength, defined as the sum of the weight of connected links of a node; (3) local clustering coefficient, represented as connection property to neighbourhood of a node; (4) global efficiency, defined as the average inverse shortest path length in the network; (5) local efficiency, defined as the global efficiency computed on node neighbourhoods; (6) node betweenness centrality, defined as the fraction of all shortest paths in the network that contain a given node.

Rich club organization

R language (3.6.0) was used to estimate rich club organization for each subject at each time point using a previously published method¹⁶. Weighted rich club coefficient was calculated by 4 steps. First, for each degree (k), a subnetwork was obtained by extracting nodes with degree more than k and links amongst them. Second, for each subnetwork, the number of links (n) and sum of weights of the links (W) was calculated. Then, the top n strongest weights of the whole network were summed. The weighted rich club coefficient of this subnetwork is subsequently calculated as follows:³⁹ where W_>k represents weights of the links more than k and sum of top n weights of links. However, nodes with lower degree in a network have lower possibilities of sharing links with each other by coincidence, even random networks generate increasing rich club coefficients as a function of increasing degree threshold k. To circumvent with this effect, 1,000 random networks with the same size and same degree distribution of the network of interest were generated. The average of the rich-club coefficients of the 1,000 random networks were calculated. The normalized rich-club coefficient of the network of interest is defined as^19,20

A subnetwork is considered as a rich club organization when ^40,41.

A node is considered as a rich club node when k = {k₁, k₂, …, k_n}, in which the highest k was called highest rich club level. Each node of a rich club organization was given a score according to their highest rich club level. And then after averaging the score of nodes from all participates at the same time point, the top 8 nodes(i.e., 10% of all nodes) were selected as rich club nodes so that number of nodes used to construct the rich club organization avoid different average degree across subjects.

Node and Connection types

There are two types of nodes in a rich club organization, namely rich club nodes and peripheral nodes. On the other hand, there are three types of connections in a rich club organization, namely rich club connections (between rich club nodes), feeder connections (between rich club nodes and non–rich club nodes), and local connections (between non–rich club nodes).⁴²

Network communication

Cost of a link was defined by product of its length and density (the count of streamlines between two brain regions). Communication cost was defined by the product of length and density based on their topological distance⁴². To calculate communication cost of three connections, first, the shortest paths between 84 nodes should be calculated. Then, each link of the shortest paths was divided into three categories (rich club, feeder and local connections). Next, calculate the communication cost of each link by length and density. Last, communication cost of three kind of connections can be obtained by summing the communication cost of links in each kind. The metric ratio of a certain connection kind was defined as the sum of metric of this certain kind divided by the total metric of the whole network. Communication cost/density ratio was defined as the communication cost ratio divided by density ratio, which was the weight of capacity⁴².

Statistical analysis

Statistical analyses were performed by R Language (version 3.6.0) with the function lmer (https://cran.r-project.org/web/packages/lme4/lme4.pdf) and blmer (https://cran.rproject.org/web/packages/blme/blme.pdf). To simplify the statistical analyses on the local structural brain networks, the left and right hemispheres were flipped so that the left hemisphere always corresponded to the ipsilesional hemisphere. Due to attrition, some of the behavioural data, as well as imaging data, were not obtained. Imputation was thus performed on the behavioural data to increase the effective sample size for subsequent statistical analyses. For patient no. 2, 13 and 15, the behavioural data at 6 months after stroke were imputed from those at 3 months. Since patient no. 11 and 16 already made full recovery at the last follow up, full behavioural scores were assumed for the missing time points. After imputation, 10 patients had the behavioural data for all 4 time points. The imputed behavioural data were underlined in Table 1.

To handle data with an unequal number of longitudinal measures, a Bayesian linear mixed model (LMM)^43,44 was used as the main model in the current study so that a posteriori estimation can be maximized in a Bayesian setting. When necessary, linear mixed models were followed by post hoc tests with Bonferroni corrections for multiple comparisons at p < 0.05. This study mainly investigated three aspects: (1) To test whether the motor outcome and network metrics changed with time. In this part, models are established with the motor outcome and network metrics as responses, time as fixed variable, each subject as random variable, and age and gender as covariates. (2) To test the correlation of metrics with motor outcome. In this part, models were established with motor assessment data as responses, network metric to be tested as fixed variable, each subject as random variable, and age and gender as covariates. (3) To find the biomarkers that may predict motor recovery. Models were established with the change of motor outcome (subtract the baseline motor outcome from motor outcome at other time points) as response, each subject as random variable, metric to be tested as fixed variable, and time, age and gender as covariates. After selecting the metrics that may be predictors of response. A linear mixed regression model was established by the method of step wise. In the whole study, a best fitting model was chosen by comparing models using Likelihood ratio test with the principle of the smaller AIC, the better the model⁴⁵. Also, the simplest model was chosen when the gap of AIC values and likelihood ratio test result were small.

Patient demographics

MRI and motor assessment were performed on 12, 16, 13 and 9 patients at within 1 week, 1, 3 and 6 months after first-time acute ischemic stroke, respectively, as shown in Table 1. Linear mixed model was used to test the fixed effect of time on motor function (UE-FM and BI) in the process of stroke recovery. UE-FM (β_time = 3.47; , p < 0.001; Figure 1 A) and BI (β_time = 7.94, , p < 0.001; Figure 1 A) significantly increased with time. Post hoc tests showed UE-FM at 1 (p = 0.01), 3 (p < 0.001) and 6 (p < 0.001) months after stroke were significantly higher than within 1 week after stroke. BI at 1, 3 and 6 (p < 0.001, all) months after stroke were significantly higher than within 1 week, and those at 3 and 6 (p < 0.001, all) months were significantly higher than 1 month.

• Download figure
• Open in new tab

Figure 1.

Box plot with error bar (standard deviance) of the upper-extremeity Fugl-Meyer motor score (A) and Barthel index (B) within 1 week (n = 12), 1 (n = 16), 3 (n = 13) and 6 (n = 9) months after first-time acute stroke. Post hoc tests with Bonferroni corrections for multiple comparisons of different time points were performed(*p <0.05, **p<0.01, ***p<0.001)

Development of rich club organization

The present research tracked the rich club organizational development of structural network of brain at 4 time points after stroke: within 1 week (Figure 2 A), 1 (Figure 2 B), 3 (Figure 2 C) and 6 months (Figure 2 D). Obvious rich-club organization was found from the structural brain network at different time points. Not each individual had normalized rich-club coefficient that was more than 1 in a wide range of degree except that at 6 months. For each time point, a list was obtained by averaging individual rich club nodes scores ranked in descending order by degree. The top 8 (10%) nodes of the list were selected to represent rich club regions at 4 time points and the other regions were identified as peripheral regions, which could ensure rich club organization established by equal numbers of nodes. The rich club regions at within 1 week included: Frontal_Sup_2_L, Frontal_Sup_2_R, Supp_Motor_Area_L, Supp_Motor_Area_R, Insula_R, Cingulate_Ant_R, Cingulate_Mid_L Cingulate_Mid_R. At 1 month after stroke, Cingulate_Ant_L appeared in rich club regions and Supp_Motor_Area_L disappeared. At 3 months after stroke, Cingulate_Ant_L disappeared in rich club regions and Putamen_R appeared. Putamen_R disappeared and Insula_L first appeared in rich club regions at 6 months after stroke when subjects were functionally recovered.

Degree of 0 < k < 34 was reported in this research according to the highest degree k more than 80% subjects had at all time points. For each degree of k, linear mixed model was used to test the effect of time on normalized rich club coefficient. Normalized rich club coefficient was affected by time when degree was in the range 28 < k < 33 (p<0.05, the logic of the likelihood ratio test). Time showed significant effect on normalized rich club coefficient when k is equal to 20 (, p = 0.007, β_time = 0.010), 21 ( = 6.87, p = 0.008, β_time = 0.011), 22 (, p = 0.008, β_time = 0.013), and 23 (, p = 0.012, β_time = 0.014).

• Download figure
• Open in new tab

Figure 2.

Mean normalized rich club coefficient curves (coloured lines) for the structural brain network at (A) within 1 week, (B) 1, (C) 3 and (D) 6 months after first-time acute stroke. Individual normalized rich-club curves are also shown in grey. Each individual and mean normalized rich club coefficients were larger than 1 in a range of k, suggesting rich club organization existed in all networks.

Development of network topological metrics

To explore the detail development of brain network, topological metrics were examined in regional scale, rich club scale, and large scale respectively. Regional scale means metric related to each node. Rich club scale means mean nodal metrics of rich club and periphery nodes, or mean link metrics of rich, feeder and local connections. Large scale means the mean metrics of nodes of the whole brain network. Topological metrics in three scales were calculated including nodal degree, nodal strength, clustering coefficient, local efficiency, and betweenness centrality through the method referred in Materials and Methods. For each node, linear mixed model was used to test fixed effect of time on nodal scale metrics with age and gender as covariates, subject as random variable. 11 regions of nodal degree changed with time including Frontal_Sup_2_R, Frontal_Inf_Oper_L, Frontal_Sup_Medial_R, OFCmed_L, OFCmed_R, OFCant_L, OFClat_L, Parietal_Sup_L, SupraMarginal_L, Pallidum_R and Heschl_R (Figure 3 A). 14 regions of nodal strength changed with time including Rolandic_Oper_R, Frontal_Med_Orb_R, OFCmed_L, OFCmed_R, OFCant_L, OFCant_R, Insula_L, ParaHippocampal_L, Amygdala_R, Calcarine_R, Cuneus_R, Lingual_L, Lingual_R, Heschl_L (Figure 3 B). 3 regions of local clustering coefficient changed with time including Frontal_Sup_2_R, Supp_Motor_Area_L, Supp_Motor_Area_R and Frontal_Sup_Medial_R (Figure 3 C) 11 regions of local efficiency changed with time including Frontal_Inf_Oper_R, OFCant_L, OFCant_R, Insula_L, Hippocampus_L, ParaHippocampal_L, Lingual_L, Lingual_R, Caudate_L, Thalamus, Temporal_Sup_R (Figure 3 D). 6 regions of nodal betweenness centrality changed with time including Frontal_Sup_Medial_R, OFClat,Cingulate_Post_L, Lingual_R, Occipital_Mid_L, Postcentral_R (Figure 3 E).

• Download figure
• Open in new tab

Figure 3.

Regional brain network measures changes in the process of stroke recovery. A. Nodal degree B. Nodal Strength C. Local clustering coefficient D. Local efficiency E. Nodal betweenness centrality. Only regions with significance were displayed. *p<0.05, **p<0.01,***p<0.001. Red * represent rich club regions.

Linear mixed model was used to test fixed effect of time on large scale and rich club scale metrics with gender and age as covariates and subject as random variable. Eight metrics changes were shown in Table 2. Time showed significant negative effect on mean degree (β_time = −0.151, p = 0.008), mean degree of peripheral nodes (β_time = 0.145, p = 0.007), mean strength of peripheral nodes (β_time = −0.119, p = 0.049), mean density of local connections (β_time = −8.934) and density ratio of local connections (β_time = −0.005, p = 0.029). While time showed positive effect on density ratio of feeder connections(β_time = 0.006, p = 0.008), Communication cost ratio of rich connections (β_time = 0.005, p = 0.020) and ccommunication cost ratio/density ratio of rich connections (β_time = 0.141, p = 0.023).

Correlations between motor behaviour and topological metrics and rich club metrics

Results of correlations between motor behaviour and topological metrics and rich club metrics were showed in Table 3. Communication cost ratio of local connection (β = −72.81, p = 0.018; β = −133.30, p = 0.005; respectively), communication cost ratio/density ratio of local connections (β = −41.21, p = 0.031; β = −64.11, p = 0.037; respectively) were found to correlated with both UE-FM and BI in the process of stroke recovery. In addition, mean clustering coefficient of peripheral nodes (β = −1326.693, p = 0.047), mean length of feeder and local connections (β = −0.354, p = 0.034; β = −0.197, p = 0.016, respectively) were found to be negative related with BI. Mean degree(β_time = 0.006, p = 0.008) and mean degree of peripheral nodes (β_time = 0.006, p = 0.008) were found to be negative with UE-FM.

Motor behaviour prediction

Linear mixed models were used to explore communication metrics of brain network on the change of motor outcome(UE-FM and BI). Data of all points were used in the models. First, the factors correlated to change of motor outcome (subtract the baseline motor outcome from motor outcome at other time points) were selected by comparing the model with the metric and without the metric with each subject as random variable, metric to be tested as fixed variable, and time, age and gender as covariates. Communication cost of rich club (β = −0.007, p = 0.031) and normalized rich club coefficient were found to correlate with the change of UE-FM. Density of local connections(β = −0.173, p = 0.019), cost ratio of feeder and local connections(β = 188.912,p = 0.016; β = −145.202, p = 0.025), communication cost ratio of feeder and local connections(β = 492.186, p = 0.01; β = −337.023, p = 0.018) and normalized rich club coefficients were found to correlate with the change of BI. Results were shown in Table 4. Variables were selected to prediction motor function as Table 5. Communication cost of rich connections and normalized rich club coefficient(k=21) can predict UE-FM change. Density of local connections, cost ratio of feeder connections, communication cost ratio of feeder connections can predict BI change.