Modeling functional network topology following stroke through graph theory: functional reorganization and motor recovery prediction

The study of functional reorganization following stroke has been steadily growing supported by advances in neuroimaging techniques, such as functional magnetic resonance imaging (fMRI). Concomitantly, graph theory has been increasingly employed in neuroscience to model the brain's functional connectivity (FC) and to investigate it in a variety of contexts. The aims of this study were: 1) to investigate the reorganization of network topology in the ipsilesional (IL) and contralesional (CL) hemispheres of stroke patients with (motor stroke group) and without (control stroke group) motor impairment, and 2) to predict motor recovery through the relationship between local topological variations of the functional network and increased motor function. We modeled the brain's FC as a graph using fMRI data, and we characterized its interactions with the following graph metrics: degree, clustering coefficient, characteristic path length, and betweenness centrality (BC). For both patient groups, BC yielded the largest variations between the two analyzed time points, especially in the motor stroke group. This group presented significant correlations (P<0.05) between average BC changes and the improvements in upper-extremity Fugl-Meyer (UE-FM) scores at the primary sensorimotor cortex and the supplementary motor area for the CL hemisphere. These regions participate in processes related to the selection, planning, and execution of movement. Generally, higher increases in average BC over these areas were related to larger improvements in UE-FM assessment. Although the sample was small, these results suggest the possibility of using BC as an indication of brain plasticity mechanisms following stroke.


Introduction
Advances in neuroimaging have enabled the study of the reorganization of brain function, which has been identified as one of the fundamental mechanisms involved in motor recovery following stroke (1,2). Recent neurologic research has emphasized the important role of distributed networks in the brain. Although the brain's structure is changed by focal damage, this change influences the function of distant brain regions (3,4). Focal brain lesions resulting from ischemic stroke may yield selective alterations in functional and structural interconnectivity of other brain circuits that are unrelated to motor function, such as recovery, although the precise biological mechanisms must still be determined (1).
Network randomization may be a final common outcome following stroke damage, resulting from a compensatory but nonoptimized outgrowth of new connections due to an impaired connection pathway (1). Wang et al. (11) showed that motor execution networks in patients were increasingly random over the course of one year of recovery. New axonal outgrowths may be partly responsible for this randomization, but they may not be the only cause. After ictus, other structural and functional changes may also contribute to the continued randomization of network configuration (11). Lee et al. (1) showed that the topological configuration of the network shifted toward a random network 3 months following ictus. The authors also found a relationship between low characteristic path length (which characterizes high network efficiency) in the ipsilesional hemispheric network just after ictus and better recovery 3 months post-stroke (1). However, the cited studies evaluated only patients with moderate to severe motor deficits, while patients with mild motor impairment or without impairment were excluded.
The aim of this work was to investigate the reorganization of network topology in both ipsilesional (IL) and contralesional (CL) hemispheres during functional recovery in two groups of stroke patients: with and without motor impairment. We also attempted to predict motor recovery from the relationship between local topological variations (more specifically, in the premotor cortex, supplementary motor area, and primary sensorimotor cortex) of the functional network and increased motor function.

Subjects and experimental design
A total of 33 patients who had their first-ever stroke were assessed for eligibility. All stroke patients were treated during the subacute phase at the Emergency Department at the University Hospital of the University of Campinas. They were enrolled in this study within the first month for the first evaluation, and they were subsequently reassessed at 3-4 months after stroke, when they were seen for follow-up at the outpatient stroke clinic at the same institution. The inclusion criteria were: 1) ischemic stroke; 2) patients up to 4 weeks from the onset of ictus; and 3) motor deficits of the contralesional upper and/or lower extremities. The exclusion criteria were: 1) patients whose symptoms completely resolved within 24 h; 2) hemorrhagic strokes or other neurological disorders; and 3) any contraindication for MRI.
Twenty patients out of 33 were excluded, and 13 patients participated in this study ( Figure 1). All patients underwent two fMRI scans: one performed within 1 month and another 3 months after stroke onset. Experiments were conducted with the understanding and written consent of each participant, and ethics approval was provided by the University of Campinas Ethics Committee (document number: 694.729; CAAE: 0841814.0.0000. 5404).

