## Abstract

Children who sustain a prenatal or perinatal brain injury in the form of a stroke develop remarkably normal cognitive functions in certain areas, with a particular strength in language skills. A dominant explanation for this is that brain regions from the contralesional hemisphere “take over” their functions, whereas the damaged areas and other ipsilesional regions play much less of a role. However, it is difficult to tease apart whether changes in neural activity after early brain injury are due to damage caused by the lesion or by processes related to postinjury reorganization. We sought to differentiate between these two causes by investigating the functional connectivity (FC) of brain areas during the resting state in human children with early brain injury using a computational model. We simulated a large-scale network consisting of realistic models of local brain areas coupled through anatomical connectivity information of healthy and injured participants. We then compared the resulting simulated FC values of healthy and injured participants with the empirical ones. We found that the empirical connectivity values, especially of the damaged areas, correlated better with simulated values of a healthy brain than those of an injured brain. This result indicates that the structural damage caused by an early brain injury is unlikely to have an adverse and sustained impact on the functional connections, albeit during the resting state, of damaged areas. Therefore, these areas could continue to play a role in the development of near-normal function in certain domains such as language in these children.

- functional connectivity
- network model
- prenatal/perinatal stroke
- resting-state fMRI
- structural connectivity

## Introduction

Early brain injury in the form of a stroke can alter the structural and functional connectivity (FC) among brain regions across and within hemispheres. These alterations can lead to deficits in many domains such as language, spatial reasoning, and motor function. In the language domain, unlike adults, who show significant impairment upon sustaining damage to left hemisphere brain areas, children with prenatal or perinatal left focal brain injury show remarkably good performance, developing skills in the low-normal to normal level ranges (Woods and Teuber, 1978; Bates and Roe, 2001; Stiles et al., 2005). A predominant theory underlying this near normal language development subsequent to the injury is via compensation or takeover of function by regions from the nondominant hemisphere (Staudt et al., 2002). However, some studies have shown that early left focal injury does not invariably result in recruitment of homologous areas in the right hemisphere (Rasmussen and Milner, 1977; Hertz-Pannier et al., 2002) and, in fact, patients with damage to classical language areas are more likely to show language-related activity in the left hemisphere (Liégeois et al., 2004), suggesting that perilesional regions from the injured hemisphere contribute to the development of normal function. Indeed, it has been shown that preservation of function in these perilesional left hemisphere regions supports improved behavioral outcomes after early stroke (Fair et al., 2010, Raja Beharelle et al., 2010). Overall, the outcome patterns of neural organization after early focal brain injury are highly variable and it is unclear how differential activity in perilesional damaged and contralateral regions plays a role in the preservation or impairment of function.

More fundamentally, it is difficult to disentangle changes in the functional architecture due to neonatal stroke itself compared with changes potentiated by postinjury reorganization. A computational model can be useful to address this issue because one can simulate the impact of a lesion caused by stroke on the structural anatomical connectivity and assess the resulting simulated FC compared with the empirical one.

In the present study, we investigated resting-state FC in 10 children who had suffered a prenatal or perinatal stroke and compared each with computationally modeled effects of the individual anatomic lesions on FC and computationally modeled healthy FC. Our modeling approach involves a mathematical model for each of the brain areas that are coupled through structural connectivity (SC) information for each participant. We obtained the SC for each participant by simulating a lesion in an average SC of a healthy adult brain. We inserted this SC data into a mean field model (Deco et al., 2014a), also available in The Virtual Brain (Sanz Leon et al., 2013), and obtained a simulated FC matrix for each individual for comparison with empirical FC. Based on previous research demonstrating the importance of maintaining a more typical pattern of neural activity in improved outcomes (Raja Beharelle et al., 2010), we hypothesized that the empirical FC after early focal injury would more closely resemble healthy FC patterns than connectivity changes that would be predicted after a lesion.

## Materials and Methods

#### Patient- and age-matched control participant characteristics

