Brain anomaly networks uncover heterogeneous functional reorganization patterns after stroke.

Stroke has a large physical, psychological, and financial burden on patients, their families, and society. Based on functional networks (FNs) constructed from resting state fMRI data, network connectivity after stroke is commonly conjectured to be more randomly reconfigured. We find that this hypothesis depends on the severity of stroke. Head movement-corrected, resting-state fMRI data were acquired from 32 patients after stroke, and 37 healthy volunteers. We constructed anomaly FNs, which combine time series information of a patient with the healthy control group. We propose data-driven techniques to automatically identify regions of interest that are stroke relevant. Graph analysis based on anomaly FNs suggests consistently that strong connections in healthy controls are broken down specifically and characteristically for brain areas that are related to sensorimotor functions and frontoparietal control systems, but new links in stroke patients are rebuilt randomly from all possible areas. Entropic measures of complexity are proposed for characterizing the functional connectivity reorganization patterns, which are correlated with hand and wrist function assessments of stroke patients and show high potential for clinical use.

Introduction

Stroke is the second largest cause of death worldwide, and it has a high burden on patients, their families, and health-care systems (Valery et al., 2014). Non-invasive techniques, such as functional MRI, play a crucial role in measuring the abnormal activity of the brain, which contains much information for prognostic evaluation and treatment planning. Given fMRI data that are collected during resting state or various cognitive tasks, complex network approaches have been able to characterize inter-relationships between different brain areas (Sporns et al., 2005). A human brain functional network (FN) can be obtained from a thresholded correlation matrix, which is derived by computing all pairwise correlations of time-resolved blood‑oxygen-level dependent (BOLD) signals from fMRI. The resulting brain networks have specific characteristics as revealed by graph analysis (Bullmore and Sporns, 2009). For instance, they consist of a number of highly connected hub nodes, show high clustering in combination with short path length yielding so-called small-worldliness, and have communities forming a hierarchical modular organization. These ubiquitous network features of healthy subjects have been observed widely in functional imaging and suggest high efficiency of human brains in processing information (van den Martijn et al., 2010; Zamora-Lopez et al., 2010; Danielle, 2017).

Deviations of FNs from optimal topology have been related to cognitive and clinical symptoms in a wide range of psychiatric or neurological disorders (Cornelis, 2014; Fornito et al., 2015; Liang et al., 2010), which result in re-organization of functional connectivity among brain regions (Aerts et al., 2016). For instance, it was shown that long-range paths and modularity are significantly reduced in Alzheimer's disease and schizophrenia (Betty et al., 2013; Aaron et al., 2010; Zac et al., 2015). Moreover, hub disruptions are conjectured to be fundamental in various brain pathologies, because hubs are vulnerable to targeted attacks (Cornelis, 2014; Achard et al., 2012; Nicolas et al., 2014). The hypothesis of hub vulnerability states that an acute severe injury could cause a significant functional damage of hub regions, resulting in rapid network fragmentation. It remains unknown, however, to what degree this random re-organization might take place, or, how this hypothesis could serve as a biomarker to collect information for clinical use.

Methodologically, it is very difficult to perform group-wide analyses across patient groups, since each patient may react differently to brain disorders, for instance, concerning the site of lesions and severity of damage. Hence, network properties based on small-worldliness and modularity structures were not found to sufficiently differentiate between healthy and patient groups (Achard et al., 2012). Indeed, we find that topological properties of brain FNs are largely conserved in our patient group. How, then, can particular features of a patient be captured in comparison to the group of healthy controls? This is a challenging question of modern graph analysis, which requires the consideration of large anatomical and functional variability among subjects.

In stroke disorders, one potential approach is to study the variations of connectivity patterns for a series of particular regions of interest (ROIs). One popular choice for ROIs is based on the Automated Anatomical Labeling (AAL) atlas (Tzourio-Mazoyer et al., 2002). This is a well-established approach, but, depends critically on one's prior knowledge of the pathological information which guides the choice of ROIs. The boundaries of ROIs can dramatically influence the conclusions of such studies (Hayasaka and Laurienti, 2010; Matthew et al., 2013). To date, no conclusive insights have been obtained on the relationship between the random reorganization hypothesis and pathological impairment.