Clinical assessment
All patients had stroke and underwent clinical assessment on the same days as the MRI scans. All patients were assessed using the modified Rankin (12) score and the Barthel index (13) to evaluate functionality and daily activities, that is, the ability to carry out everyday tasks. Patients were assigned to one of the two groups depending on their first clinical assessment: with (motor stroke group) or without (control stroke group) motor impairment. Patients without impaired functionality had a Rankin score of 0 or 1 and a Barthel score of 100. Patients with impaired functionality had Rankin scores X1 and Barthel scores o100 and were also evaluated using the Fugl-Meyer assessment for motor function of the upper and lower extremities (14). The Fugl-Meyer assessment evaluates sensitivity, speed, coordination, and motor function and is scored between 0 and 2 (0: unable, 1: partly able; 2: fully able to complete movement) with a total score range of 0-66 for upper extremities (UE-FM) and 0-34 for lower extremities (LM-FM) (14). To measure stroke severity, we used the National Institutes of Health Stroke Scale (NIHSS) (15).
We preprocessed the data with the UF 2 C toolbox (16), a MATLAB suite (mathworks.com) for MRI data preprocessing specifically focused on functional connectivity analysis. The preprocessing included functional image realignment, functional-structural images coregistration, and tissues segmentation, normalization, and smoothing. The images origins were manually set at the anterior commissure for all subjects to minimize distortions resulting from the coregistration and normalization steps. In addition, subjects who had more than 30 discarded volumes due to movement constraints were excluded from this study (step known as ''scrubbing'', see Power et al. (17) for a review).
Finally, the BOLD time series were pre-whitened following the approach suggested by Santosa et al. (18) to remove the signal autocorrelation. This enables the calculation of Pearson correlation coefficients for estimating the functional networks while avoiding inflated correlation values, yielding more reliable and meaningful results (19) and, hence, a more accurate graph representation of the brain's functional connectivity (FC). All images were flipped to constrain the lesion's location to the right brain hemisphere (i.e., all patients' lesions were constrained to be on the positive MNI x-coordinates by simply inverting the signal of the voxels along the x-axis when the lesion was located on the left hemisphere).

Data analysis and functional networks
Functional networks were modeled as a graph, a mathematical tool composed of a set of nodes (or vertices), which in this case represented brain regions, and links (or edges), structures that dictate interactions between the nodes. Each node was defined according to the functional atlas proposed by Power et al. (20), and the strength of the links between nodes was calculated using the Pearson correlation, yielding a 264 Â 264 connectivity matrix. Each node was centered at the coordinates suggested by the functional atlas, encompassing all of its nearest neighbors within a three-voxel radius. Hence, each node's time series consisted of the average signal of the aforementioned voxel neighborhood.
To further avoid spurious correlations, following common trends for thresholding the connectivity matrix (21), only the 20% strongest connections were maintained in the graph's adjacency matrix (A). This approach also ensured that the graphs of all subjects had the same total number of functional links. Therefore, each entry a ij of A could take the value of either '1' or '0', indicating whether there was a connection between nodes i and j; that is: The functional connectivity analysis consisted of assessing how specific metrics varied across fMRI acquisitions for the two patient groups. The chosen metrics were degree, clustering coefficient, characteristic path length, and betweenness centrality. The degree and the clustering coefficient were selected for their easy-tointerpret yet powerful meaning. The characteristic path length was chosen to provide complementary information regarding the efficiency of connections within the network, and the betweenness centrality was selected to provide a notion of the importance of a node in the sense of information flow within the network. A brief description of these metrics is provided below, as they will be needed for interpreting the results. For further details, the interested reader may refer to previous studies (22,23).
The degree for a node i (D i ) represents its number of connections, that is (22): where N is the total number of nodes in the network (in this case, N=264). The higher the degree of a node, the greater the number of nodes it connects to. The clustering coefficient (CC) provides an idea of clustering of neighboring nodes: given any three nodes i, j, and k, if i and j are connected, as well as j and k, the CC can be thought of as the probability of i and k also being connected. Mathematically, the CC for a node i can be calculated as (22): The shortest path length represents the number of links that compose the shortest path between nodes i and j (l ij ). If there is no connection between i and j, then l ij ¼ N.
The characteristic path length of node i, then, is simply the average value of l ij , that is (22): The BC is a metric that describes the importance of a node regarding how significantly it acts as a bridge between two other nodes within the network. For node i, this metric can be calculated as (22): in which l jk represents the number of shortest paths that go from j to k, and l jk (i) represents those that specifically pass through node i.
Finally, changes in metrics between the two fMRI scans were calculated as follows: in which M refers to one of the aforementioned metrics, and indices 1 and 2 represent M's values on the first and second scans, respectively. Therefore, DM indicates a relative value for the change in this metric regarding its magnitude in the first fMRI evaluation. Metric variations were attributed to standardized anatomical regions of the AAL atlas by comparing its coordinates to those of the functional atlas used in this study. The regions of interest found in our analyses, which are related to sites where DM variations were more prominent, are reported in Table 1.
Finally, the last analysis involved seeking correlations between variations in the average BC and the UE-FM scale on the primary sensorimotor cortex, premotor area, and supplementary motor area.

Results
Of the 13 participating stroke patients, six presented impaired motor function (motor stroke group, mean age 63±9 years, range 53-80), and seven did not present impaired motor function (control stroke group, mean age 58 ± 9 years, range 49-72). Table 2 presents the baseline characteristics for all subjects. There were no statistically significant differences between the groups regarding age, sex, time post-stroke, or stroke hemisphere. The patients in the motor stroke group had statistically worse scores on the modified Rankin score and Barthel index (see Table 2).

Group average results
Graphs of metric variations are displayed in color maps and in relative (percentage) scales. Each metric will be discussed separately, as each one carries a specific and unique meaning regarding its interpretation within the functional network. Figures 2 and 3 show variation maps for the studied metric between both fMRI acquisitions for patients of the motor stroke group and control stroke group, respectively. The behavior shown in the figures was estimated by taking the average connectivity matrix for each patient group. Areas closer to yellow indicate increases, while those closer to blue indicate decreases. In these maps, the white areas indicate values considered too small to be displayed, that is, the white regions in Figures 2 and 3 indicate unplotted values. Only variations that exceeded half the average DM across the whole brain are shown. Finally, we reiterate that the images of some patients were flipped to restrict the lesions to the right hemisphere in all subjects.
To further investigate the patterns of the regions in Figures 2 and 3, we computed average values for the metrics of the cortical and subcortical areas with the largest variations. The results are shown in Table 3, with motor stroke and control stroke groups separated by graph metric and hemispheres, as well as for increase or decrease of the corresponding metric, which are indicated as upwards and downwards arrows, respectively. Note that some regions presented both increases and decreases, and thus both arrows are displayed. We interpreted the presence of both increasing and decreasing trends as changes belonging to distinct portions (i.e., distinct graph nodes) of the same anatomical region.

Individual BC results
Based on the results in Figures 2 and 3 indicating that BC was the metric with the largest changes in fMRI acquisitions, we further explored its variations for each subject. These results are shown in Figure 4 for all participants in both the motor impairment (A) and control stroke (B) groups.
In addition, we also investigated whether this metric correlated with the UE-FM by studying how the average BC change for a given subject was related to the Fugl-Meyer score change ( Figure 5). For this analysis, three areas of interest were explored: the primary sensorimotor cortex, the supplementary motor area, and the premotor cortex for both the CL (upper panels) and IL (lower panels) hemispheres. The correlation strength between the two variables, i.e., the oDBC4 and the corresponding changes in UE-FM (indicated by DUE-FM), are shown as r, that indicates the value of the correlation coefficient. When the p corresponding to that rho is less than 0.05, the correlation is significant, with an asterisk indicating statistical significance (Po0.05).

Discussion
Motor stroke group Group average results. Several studies have reported alterations in brain topological organization and disruptions in functional connections in patients following stroke when compared to healthy controls. As expected, there is a decrease in global efficiency, indicating a reduced capacity for information transfer across the entire stroke brain (9,24). In this work, to better understand brain changes after stroke, we compared functional network changes between groups of stroke patients with different stroke outcomes, namely, with and without motor impairment. Interestingly, robust changes in global integration, including alterations in strength, clustering coefficient, characteristic path length, and betweenness centrality, were identified in both groups, as well as a correlation of these network changes with clinical variables that assess motor impairment (for the motor stroke group).
As widespread brain regions and extensive networks may be damaged in stroke patients, studies investigating the whole network of functionally interacting brain regions may be more valuable for a better understanding of the pathological mechanisms of stroke than studies investigating local connections (9). In this context, we found several connectivity alterations over time, not only in the lesion area but also in the frontal and temporal regions, the parietal gyrus, and the basal ganglia, in both hemispheres. Although it is tempting to attribute these changes to brain reorganization as a mechanism to suppress further impairments due to the stroke, it is important to be cautious with these conjectures, given that the analyzed patient sample was small, and the patients had a great variability regarding the damaged brain areas.
The largest increases in metrics on the IL hemisphere ( Figure 2A) were observed in the frontal and temporal areas. In the CL hemisphere, the largest increases were observed in the parietal cortex (Brodmann 40), and the largest decreases occurred in the superior temporal cortex. Some gyri were common to both hemispheres, presenting the highest variations (Table 3), but the degree tended to increase more prominently in the CL hemisphere. In other words, the number of connections tended to increase more significantly contralateral to the lesion. An increased degree (i.e., increase in the number of functional connections) for specific regions implies that such regions display more synchronicity of their recorded fMRI time series at a global level (that is, considering every other region of interest in the brain). For the clustering coefficient ( Figure 2B), the CL hemisphere tended to present the most significant cortical increases, especially in areas such as the superior and inferior temporal gyri and the postcentral gyrus. Ipsilateral to the lesion, most displayed locations remained at approximately the same values for this metric, with the largest (yet mild) clustering coefficient increases being observed on the precentral and medial frontal gyri (Table 3). In addition, altered areas for the CC were similar to those involving the degree (Figure 2A and B); relative variations, however, were greater for the latter metric (CC highest changes were approximately 5%, while maximum degree changes were approximately 10%). Since the CC considers functional interactions at a more local scale, as it accounts for elements of the adjacency matrix that are nearest neighbors to each another [see Equation 3], when analyzed simultaneously with degree, it can further indicate whether, to any extent, the variations in the number of functional connections of a given region of interest were spatially restricted to it. In this case, a qualitative visual inspection of Figure 2A and B suggests that this seemed to be the case for the motor stroke group.
The characteristic path length ( Figure 2C) followed the same principle as the other metrics, with increased values in the CL hemisphere. Increases in the CPL imply that the information must travel through a larger number of functional links from origin to destination. Hence, regions with increased (or decreased) CPLs can be regarded as requiring a higher (or smaller) number of functional connections to transfer or receive information, indicating a loss (or gain) of information transfer efficiency. To summarize, when concomitantly analyzing the findings for the three aforementioned metrics, the CL (or IL) increased (or decreased) for the degree, and the CC indicated a higher (or lower) number of connections and efficiency in local communication of a node, resulting in a smaller (or larger) CPL. In other words, there seems to be a tendency for the CL hemisphere to become more efficient in information transfer, perhaps due to compensatory mechanisms following a stroke in the other hemisphere.
Greater connectivity and efficiency of the CL hemisphere may occur to compensate for loss of connectivity in the ischemic region. Functional neuroimaging studies suggest that activity in the sensorimotor network or ipsilesional motor cortex, is most abnormal in the early phase after hemiparetic stroke, and that motor recovery is related to normalization of activity (25). In a previous cross-sectional study, we had already demonstrated increased connectivity on the CL hemisphere in stroke patients with impaired functionality (6). It is possible that the normalization of brain activity occurs later. Lee et al. (1) suggest that the damaged brain enters a chronic stage and further reorganization might not occur because most network reorganization usually happens during the first 3 months after stroke.
Finally, the BC ( Figure 2D) displayed the most dramatic variations, with increases of more than 50%. Additionally, no significant decreases were observed. The regions with high values of BC are considered network hubs, and group differences in BC of nodes reflect effects of the disease on the global roles of regions in the network (26,27). Similar to the works by Zhang et al. (9) and Yin et al. (24), we verified an increased nodal BC in the CL superior frontal gyri (identified as supplementary motor area) and the CL inferior parietal lobule, which may participate in processes related to movement selection and movement planning. These areas may be viewed as reflecting the selective neuroplastic recruitment of the unaffected motor network to compensate for the damage induced by the lesion. Furthermore, increases occurred in subcortical areas, such as the basal ganglia and thalamus, structures involved in known motor circuits.
Individual BC results. When analyzing BC variations individually, it is noticeable how, despite existing common trends among all subjects, there were many particularities. However, overall, frontal, and parietal areas were common sites of the largest alterations (although the spatial extent of the involved sites is very subject-dependent). BC can be an important parameter to identify changes in functional reorganization since it increased in the frontal and parietal regions, which are responsible for motor planning and integration of sensory input and motor output. This reflects the clinical evolution of patients Table 3. Increase (m) or decrease (k) of graph metrics in the motor stroke and control stroke groups.

Graph metric
Control stroke group Motor stroke group  Table 1. with motor deficits over time, since there was an improvement in functionality and greater independence in daily activities, shown by the Rankin and Barthel scales and by the Fugl-Meyer assessment. However, the sample was too small for further inferences about statistical significance.
Furthermore, the results in Figure 5 showed that the average changes in the BC and UE-FM assessments were significantly correlated for the primary sensorimotor cortex and the supplementary motor area for the CL hemisphere. On the other hand, no significant correlation was found for the same areas on the affected hemisphere, nor for both hemispheres for the premotor cortex. This means that both the increase in importance of the primary sensorimotor area and the supplementary motor area regarding information flow through the brainas reflected by the BC metricwere actually linked to clinical improvement of patients over time. This finding suggested that the BC in these areas might be an indicator of neural plasticity related to motor recovery. Moreover, the fact that these changes were located on the CL hemisphere may suggest some type of compensatory mechanism by the unaffected areas. Nevertheless, further analyses are still necessary at this point, as our small number of subjects limits the generalizability of our results.

Control stroke group
Group average results. Most change patterns were similar to those found for the motor-impaired patients. The degree and clustering coefficient presented the most significant increases in the CL hemisphere. As with the motor impairment patients, most sites with increased degree (such as the frontal and temporal gyrus) also roughly exhibited larger clustering coefficients. The metric's changes were restricted to a small range, from 2 to 5% ( Figure 3A and B).
As observed for the motor stroke group, these patients displayed decreases (or increases) in the CPL in the proximities of sites that exhibited increases (or decreases) in the degree and CC ( Figure 3C). However, locations encompassing structures such as the superior frontal, superior temporal, middle frontal, precentral, and postcentral gyri showed considerable decreases (although these variations were small compared to the other metric, of approximately 1%).
The BC alterations were also the most considerable in this group, reaching values up to 40%. The main areas of increase in the CL included the superior frontal, superior temporal, and inferior temporal gyri, whereas IL increases mainly occurred in the middle frontal, inferior frontal, precentral, and postcentral gyri. Unlike the motor-impaired patients, the average for the control group also displayed considerable BC increases in the proximities of the cerebellum ( Figure 3D).
Indeed, several studies have shown an increase in cerebellar activity in stroke patients compared with healthy individuals (11,24). In the present study, patients without motor deficits also showed increased cerebellar connectivity over time. Most likely, the greater motor performance of these patients compared with patients with motor deficits may have been influenced by activation of the cerebellum, since this region interferes with coordination and motor learning.
In a previous work, we found higher connectivity in some non-motor networks in healthy controls than in stroke patients with impaired functionality (6). Now, patients in the control group had more regions with variations in connectivity over time than the motor deficit group, either increasing or decreasing. One hypothesis is that these modifications and the alterations of other areas occurred because it is not only the motor circuit that can help restore the compromised side. However, it is still unclear how all these changes, even with reduced connectivity, can interfere with the process of brain reorganization and functional recovery.
Individual BC results. Overall, patients in the control stroke group presented a tendency to display the most significant increases in the CL hemisphere ( Figure 4B; see, e.g., Subjects 1, 2, 3, 4, and 5). Moreover, occipital areas and sites close to the cerebellum were more prominent regarding BC changes for this group than for the motor stroke group ( Figure 4B). Subject 7, in particular, presented the highest and most noticeable variations.
Although some studies indicate that patients with good functional rehabilitation outcome usually have activation in the perilesional area during functional activities (28), other findings indicate CL activation, which is more in agreement with our present findings. Indeed, Calautti et al. (29), Fujii et al. (30), and Almeida et al. (6) suggest that increased activity in CL sensorimotor cortices is an available mechanism for compensating, at least partially, Figure 5. Relationship between upper-extremity Fugl-Meyer (UE-FM) scale changes and average betweenness centrality (BC) variation for the motor stroke group. In this plot, each point corresponds to the oDBC4 and DUE-FM values of a subject. Plots in the upper and lower panels correspond to the contralesional (CL) and ipsilesional (IL) hemispheres, respectively. Each column designates oDBC4 values gathered at specific areas. The correlation strength between the two variables, i.e., the oDBC4 and the corresponding changes in UE-FM (indicated by DUE-FM), are shown as r, that indicates the value of the correlation coefficient. When the P corresponding to that rho is less than 0.05, the correlation is significant. The correlation strength between the two variables is shown with an asterisk (*) when significant (Po0.05).
for stroke-induced motor impairments. Also, Xu et al. (31) and Swayne et al. (32) reported that while it is unclear how the CL activity may influence motor recovery, the reorganization of residual tissue to re-enable motor function likely depends on some degree of intracortical disinhibition to allow access to additional networks.

Conclusions
In this study, we investigated how graph metrics vary over time between two distinct fMRI acquisitions for patients with different types of impairment following stroke. We found that most metrics displayed only slight changes in all patients with stroke, with the exception of betweenness centrality, which showed changes of up to 50% for motor-impaired patients and 40% for control patients.
On further investigation of BC alterations, we found significant correlations between average change in BC and alterations in the UE-FM score for the CL hemisphere in the primary sensorimotor cortex and in the supplementary motor area for the motor impaired group, up to 3-4 months after stroke. However, additional investigations with a larger number of subjects are necessary to verify whether this relationship can be generalized. Therefore, at this point, we can only speculate that BC may reflect aspects underlying brain plasticity mechanisms after stroke. These results should be explored in future studies in the field.