Thirty-one typically developing children (age 7–10 years, mean ± SD: 9 years 0 months ± 8 months) participated in the study as age-matched control (AMC) participants. All but four of the control participants were right-handed as determined by the Edinburgh Handedness Inventory (mean ± SD handedness of left handers: −17 ± 0.8, mean ± SD; mean handedness of right handers: 16 ± 0.7; Oldfield, 1971). Formal handedness information was not available for one control participant, who did, however, report right-handedness. Upon inspection, the brain activity pattern for each left-handed control participant did not differ qualitatively from that of the right-handed control participants and they were therefore included in all analyses. Eleven children with prenatal or perinatal unilateral focal brain injury due to stroke (age 7–9, mean ± SD = 8 years 8 months ± 7 months) participated in the study. All of the children with early left hemisphere injury were left handed and the two children with right focal brain injury were right handed as determined by the Edinburgh Handedness Inventory (mean ± SD handedness of left handers: −19 ± 0.8, mean ± SD handedness of right handers: 16 ± 1.5; Oldfield, 1971). All participants were recruited regardless of race, sex, hemiparesis, or lesion site and were screened to exclude individuals with additional preexisting neurological and psychiatric disorders. The parents or guardians of all participants gave formal consent; the participants gave their assent in accordance with the Institutional Review Board of the Biological Sciences Division of The University of Chicago, which approved the study. The study was performed at The University of Chicago when three of the authors (A.R.B., A.S., and S.L.S.) were there. Only children with prenatal or perinatal stroke, which is defined by the National Institute of Child Health and Human Development and the National Institute of Neurological Disorders and Stroke of the National Institutes of Health of the United States as occurring between the 20^{th} week of gestation and the 28^{th} postnatal day (Raju et al., 2007), were included in the study. Early brain injury participants were recruited through local and national support groups for childhood “hemiplegia” (hemiparesis). as well as by referrals from area neurologists and psychiatrists. One child was excluded from the analyses because the lesion was too small to be identified radiologically. Lesion etiology was determined as either having occurred due to periventricular hemorrhagic infarction, mainly occurring in the terminal vein (“periventricular”; *n* = 2), or ischemic infarction of the middle cerebral artery (“vascular”; *n* = 9; Kirton et al., 2008). Participant demographic and lesion etiology information is given in Table 1.

#### Empirical FC

##### fMRI data acquisition.

Resting-state data were collected for 5 minutes. Participants were instructed to keep their eyes closed for the duration of the run. fMRI was recorded on a Siemens 3.0 T TIM Trio MRI scanner with a 32-channel head coil. The fMRI BOLD signal was measured with a susceptibility-weighted (T2* gradient echo) echoplanar imaging (EPI) sequence. The following parameters were used: repetition time (TR) = 2000 ms, echo time (TE) = 20 ms, flip angle = 75°, matrix size = 128 × 120, slice thickness = 3.4 mm with a 1 mm gap, number of axial slices = 32. After functional image acquisition, a high-resolution MPRAGE T_{1}-weighted 3D structural image was acquired in an oblique sagittal plane for each participant (TR = 2300 ms, TE = 2.99 ms, flip angle = 9°, matrix size = 256 × 256, slice thickness = 1 mm, number of slices = 256).

##### fMRI preprocessing.

Anatomical scans were transformed from an oblique dataset to a cardinal dataset at the angle of the functional scans with AFNI. The time series of functional scans were then despiked and slice-timing correction was performed. The anatomical T_{1}-weighted scan was aligned to a base functional EPI run that was the sixth image of the run. The first five images were not used as the base functional run because the scanner signal was not at full strength for those images. The remaining functional scans were then coregistered to the base run using six-parameter rigid-body transformation to correct for head motion throughout the scan.

Because studies examining resting-state data have shown that head motion can have an effect on FC measures (Van Dijk et al., 2012, Yan et al., 2013, Power et al., 2012, 2014, 2013), we assessed whether there were differences between patients and healthy controls on the six motion parameters for rigid motion correction as well as frame-by-frame displacement (Power et al., 2012) using a two-sample *t* test on the mean of the head displacement parameters. The framewise displacement parameter was obtained using the “dvars” metric in FSL.

The nuisance regression was performed by including the motion parameters and the two nuisance regressors extracted from CSF and white matter as the only terms into AFNI's 3dDeconvolve function and then outputting the error term to use as “cleaned” data for the FC analysis. We did not include the nuisance regressors as covariates of no interest in the FC analysis.

We performed additional preprocessing steps to “clean” the functional time series of sources of spurious variance. Small regions of interest (ROIs; 5 mm^{3}) were drawn in the white matter and ventricles for each participant and the average signal from each of these regions was regressed from the functional time series (Fox et al., 2005). The rationale for choosing small ROIs was to ensure that the ROIs would not contain any partial volumes from gray matter. Data were then spatially smoothed with a 6 mm full width at half maximum Gaussian kernel.

##### Lesion mapping.