There are three main objectives in this study. First, we propose a framework for constructing spatial anomaly networks (SAN) in stroke patients. This approach combines time series information of each individual patient with the healthy control population. Specifically, given the BOLD signal of one voxel from both groups (one patient vs. the healthy control ensemble) at the same position, we compute the associated anomaly time series which characterizes the abnormal variability (deviation) of the voxel in comparison to the group of healthy volunteers (see Materials and Methods). The underlying assumption is that there is a representative brain network for the matched healthy control group. Second, to overcome the subjective criteria required to choose ROIs, we propose a data-driven algorithm to identify several fundamental ROIs that are most relevant to stroke pathology. Third, we assess the re-organization mechanism, showing how strong functional connections of ROIs of the healthy group are destroyed in stroke patients and how strong links in patients are newly, and randomly, created. Finally, we propose entropy as a complexity measure for characterizing connectivity re-organization, which shows a high correlation to clinical assessments of severe motor impairments.

Materials and methods

Experimental design and data acquisition

We collected data from 32 sub-cortical stroke patients (after head motion correction) with left motor pathway damage and 37 age-, gender-, and handedness-matched healthy controls. Both groups had corrected head movement. The healthy volunteers were denoted by H = [1, 37], while the patients were indexed by P = [38, 69]. Inclusion criteria of patients were as follows: (1) first-onset stroke, (2) pure motor deficits, (3) right-handedness, (4) sufficient cognitive abilities (Mini-Mental State Examination, MMSE ≥27), (5) examination time more than 3 months from stroke onset. See Supplementary Materials (SM-I) for a full description on data acquisition and inclusion/exclusion criteria. The detailed clinical and demographic information of the patients, including the location of lesions and lesion volumes, is discussed in the SM (see SM-I, Table S1, S2, Figs. S1, S2).

All 32 post-stroke patients were evaluated by the Fugl-Meyer Assessment (FMA) (Gladstone et al., 2002). FMA is a well-designed, pragmatic and efficient test that has been widely applied clinically and in research to determine disease severity, describe motor recovery, and plan and assess treatment. A smaller FMA value means that the patient is more severely affected by stroke. In the clinic, many stroke patients that appear to have similar locations and extent of lesions may actually have extremely different outcomes in hand functions. Some have regained certain level of practical abilities in hand function after the stroke, whereas others have lost all functional capacity of their hand. Therefore, following our previous studies (Yin et al., 2012; Yin et al., 2014), we focus on FMA scores that test hand and wrist function.

The fMRI data acquisition steps are described in the SM (see SM-I). The resting state fMRI from all subjects were acquired on a Siemens Trio 3.0 Tesla MRI scanner (Siemens, Erlangen, Germany) at the Shanghai Key Laboratory of Magnetic Resonance, East China Normal University (ECNU). The protocol for this study was approved by the Institutional Ethics Committee of ECNU (Shanghai, China), and all participants or their guardians signed informed written consent. Resting state fMRI data of the whole brain were acquired using an echo-planar imaging sequence (EPI). Both T1 and T2 weighted images were also collected (see SM-I for all parameter settings for EPI, T1 and T2 acquisition). During resting-state fMRI data acquisition, the participants were instructed to remain awake, relaxed with their eyes closed, and motionless without thinking about anything in particular. Each scan lasted for 8 min and 6 s; however, the first 6 s was consumed by a dummy scan. In total, we collected 240 image volumes for a subject.

Data preprocessing

Preprocessing of the resting-state fMRI data was performed using Statistical Parametric Mapping (SPM8, http://www.fil.ion.ucl.ac.uk/spm). We discarded the first 10 volumes of the dataset for each participant to allow for magnetization equilibrium, leaving 230 volumes for further analysis (see SM-II).

In constructing brain FNs, one traditional way to acquire time series data is based on parceling the functional images into cortical and sub-cortical regions of interest (ROIs), for instance, in the so-called AAL template of 90 ROIs (Tzourio-Mazoyer et al., 2002). Then, the representative time series of each ROI is extracted by averaging the BOLD time series across all voxels within that region. This coarse graining method has drawbacks, since each ROI defined by the AAL template occupies in general a different brain area (Hayasaka and Laurienti, 2010). In order to account for both noise effects and the computational power available to an ordinary PC, we subdivided each AAL region into 40 subregions and extracted one representative time series by averaging over the particular subregion, which yielded 3600 subregions (subR) for the entire brain (see SM-III). Our up-sampling procedure does not cross the relevant anatomical boundaries specified by the original AAL atlas, as often required (Hermundstad et al., 2013). Note that consistent results were obtained if we performed similar analyses for the whole brain at a voxel-wise level (see SM-VII).

Anomaly time series

Considering the healthy controls as a reference group, we obtained anomaly time series for each subR of each subject by the following steps: (1) Let us assume that the healthy controls behave in a similar manner. Hence, we first quantified the averaged characteristics and variance for the healthy group. Denoting the time series of i-th subR as x(i, t), i ∈ [1, 3600], t = [1,230], and H ∈ [1, 37] as the index for a healthy subject. We computed 〈x〉(i, t), where 〈⋅〉 was the ensemble average over all 37 healthy subjects. Respectively, the variance over the healthy group was denoted by the standard deviation σ(i, t). (2) For the time series of each subR of a patient x(i, t), P ∈ [38, 69], we extracted the anomaly time series y(i, t) aswhich captures the deviations of the BOLD signal of the i-th subR from the healthy group. Note that σ(i, t) characterizes the variance of the i-th subR across the ensemble of the healthy subjects. Further improvement could consider the variance across all brain regions as well.

Spatial abnormal networks

For each patient P, we constructed a network from the anomalies. First, we obtained the correlation matrix ℛ by calculating the Pearson correlation coefficient between any pair of subRs i and j. Then, we defined that two subRs are functionally connected if their temporal correlation exceeds a predetermined value (p-value <0.05). More specifically, from the correlation matrix ℛ, a threshold value was chosen such that the resulted functional network had a specified density of links ρ (see SM-IV A).

Identifying fundamental ROIs

There were four steps in the following data-driven procedure, which was proposed to identify subRs that behave substantially different with respect to the healthy ensemble. (1) Given time series of i-subR x(i, t) and its associated anomaly y(i, t), the sensitivity of i-subR to the brain disorder was characterized by the linear correlation coefficient. More specifically, we computed R(y(i, t), 〈x〉(i, t)), i ∈ [1, 3600]. The smaller the correlations, the larger the deviations of the subRs from the healthy conditions. Our results do not change if we replace 〈x〉(i, t) with x(i, t). (2) We sorted R in an ascending order and chose the bottom 5% of most different subRs as the representative for the patient, which yielded 3600×5% = 180 subRs. (3) We repeated the above two steps for all patient subjects. A subR was considered to be substantially affected by stroke if it appeared consistently (frequently) in 2/3 of the patient population. Only 15 (out of 3600) subRs were identified as meeting these criteria, which are therefore considered to be most-relevant to stroke. (4) We denoted 15 subRs as fundamental by ROI(i, t), i ∈ [1, 15]. These ROIs are preserved when using different thresholds, for instance, 3/4 or 1/2 of the patient population size (see SM-V for discussions on the effects of different thresholds).

Results

Global network properties are not sufficient to characterize the difference between patients and healthy subjects

It has been demonstrated that the degree distribution, P(k), of brain FNs for healthy subjects often follows power-laws, (Eguíluz et al., 2005), or exponentially truncated power-laws (Bullmore and Sporns, 2009). In patients, one might expect some sort of deviation from a power-law distribution. Our numerical results, however, do not show evidence in support of this hypothesis, since P(k) does not show a clear qualitative difference between healthy and patient groups (see Fig. S3, SM-IV B).

Further, we considered various network structural measures (Bullmore and Sporns, 2009; Newman, 2003) as potential candidates for discriminatory statistics: clustering coefficient ; average shortest path length ℒ; modularity ; transitivity ; assortativity ℛ; and entropy of degree randomness (Fig. S4, SM-IV C). By two-sample t-tests based on global network measures, we concluded that no significant differences were captured, although pronounced clinical differences were observed between these two groups. Instead, all results suggested consistently that both healthy and patient groups share some general network properties, such as small-worldliness and modularity. Furthermore, we establish the statistical significance on the network level through the connection rewiring method (Newman, 2003) (see Fig. S4, SM-IV D for statistical significance). Thus, these general network measures do not indicate which brain areas play a decisive role in the reorganization of functional connections. This finding calls for more sensitive network markers of brain disorders, focusing on local properties.

A data-driven algorithm identifies disorder-relevant representative ROIs

One practical solution for understanding the reorganization mechanisms is to study the connectivity changes of particular ROIs, for instance, the bilateral precentral gyrus and other movement-related brain regions, since these regions are believed to be crucially involved in the development of stroke disorders. One drawback of this approach is that it relies heavily on prior knowledge. Here we propose a data-driven approach for identifying critical ROIs.

Taking the healthy controls as a reference group, we propose to compute anomaly time series from BOLD for each sampled subR (see Materials and Methods). By anomaly we mean a departure from a reference value which is provided by the control group. In this work, we considered the ensemble average as the reference value, and examined the variance of the reference group. A positive anomaly indicates that the observed BOLD signal has higher amplitudes than the reference value, while a negative anomaly indicates that the observed signal has smaller amplitudes than the reference value. The computation of anomalies for all subRs provides a global survey of the abnormalities of a patient compared to the averaged reference group.

Every patient responds differently to a stroke. Consequently, strong variations within the patient group prevented us from consistently selecting subRs to reflect the stroke-relevant features. We considered a subR to be affected by stroke if it behaves abnormally over 2/3 of the patient population, as shown in Fig. 1. Discussions on the effects of different thresholds such as 3/4 or 1/2 of the patient population are included in the SM (see SM-V, Tables S4 and S5). Remarkably, only 15 (out of 3600) subRs were identified as meeting this criterion. These affected subRs represent 11 ROIs in terms of the AAL template. All 15 subRs are visualized in Fig. 1 (see SM-V, Table S3). The following brain regions were identified as crucial in stroke: bilateral precentral gyrus (PreCG), paracentral lobule (PCL), supplementary motor area (SMA), precuneus (PCUN), the right postcentral gyrus (PoCG), the left dorsolateral superior frontal gyrus (SFGdor) and the right superior parietal gyrus (SPG). Some notes on the roles of these brain regions in stroke patients are included in the Discussion.

Fig. 1

Identifying most-stroke-relevant fundamental subRs. The shaded area corresponds to subRs, as highlighted by red stars in the upper panel, that behave abnormally over 2/3 of the patient population size. Surface illustrations for these 15 affected subRs.

Strong connections are broken down specifically, but new links are built randomly

Functional networks of the patient group largely preserve global network properties, that is, the existence of hub nodes, small-worldliness and modularity. Global network measures were not able to disclose the hidden connectivity reorganization (see SM-IV, Fig. S4). Therefore, we analyzed rewiring patterns of strong functional connectivity to each subR (ROI) as identified above. We focused on quantifying strong connections, since they are more likely to be linked to hubs in the resulting network. In particular, there are two crucial questions that need to be addressed: (i) How are strong connections in the healthy controls broken down? That is, which particular strong connections disappear in patients? (ii) Where do newly created strong connections arise in patients? We show subject-specific properties in comparison to the healthy control group in Fig. 2, Fig. 3.

Fig. 2

Comparison of strong connections of the left precentral gyrus (the first ROI, PreCG.L) between patient subject P = 53 and averaged healthy controls. (A) scatter plot of the correlation coefficients (r vs. 〈r〉). The difference sets D, E are denoted by red or blue color respectively, while the common set C is highlighted by green. The non-vanishing set C suggests that the strongest links (about 35% of top 100 links) are preserved. (B) Pie plot for the set A of the healthy controls. The strong connections of the PreCG.L are the following: 17% of the connections are within the PreCG.L region (denoted as ‘AAL1’); 17% are inter-hemispheric connections to PreCG.R region (‘AAL2’); 28% are to the left postcentral gyrus (PoCG.L) region (‘AAL57’); 22% are to the right postcentral gyrus (PoCG.R) region (‘AAL58’); 16% spread over the other 6 brain regions (‘AALX=0.16/6’). Strong connections are concentrated on some particular brain areas for the healthy volunteers, yielding a smaller entropy value ℋ〈(1, A) = 2.65. (C) Pie plot for the set B of the patient: about 57% of the strong connections of PreCG.L are distributed across 24 different brain areas. The inter-hemispheric connections to the right precentral gyrus area (PreCG.R) have been greatly reduced from 17% to negligible (B, C). The more homogeneous connectivity distribution in the patient leads to a larger entropy value ℋ53(1, B) = 4.27. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3

Illustrations of the difference sets D and E. (A, B) the healthy control, and (C, D) the patient (P = 53). (A, B) The intermediate strong connections of the PreCG.L area are destroyed as follows: 15.4% are inside the PreCG.L region (denoted as ‘AAL1’); 20% are inter-hemispheric connections to the PreCG.R region (‘AAL2’); 18.5% are to the left postcentral gyrus (PoCG.L) region (‘AAL57’); 24.6% are to the right postcentral gyrus (PoCG.R) region (‘AAL58’); the remaining 22% are distributed across the other 6 brain regions (‘AALX=0.22/6’). One obtains a smaller entropy value ℋ53(1, D) = 2.81. (C, D) In a patient, about 65% of newly created strong connections are homogeneously distributed across 22 brain regions, yielding a large entropy value ℋ53(1, E) = 4.40.

Excluding self-correlation, for all Pearson correlation coefficients r, j ∈ [1, 3599] of the i-th subR, we compared each patient to the healthy group, namely, r versus 〈r〉, where P ∈ [38, 69]. We illustrate the approach with the connectivity patterns of the left precentral gyrus (PreCG.L as represented by the first subR) of patient subject P = 53 as shown in Fig. 2A. From the scatter plot (Fig. 2A), we find that the reorganization of hubs is captured by top correlation values, focusing on the 100 largest values. Specifically, we performed the following partitions of the correlation space to uncover the connectivity reorganization patterns of the i-th subR:

A: set of the top 100 strongest correlation coefficients of the healthy group,

B: set of the top 100 strongest correlation coefficients of a stroke patient,

C = A ∩ B: common set of strong connections shared by both healthy control and a patient,

D = A − B: difference set of strong connections that are in set A of healthy controls, but are rewired in set B of a patient.

E = B − A: difference set of strong connections that are in set B of a patient, but are not in set A of the healthy controls.

We find that the number of subRs in set C is substantially positive across all patient subjects, indicating that the strongest links are not destroyed. The set D quantifies how strong connections of healthy subjects are broken down in a patient, in contrast, the set E characterizes the newly built connections of a patient. For each correlation partition set, we addressed the connectivity heterogeneity over 90 AAL brain areas by computing the entropy (Eq. (2)).

Fig. 4 shows the reorganization behavior of the top 100 strong connections to PreCG.L, which is averaged over the patient group: the strong connections of the healthy group are destroyed in a specific manner from a limited number of brain areas, however, the newly created strong connections are more homogeneously distributed across all brain areas. Consistent results are obtained for each subR (see SM-VI, Fig. S9-S23). Furthermore, inter-hemispheric connections (PoCG.R to the seed ROI of PreCG.L) appear to break down significantly, which confirms recently reported results (Anne and Grefkes, 2013).

Fig. 4

Distributions of connections associated with the representative subR of PreCG.L. across 90 AAL brain areas. (A) 100 strong connections in the healthy group. (B) 100 strong connections in the patient group. (C) strong connections of healthy controls that are broken down in patients for a limited number of brain areas. (D) strong connections that are newly built in patients across all brain areas. These results are further quantitatively characterized in Fig. 1, Fig. 2.

Reorganization of the connectivity of patients is captured by complexity

Fig. 2, Fig. 3 show that strong connections are heterogeneously concentrated on several specific brain areas, instead of uniformly spreading across the entire brain (shown in Fig. 3). We use 90 AAL brain areas as references for our discussion below. Note that the pairwise correlation has been computed for all 3599 subRs to the i-th subR. The heterogeneity of the strong connectivity of i-th subR is conveniently described by the complexity measure of Shannon entropywhere p is the frequency of connections that belong to brain area j. We used a histogram to estimate the frequency p and hence 90 AAL brain areas provide a natural choice of bins for the estimator. The advantage of this choice of bins is that it does not cross the important anatomical boundaries of the original AAL atlas, which helps to interpret the connectivity reorganization patterns. We note that the entropy is a traditional tool for characterizing the connectivity heterogeneity of neural networks (Tononi et al., 1994). ℋ is small if connections are concentrated on a small number of brain areas forming clusters of connectivity in a resulting network. In contrast, ℋ is large if connections are more homogeneously distributed. Simply speaking, the entropy of each partition of the correlation space is computed as ℋ(i, A), ℋ(i, B), ℋ(i, C), ℋ(i, D), and ℋ(i, E).

For strong connections of the i-th subR of a patient, we propose a complexity-based distance function that characterizes the deviation of the clustering of strong correlations of a patient from the averaged group of healthy controls. More specifically, taking the set B of the strong connections of the i-th subR as an example, we compute the square distance of the complexity measure of a patient in comparison to the averaged healthy group aswhere ℋ(i, B) is the entropy value of patient subject P, and respectively, ℋ〈(i, B) is for the healthy group. Analogously, one can compute ε2(i, D) and ε2(i, E) to capture the reorganization patterns. We find that the patient group has larger entropy values than the healthy controls (see SM-VI, Fig. S5). In other words, the strongest connections of a patient are spread across the entire brain more uniformly than in healthy subjects.

Complexity distance functions are correlated to clinical variables

As discussed above, the reorganization of strong connections is largely captured by the two sets D and E, which represent broken strong links of the healthy controls and newly built connections of a patient, respectively. Therefore, we propose the following characteristic distance function as an indicator of the connectivity reorganization patterns of a patient.

The correlation of the characteristic complexity distance function (Eq. (4)) with the clinical variable FMA is shown in Fig. 5. Connectivity re-organization is much more significantly expressed if a patient is severely affected by stroke (Fig. 5B). For instance, for patients with severe motor impairments (i.e., FMA ≤10), the corresponding linear correlation value −0.36 (with p1-value <0.05) suggests that such patients experience a greater amount of rewiring of the strong connections. Here the term reorganization includes two aspects: (i) strong links that would be expected for healthy controls, and (ii) strong links that have been newly created in patients in comparison to the averaged healthy controls. The correlation between Eq. (4) and FMA has been further confirmed by voxel-wise based regional homogeneity statistical analysis (see SM-VII).

Fig. 5

Correlation of the characteristic complexity distance function ε2 (Eq. (4)) to the clinical variable of FMA. (A) bar charts of FMA, ε2 for each patient subject, P ∈ [38, 69]. (B) scatter plot between FMA and ε2. The dashed vertical line corresponds to the empirical value of FMA =10. The linear fit (ρ1 =  − 0.36) for patients of FMA ≤10 suggests that ε2 is correlated to FMA (p1-value <0.05). By contrast, no significant correlations to FMA have been found for patients of FMA >10 (p2-value >0.1).

We further perform a two-sample t-test for lesion volume (ml) between the patients with FMA >10 and those with FMA ≤10, not showing significant difference between these two subgroups (Mean ± Std: 7.85 ± 5.86 ml for FMA ≤10); 9.32 ± 5.38 ml for FMA >10; p = 0.44). We also performed a voxel-based lesion-symptom mapping (VLSM) analysis (Bates et al., 2003; Liu et al., 2018). The VLSM results suggested that stroke patients do not show significant correlations (at a voxel-wise threshold of p < 0.001, uncorrected) between lesion location at chronic stage and both FMA scores and entropy values. The details of the VLSM analysis are presented in SM-VIII. Additionally, the lesion overlap of the two subgroups was respectively shown for these two sub-groups (see SM-I, Fig. S2). From the visual inspection, the lesion locations were qualitatively similar in both sub-groups. Therefore, the complexity distance function captures the functional connectivity reorganization patterns, providing insights for discriminating patients that are severely affected by stroke.