Lesion masks were manually traced in native space according to visible damage on a T1-weighted anatomical scan. The peri-infarct region was considered outside of the mask. If the lesion impinged on any ROIs in the participants, then the fraction of overlap between the volume of the lesion and the total volume of the ROI was calculated. To calculate this fraction of overlap between lesion and standardized ROI, lesion masks were transformed spatiallyinto Talairach space (Talairach and Tournoux, 1988) using a procedure called Virtual Brain Transplantation (VBT; Solodkin et al., 2010). VBT consists of filling in the lesioned area with “transplanted virtual tissue” from the nonstroke hemisphere, which provides anatomical landmarks for the standard alignment algorithms that would be missing due to the tissue damage caused by the lesion. In addition, VBT takes into account that the ipsilesional hemisphere is likely atrophic and distorted (Bucy, 1964; Warabi et al., 1990; Graham and Landtos, 2002; Kabir et al., 2007) and therefore scales proportionally the hemisphere containing the “transplanted virtual tissue” (taken from the nonstroke hemisphere) to match the size of the stroke hemisphere before filling it in. VBT was used in participants in whom the cortical stroke was large enough to affect the gray matter (Patients 06–10) to minimize distortion during spatial normalization. The total volume of the ROI was taken as the total volume of the ROI in the standardized parcellation template used for the FC analyses (described in further detail below; Desikan et al., 2006). Lesion masks drawn in native space for each participant are shown in Figure 1.

##### Generation of FC matrices.

The gray matter was parcellated into 34 regions per hemisphere (68 areas total) according to anatomical landmarks (Desikan et al., 2006). These regions were then more finely parcellated into 998 ROIs of approximately equal volume. The average time series of each ROI was calculated and FC was examined by calculating pairwise correlations between the BOLD signal time series in the ROIs, which generated a correlation matrix for each child and AMC participant. The AMC FC values between 68 brain areas were obtained by averaging across 27 AMC participants; those between 998 ROIs were obtained by averaging across 17 AMC participants. The reduced number of healthy participants at the fine-grained parcellation was due to the lack of complete overlap between the functional data of some participants with the structural template.

#### Language tests

Outside of the fMRI scanning session, child participants were given the Letter–word Identification subtest of the Woodcock-Johnson III testing battery (Woodcock et al., 2001), which measures reading aloud of letters, single words, and nonsense words. Participants were also given the Word Attack subtest of the same testing battery, which requires participants to read pronounceable nonwords (pseudowords). The Word Attack subtest will hereafter be referred to as the nonword reading test, which is more explanatory of its content.

#### SC

Neuroanatomical SC for the healthy (control) comparison was obtained using tractography estimates from diffusion spectrum imaging (DSI) data in 10 healthy right-handed adult male human participants (mean age 22.5 ± 1.8 years). DSI and MPRAGE T1-weighted acquisitions were performed on a Siemens 3.0 T TIM Trio MRI scanner equipped with a 32 channel head coil. The DSI sequence included 128 diffusion-weighted images with a maximum *b*-value of 8000 s/mm^{2} and one *b*_{0} reference image. The acquisition volume contained 96 × 96 × 34 voxels with 2.2 × 2.2 × 3 mm resolution. TR and TE were 6100 and 144 ms, respectively. The MPRAGE acquisition had a 1 mm in-plane resolution and 1.2 mm slice thickness, covering 240 × 257 × 160 voxels. The TR, TE, and TI were 2300, 2.98, and 900 ms, respectively.