Discussion

Practically, in brain functional network analysis, there exists no unique criterion for choosing what should represent a network node. This fact has resulted in a wide range of analyzed network sizes, ranging from 90 to 14,000 nodes (Bullmore and Sporns, 2009; Hayasaka and Laurienti, 2010; Matthew et al., 2013; Eguíluz et al., 2005). Our up-sampling procedure based on the AAL atlas is a compromise between noise and area averaging effects when obtaining time series. This sampling technique does not cross the important anatomical boundaries specified by the original AAL atlas, which is widely used (Hermundstad et al., 2013). More importantly, unlike most existing studies, we propose to extract anomaly time series from BOLD signals. The resulting anomaly networks focus on characterizing the deviation of each individual subject from the averaged behavior of the healthy control group, reflecting the abnormal features that are crucial to stroke.

Small-worldliness and a pronounced modular architecture have been observed as general properties of the human brain (Bullmore and Sporns, 2009; Goñi et al., 2014), which are correlated with better cognitive performance (van den Martijn et al., 2010). In the particular case of stroke lesions, altered topological properties of brain networks have been observed across different spatial scales, such as, reduced inter-hemispheric functional connectivity between cortical motor areas, which correlates with the severity of motor deficits (Anne and Grefkes, 2013; Fallani et al., 2013); and reduced network efficiency even in patients with good clinical recovery (Anne and Grefkes, 2013; Bullmore and Sporns, 2012). Furthermore, hubs of brain networks have been shown to be central in clinical disorders (Zac et al., 2015; Baronchelli et al., 2013). All these results are in favor of the random reorganization hypothesis. In contrast, we find that global properties (including small-worldliness, modularity etc.) are conserved in the patient group, challenging the random re-organization hypothesis.

Our data-driven algorithm identifies 15 important subRs as ROIs. With reference to the AAL template, these subRs mainly belong to four crucial brain areas, left precentral gyrus (PreCG.L), right precentral gyrus (PreCG.R), right postcentral gyrus (PoCG.R) and left postcentral gyrus (PoCG.L). All these brain regions are related to sensorimotor functions and frontoparietal control systems. The sensorimotor system is prevalently activated during motor execution and somatosensory information processing (Lacourse et al., 2005; Yang et al., 2007). In contrast, the frontoparietal system plays a crucial role in high-order motor-cognitive processing, such as motor imagery and control of goal-directed movement (Yang et al., 2007; Hanakawa et al., 2008; Haaland et al., 2000; Filimon, 2010). Both systems are engaged in cerebral reorganization following stroke with abnormal task-evoked brain activation (Calautti and Baron, 2003; Ward et al., 2003), as well as disrupted resting-state brain connectivity (Park et al., 2011). However, task-based fMRI techniques cannot easily distinguish whether abnormality of brain activation is induced by a task or is intrinsically disturbed even in resting state. In addition, the variations of global network measures obtained by graph analysis from resting-state data cannot identify to which brain region the abnormality should be attributed, since the correlation coefficient is a symmetric measure. Therefore, our proposed method rectifies the short comings of current methods, and also provides new insights into functional reorganization after stroke.