Subject-wise cortical ROIs were generated by registering two distinct parcellation templates to the T1 participant space. The following templates were used: the Desikan–Killiany anatomical atlas, composed of 68 brain areas (Desikan et al., 2006), and the 1000-ROIs template defined by Cammoun et al. (2012). The surface-based registrations were performed using FreeSurfer software (http://surfer.nmr.mgh.harvard.edu). The 1000-ROIs template was generated (Cammoun et al., 2012) by means of a recursive region-growing process within the 68 brain areas of the Desikan–Killiany atlas. Each one of the 1000 ROIs can therefore by associated to a specific anatomical parcel of the Desikan–Killiany atlas. We evaluated the distribution of sizes for the 1000 ROIs of the 10 healthy participants included in this study and used them for the SC estimation. The average ROI volume was 516.66 mm^{3} (SD 187.41 mm^{3}; Fig. 2).

Deterministic streamline tractography within the white matter was used to estimate the fiber tract density connecting each pair of ROIs. Anatomical connectivity strength between every pair of ROIs was calculated by dividing the number of connecting fibers by the average area of the two ROIs and by average fiber length between the two ROIs (Hagmann et al., 2008; Honey et al., 2009). The 10 weighted connectivity matrices were then averaged across participants.

Because the dynamical model of one region already takes into account the effect of its internal connectivity (see below), the connection of a region to itself was set to 0 in the connectivity matrix for the simulations.

#### Lesion modeling

We modeled the effects of lesion using three different approaches. The information about the extent of the brain lesion for each affected child was computed as a fraction of the voxels belonging to each ROI that were damaged. In the first two approaches, we used this information to change the structural connections of the damaged ROIs. First, because the FC data were parcellated into 998 ROIs, we chose the SC between corresponding ROIs. In the first approach, for each ROI that had at least 10% damage, we decreased the strength of all its existing connections in proportion to the damage, resulting in a 998 dimensional “lesioned” SC. In the second approach, for each ROI that had at least 10% damage, we set the weakest number in proportion to the damage of all its existing connections to zero. In other words, we pruned the connections of the damaged ROIs in this approach and obtained the corresponding 998 dimensional lesioned SCs. In each of the two approaches, we used this fine parcellation SC matrix to calculate the corresponding lesioned SC at the coarse parcellation, consisting of 68 brain areas, by aggregating the connectivity strengths between all ROIs belonging to every pair of brain areas in question.

In the third approach, we did not change the structural connections of any of the ROIs, but decreased the mean activity of excitatory population in the local model (described below) of each damaged brain area in proportion to the damage.

#### Dynamic mean field model

To obtain simulated FC among different brain areas, we first computed the global dynamics of each area using a reduced dynamic mean field (DMF) model (Wong and Wang, 2006). The DMF model is an approximation of a spiking network model (Deco et al., 2012, 2013) consisting of populations of excitatory and inhibitory neurons and expresses the averaged activity of each population using a single variable, namely the average firing rate. In the DMF approach, each population firing rate depends on the input currents into that population, which in turn depend on the firing rates. Therefore, population firing rates can be determined self-consistently by a reduced system of coupled nonlinear differential equations.

In our case, we represent each brain area by a network of one excitatory population and one inhibitory population. Inter-area connections are modeled as long-range synaptic connections between the excitatory pools in those areas. They are weighted by the strengths specified in the neuroanatomical matrix SC, as obtained by DSI described above, and a global scaling factor *G*. We use this scaling factor as a control parameter to identify the regime in which the simulated FC correlates best with the empirical one. The equations describing the population activities of each brain area are as follows:
Where *r _{i}*

^{(E,I)}denotes the population firing rate of the excitatory (

*E*) or inhibitory (

*I*) population in the brain area

*I*;

*S*

_{i}^{(E,I)}denotes the average excitatory or inhibitory synaptic gating variable at the local area

*I*; the input currents to the excitatory or inhibitory population

*i*are given by

*I*

_{i}^{(E,I)}

*I*

_{external}encodes external stimulation for simulating task-evoked activity;

*C*

_{ij}denotes the strength of the SC, that is, the fiber density as measured by DSI, between areas

*i*and

*j*;

*w*is the strength of local excitatory recurrent connections; and

*u*

_{i}is uncorrelated standard Gaussian noise. Values of all constants used in these equations are listed in Table 2.

*J*_{NMDA} denotes the local and long-range excitatory synaptic coupling constant and *J*_{i} denotes the local feedback inhibitory synaptic coupling. As in Deco et al., 2014b, we adjust this local feedback inhibition for every value of the global coupling parameter *G* so that every brain area network is in a state of balanced excitation and inhibition. In a balanced spontaneous state, the average firing rate of the excitatory population is confined ∼3 Hz, which constrains the input current to the excitatory population to satisfy the condition: *I _{i}^{E}* −

*b*/

_{E}*a*= − 0.026. We carry out the feedback inhibition control for each value of the control parameter

_{E}*G*by adjusting the value of

*J*

_{i}so that this condition is satisfied. Specifically, we simulate the system of Equations 1–6 for 10 s and compute the time average of

*I*

_{i}^{(E)}. If

*I*−

_{i}^{E}*b*> − 0.026, then we upregulate the corresponding local feedback inhibition

_{E}/a_{E}*J*=

_{i}*J*+ Δ; otherwise, we downregulate

_{i}*J*=

_{i}*J*− Δ. We repeat this procedure recursively until the constraint on the input to the local excitatory pool is fulfilled in every brain area.

_{i}Using this feedback inhibition control, we simulate the system of stochastic differential Equations 1–6 for 10 min to obtain the simulated neuronal activity of each brain area. We then transform it to BOLD signal using the Balloon–Windkessel hemodynamic model (Friston et al., 2003). Finally, we calculate the simulated FC matrix in which each entry is a correlation between simulated BOLD signals of the corresponding pair of brain areas.

#### Calculation of simulated FC using statistical moments method

We use the DMF model described in the previous section to calculate the simulated FC matrix at a coarse parcellation of 68 brain areas. However, it is computationally very expensive to calculate the FC at the fine parcellation of 998 brain areas using the DMF. Therefore, to calculate the FC at such a high resolution, we follow Deco et al. (2013, 2014b), wherein deterministic dynamical equations for statistical moments of the network are derived from the coupled equations of the DMF. Specifically, we obtain the equations for the expected mean gating variable μ_{i}^{(m)} of each local neural population (where *m* = *E* or *I*) of the cortical area *i* and the covariance *P*_{ij}^{(mn)} between gating variables of neural populations of type *m* and *n* of local cortical areas *i* and *j*, respectively. The moments are defined as follows:
Where the angular brackets <.> denote the average over realizations. The equations of motion for the mean values of gating variables take the following form (Deco et al., 2014b):
Where *u*_{i}^{(m)} is the mean input current to the neural population *m* = *E*,*I* of cortical area *i*, defined as follows:
The equations for the second order moment, the covariance, can be written in matrix form. If *P* is the covariance matrix between gating variables, defined as follows:

**P** = 〈*S⃗S⃗ ^{T}*〉 =

*T*is the transpose, then the equation of motion for the covariance matrix can be written as follows (Deco et al., 2014b): Where

*Q*

_{n}is the covariance matrix of the noise (which is diagonal for uncorrelated white noise) and

*A*is the Jacobian matrix as follows:

These derivatives in the *A* matrix can be written as follows:
Where **K** = *K*_{ij}^{(EI)} = − *J _{i}*δ

*,*

_{ij}*K*

^{(IE)}=

*J*

_{NMDA}**I**,

*K*

^{(II)}= −

**I**, and

*K*

^{(EE)}=

*w*

_{+}

*J*

_{NMDA}**I**+

*GJ*

_{NMDA}**C**, where

*I*and

*C*are the

*N*-by-

*N*identity matrix and anatomical connectivity matrix, respectively.

Equations 13–14 indicate that the covariance of gating variables is determined by both the underlying connectivity and the dynamics. In the stationary regime, the covariance matrix of fluctuations around the spontaneous state is given by the following algebraic equation:
Which can be solved using the eigen decomposition of the Jacobian matrix evaluated at the fixed points of Equations 9–10. In particular, *A* = *LDL*^{−1}, where *D* is a diagonal matrix containing the eigenvalues of *A*, noted λ_{i}, and the columns of matrix *L* are the eigenvectors of *A*. Multiplying Equation 15 by *L*^{−1} from the left and by *L*^{−}† from the right († being the conjugate transpose), we get the following:
Where *M* is given by the following: *M _{ij}* = −Q̃

_{ïj}/(λ

*+ λ*

_{i}_{j}*) and

**Q̃ = L**.

^{−1}Q_{n}L^{−†}Therefore, the covariance matrix of spontaneous fluctuations is determined by the eigenvalues of the Jacobian matrix, which in turn is related to the connectivity matrix and the dynamics. Therefore, Equation 15 provides a direct link among the correlation structure, the underlying connectivity, and dynamics. Indeed, the interpretation of Equation 15 is that the input covariance (*M*) is propagated through the dynamical system and is mapped to its approximated output (*P*).

Finally, to compare the results from the above moments' method with the empirical FC, we estimate the correlation matrix, noted *Q*, between the gating variables of the excitatory populations, defined as follows: *Q _{ij}* =

*P*

_{ij}

^{(EE)}/

## Results

In the first step, we compared the performance of the participants on language tasks with those of age-matched typically developing children. Patients scored an average of 44 ± 8 (of 76; 102 ± 11 for normative score) on the letter–word identification test and 15 ± 7 (of 32; 100 ± 8 for normative score) on the nonword reading test. Data for both of these tests were not available for one patient. Patients had significantly lower performance on letter–word identification (*t*_{(37)} = −2.22, *p* = 0.032) and nonword reading (*t*_{(37)} = −2.32, *p* = 0.026) compared with controls (Raja Beharelle et al., 2012). Next, we compared the FC values for every patient with those averaged across AMC participants. The comparisons were made by grouping together the FC values of all and only the damaged among the 68 (998) brain ROIs with all other ROIs. Figure 3*a* displays these two comparisons for the case of coarse-grained parcellation of 68 areas; Figure 3*b* displays them for the case of fine-grained parcellation of 998 ROIs. The correlation values for both comparisons were >0.4 in the case of coarse-grained parcellation and were much smaller in the case of fine-grained parcellation.

We found no significant differences between healthy, typically developing children and patients in the means of any of the six motion parameters taken across the functional scans (mean motion1: *t*_{(36)} = 0.3564, *p* = 0.3564; mean motion2: *t*_{(36)} = 1.5785 *p* = 0.1232; mean motion3: *t*_{(36)} = −0.1082 *p* = 0.9145; mean motion4: *t*_{(36)} = −1.8920, *p* = 0.0666; mean motion5 = *t*_{(36)} = 0.1317, *p* = 0.8960; mean motion6 = *t*_{(37)} = −1.7088, *p* = 0.0961). In addition, we found no significant differences in the mean framewise displacement (*t*_{(36)} = 0.7999, *p* = 0.4290). Together, these results suggest that any group effects in the FC results are not attributable to group differences in head motion.

Next, we sought to obtain model FC matrices using a healthy adult SC matrix averaged across 10 participants and a lesioned SC matrix for each child participant. Here, we used the first approach as described in the Materials and Methods to simulate the impact of lesion by decreasing the structural connections of each damaged area to all other areas in proportion to the damage percentage suffered by each such area. At the coarse parcellation of 68 brain areas, we inserted the healthy SC into the stochastic DMF model and simulated it to obtain model FCs at several values of the global coupling parameter *G*. We then compared, for every *G*, these simulated FC values of all brain area pairs with the corresponding empirical FC values averaged across 27 AMC participants using Pearson's correlation measure. The value of *G* at which the correlation coefficient was a maximum was identified as the optimal working point of the model (Deco et al., 2014a, 2014b) and the corresponding simulated FC matrix was used as the healthy model FC matrix. We repeated this process for each of the brain-injured children by using the corresponding lesioned SC matrix in the model and comparing the model FCs at several values of *G* with the empirical FC of that child. The simulated FC matrix corresponding to the value of *G* at which the correlation was maximum was identified as the model FC matrix for that child. Figure 4*a* plots the correlation fit between simulated and empirical values of FC for the healthy case (blue) and for the case of Participant 01 (red). The optimal working point for each case was different. In addition, the maximum correlation between model FCs and the empirical FC averaged across AMC participants was higher than for Participant 01. Figure 4, *b* and *c*, shows the corresponding scatter plots between values of the optimal model FC and the empirical FC averaged across AMC participants and of Participant 10, respectively.

Next, we tested whether the resting-state FC for each of the brain-injured children resembled more the resting-state FC of a healthy participant or that of a participant with a simulated lesion from stroke. Such comparisons can be made empirically in patients who suffered stroke in adulthood because functional neuroimaging data can be obtained at several points of time after the onset of stroke in them. However, in the case of children who suffered stroke during the neonatal period, only the first comparison can be made (Fig. 3) empirically. This is precisely the advantage of a measure based on a computational model in which the effects of a lesion can be simulated.

Therefore, we compared, for every child, the maximum correlation between empirical FC and the optimal model FC obtained using the corresponding lesioned SC with the correlation between the empirical FC and the optimal model FC obtained using an average healthy SC matrix. Figure 5*a* displays the comparison between the values of all brain area pairs for each of the 10 children. Figure 5*b* compares the two correlations averaged across all affected children. Figure 5*c* compares the two correlations averaged across only the most affected children. The most affected children (Participants 01, 06, 07, 08, 09, and 10) had a total damage percentage >25%, as can also be seen in Table 1. As the figure shows, the difference is marginal when the connectivity values of all brain area pairs are together.

Next, we made the comparisons using the FC values of only the damaged areas with all other areas. As before, Figure 6*a* displays the comparison for each child and Figure 6, *b* and *c*, displays the comparison of correlations averaged across all and the most affected children, respectively. Here, we found that the empirical FC values of damaged areas with all other areas in most affected children displayed greater correlation, on average, with their optimal simulated FC values in a healthy adult brain than those in a lesioned brain and the difference was statistically significant.

We also used the FC values of only the contralesional areas with all other areas for further comparisons. As shown in Figure 7, there was no statistically significant difference between the correlation of the empirical FC with the simulated FC for an average healthy adult and with that for an average child with stroke.

We repeated these investigations at the fine-grained parcellation of 998 ROIs. As mentioned in the description of the model in the Materials and Methods, we used a further approximation of the DMF model in this case that allows calculation of the moments of the gating variables of each population in each region. Figure 8*a* displays the comparison, for every child, between intrahemispheric and interhemispheric connectivity values of only the damaged ROIs. Figure 8, *b* and *c*, displays the comparison of correlations averaged over all and most affected children, respectively. Again, the empirical FC values of only the damaged ROIs, especially in the most severely affected children, correlated better, on average, with their optimal simulated FC values in a healthy adult brain than in a lesioned brain. The difference in correlations was statistically significant.

To assess whether the results we obtained so far were related to the specific manner in which we modeled the structural impact of the lesion, we tested two other lesion models described in the Materials and Methods. In the second approach, we constructed the SC matrices for the brain-injured children by pruning the weakest of the connections of damaged areas in proportion to the damage suffered by each area and then simulated the model to obtain the model FC. Figure 9 displays the comparison of correlations averaged across all children (Fig. 9*a–c*, top) and the most affected children (Fig. 9*d–f*, bottom). The correlations between empirical and simulated FC were calculated using all FC values (Fig. 9*a*,*d*, left), FC values of only damaged areas with all other areas (Fig. 9*b*,*e*, middle), and finally using the FC values of only the contralesional areas with all other areas (Fig. 9*c*,*f*, right). We found that there was a significantly higher correlation between empirical FC values of only the damaged areas with all other areas and the corresponding healthy model FC values than that between empirical FC values and the lesioned model FC values.

We tested another possible model for the lesion impact. In this approach, we used the SC matrix, averaged across healthy adults, for each of the brain-injured children in our computational network model, but decreased, in proportion to the damage, the mean simulated activity of the local population of excitatory neurons in each of the damaged brain area. As Figure 10, *b* and *e*, illustrates, the correlation between empirical FC values of only the damaged areas with all other areas and the corresponding lesioned simulated FC values, whether averaged across all participants or only most affected participants, was significantly lower than that between the empirical FC values and healthy simulated FC values.

These results demonstrate that the effects of structural damage of brain regions caused by a prenatal or perinatal stroke on their FC poorly explain the evolved connectivity values observed at a later age, assuming that the lesion does not resolve.

## Discussion

Here, we sought to disentangle the effects of an early brain lesion on the resting-state FC in children with prenatal or perinatal stroke from the changes due to postinjury reorganization. These children displayed a significantly lower performance on the letter–word identification and non-word reading tasks (Raja Beharelle et al., 2012). The resting-state FC values for these children displayed, on average, a good correlation with the corresponding values for typically developing children (Fig. 3). In this study, we used a computational model to simulate a structurally healthy brain and a structurally damaged brain representative of the extent of damage suffered by each child. Correlations between the resulting activities of brain areas yielded simulated FC at a participant-specific optimal working point of the model for every affected child and the average healthy participant.

We found that the empirical FC values of all brain area pairs together did not display, in terms of correlation with their simulated counterparts, a statistically significant difference between healthy and brain-injured cases. However, when we considered the functional connections of only the damaged areas, the empirical values correlated significantly better with the healthy values than with the damaged ones. This result was observed at both the coarse parcellation of 68 brain areas as well as the fine parcellation of 998 ROIs. We emphasize here that this observation is not an artifact of the construction of the lesioned structural matrix because we calibrated the model for each SC, whether healthy or lesioned, using the corresponding FC of the participant. Essentially, for every participant, we find an optimal value of the global coupling parameter *G*, multiplying the corresponding SC at which the simulated FC correlates maximally with the empirical FC (Deco et al., 2014a). This participant-specific optimal working point of the model, at which we obtain the simulated FC for comparisons, is found by correlating the simulated connectivity values of all areas, and not the damaged areas alone, with the corresponding empirical values. We found similar results with three different models of lesion impacts—the first two modeled the impact of lesions on the structural connections of lesioned areas and the third modeled the impact on the local neuronal activity of lesioned areas. The results were stronger in the third approach due to probable loss local excitation–inhibition balance in the model. Therefore, our results are robust against three different possible ways in which a lesion could be simulated in a computational model.

Given that the areas we examined are affected by the lesion, one would expect that their FC patterns should resemble those patterns predicted by a simulated lesion. However, because neurodevelopmental and lesion-related compensatory processes occurred after the time of injury until the time of assessment, our results can be seen as a marker of physiological normalcy of the damaged areas.

Prenatal or perinatal stroke can cause significant alterations in the development of normal anatomical pathways. This structural reorganization can impair the normal function in certain domains while preserving it in others. Unilateral brain damage occurring before or during the time of synaptogenesis of corticospinal motor projections can disturb the normal strengthening of contralateral projections in favor of ipsilateral projections from the affected hemisphere, resulting in the abnormal function of the paretic hand (Staudt et al., 2002, 2004). Conversely, developing thalamocortical somatosensory projections can bypass even large periventricular brain lesions to reach their original cortical targets, thereby preserving the somatosensory function (Staudt et al., 2006). In the language domain, compensatory activity in the nondominant hemisphere for language-related tasks has been proposed as the central mechanism for the preservation of language function (Lidzba and Staudt, 2008). Our findings suggest that, after prenatal or perinatal stroke, resting-state connectivity of perilesional areas most resembles neurologically normal participants. Although there may be strong interindividual differences, our results suggest that development of cortical networks may not be entirely equipotential (Lashley, 1951, Lenneberg, 1967), as studies of extensive cases of early cortical injury such as hemispherectomy would suggest (Muckli et al., 2009). Further, when extreme adaptation of neural networks is not required, there may be some predisposition toward certain types of configurations in the functional architecture that are present in the healthy brain despite an intervening injury. This may be interpreted as being consistent with findings in the language domain showing a very early predisposition toward left frontotemporal processing of speech sounds similar to adults (Dehaene-Lambertz, Dehaene, and Hertz-Pannier, 2002; Dehaene-Lambert et al., 2006) and in the visual domain, where formation of cortical circuits cannot be explained by experience alone and many features of visual system development are instead related to self-generated or spontaneous neuronal activity (Katz and Shatz, 1996). Examining the role of how these endogenously generated patterns of activity help to sculpt the development of functional circuits may provide intriguing clues into the mechanisms underlying the robustness of specific patterns of connectivity to early brain injury.

We also assessed the relation between empirical FC values of all contralesional areas with all other areas and the corresponding values from the simulated healthy and lesioned matrices. The motivation for this comparison lies in the fact that several functional imaging studies, especially those involving language tasks, show higher activations in the regions that are contralateral to the damaged areas (Staudt et al., 2002). Therefore, we wanted to test whether it is also reflected in terms of the resting-state connectivity values of contralesional areas. We found that, although connectivity values of these areas during the resting state displayed, on average, a higher correlation with the simulated lesioned values than with the simulated healthy values, the difference was not statistically significant.

One limitation of our study is the fact that we did not have access to the SC information for a healthy child and/or a child with stroke. Therefore, we had to use the SC averaged across 10 healthy adult participants to obtain an optimal healthy model FC and to construct the lesioned SC for finding the model FC corresponding to each of the affected children. Although there is little difference in the white matter wiring pattern in a healthy human brain after the first 2 years of life (Hagmann et al., 2010), significant increases in the axonal diameter and myelin thickness are known to occur during maturation of these connections (Fields, 2005). These changes have been linked with significant increases in average node strengths and efficiency and a significant decrease in clustering coefficient of the SC with age during first two decades of life (Hagmann et al., 2010). Strengthening of long-range and weakening of short-range connectivity in functional networks of adults compared with children reflects these structural changes (Supekar et al., 2009, Fair et al., 2009). Therefore, this limitation affected the value of maximum correlation between our model FCs and the FC averaged across AMC participants (Fig. 4*a*, blue curve). Second, we did not determine precisely which pathways may be affected in constructing the lesioned SCs by uniformly decreasing the strength of connections to a damaged area in proportion to the percentage of damaged voxels in it. Both of these approximations most probably contributed to the values of correlation between empirical and optimal simulated FC values for each of the affected children, being typically between 0.3 and 0.4.

Another limitation of the study is that we did not model alternate patterns of FC to test whether they provide a better fit to the child empirical FC than those from the healthy, mature architecture. One reason for this is due to the large subset of alternate configurations that could potentially be modeled given that no research has yet presented testable candidates for atypical FC patterns that could arise in compensation for prenatal or perinatal stroke. Therefore, we cannot entirely rule out a stronger correlation occurring with a certain atypical FC pattern. However, our results do demonstrate that functional reorganization after early focal injury, in part during the resting state, involves shifts in connections toward healthy configurations, particularly in the areas affected by the lesion.

Our results demonstrate that computational modeling of the fMRI signals during the resting state can be used effectively to differentiate between effects of a lesion and postinjury reorganization in individuals suffering from stroke. Specifically, in the case of children suffering from prenatal or perinatal stroke, such a model-based approach allows comparison of empirical FC with FC simulated using a structurally damaged brain, which is impossible to do empirically. Finally, and most importantly, our results argue against the hypothesis that FC of brain regions damaged during an early brain injury will evolve significantly differently from their evolution during a typical brain development. In this sense, they indicate physiological normalcy of areas damaged by stroke. This result, although obtained using resting-state activity, lends support to the hypothesis that the damaged areas continue to support the performance on tasks involving them and that maintenance of typical functional connections of these regions plays an important role in the preservation of normalcy after early injury.

## Footnotes

M.H.A. and G.D. were supported by the European Research Council (Advanced Grant DYSTRUCTURE 295129). The imaging data analyzed in this study were obtained as part of a project at The University of Chicago in language development that was supported by the National Institute of Child Health and Human Development–National Institutes of Health (Grant P01 HD040605). Their support is gratefully acknowledged, as is the help of our coinvestigators on that grant, Susan Goldin-Meadow and Susan Cohen Levine, and our research assistant, Victoria Li.

The authors declare no competing financial interests.

- Correspondence should be addressed to Mohit H. Adhikari, Center for Brain and Cognition, University of Pompeu Fabra, Calle Roc Boronat, 138, 08018 Barcelona, Spain. mohit.adhikari{at}upf.edu