When a brain is affected by stroke, the connectivity reorganization includes two aspects: a break-down and rewiring of the existing links, and creation of new links. We show that the strongest links in the healthy controls are conserved in patients, while links of intermediate strength are destroyed in a non-random manner. In addition, the newly built strong connections in patients result from random establishment of links across all possible brain regions. These connectivity reorganization mechanisms are quantified by the characteristic complexity distance function ε2. Further work should focus on characterizing spatial anomaly patterns in both functional and structural connectivity (e.g. from diffusion tension image data), and in task-based fMRI data (Hermundstad et al., 2013; Honey et al., 2009). In combination with perfusion techniques, we need to uncover specific neural dynamic changes after stroke by investigating hemodynamic response function models.

Various rehabilitation approaches have been suggested for patients with stroke. Longitudinal studies may provide insights for brain networks at different stages during the recovery processes (Liang et al., 2010). It is obvious that one has to construct dynamic temporal networks, as environmental contingencies vary in the re-organization processes associated with treatments of stroke (Vértes and Bullmore, 2015). It would be expected that FN analysis can shed light on the distinctive roles of various therapies for rehabilitation. It may also be expected that complexity in connectivity re-organization will prove to be a dynamically characteristic variable over longer time scales and useful for lesion inference (Zavaglia et al., 2015). Additionally, it is necessary to construct a database of healthy subjects as a reference and then compare the patients to this reference. In this respect, computational models may provide insightful predictions for the specific profile of neural or behavioral changes following stroke (Forkert et al., 2015).

It is known that distinct functional reorganizations are involved in different activation patterns following cortical and subcortical strokes, presenting highly heterogeneous lesion locations and sizes. In this work, we focused on brain reorganization in subcortical patients only. Therefore, further studies with our approach in patients with cortical lesions are required in order to check whether the functional connectivity reconfiguration mechanisms as presented in this work hold for more general interpretations in other patient group. In addition, one would need a larger sample size of stroke patients to increase the statistical power of the correlation between the complexity measures and the clinical FMA scores. In particular we plan to check whether similar signatures of connectivity reorganizations hold for moderate stroke patients of mild FMA scores. For moderately impaired patients, the statistical evaluation may even challenge the connectivity reorganization mechanism since moderately impaired patients could benefit more readily from re-adaption after stroke. Although the present study focused on patients in the chronic state, the method shows potential for exploring the prediction of training and rehabilitation in acute state of stroke patients.

In summary, we have proposed a unified framework for assessing the hypothesis of random reorganization of FNs following stroke. Based on individual anomaly brain functional networks, we identified regions of interest which are most stroke relevant using a data-driven algorithm. We suggest that entropy complexity measures can be used to assess the random reorganization mechanisms, which are significantly correlated with the clinical assessments.