An Evolutionary Computation Approach to Examine Functional Brain Plasticity

One common research goal in systems neurosciences is to understand how the functional relationship between a pair of regions of interest (ROIs) evolves over time. Examining neural connectivity in this way is well-suited for the study of developmental processes, learning, and even in recovery or treatment designs in response to injury. For most fMRI based studies, the strength of the functional relationship between two ROIs is defined as the correlation between the average signal representing each region. The drawback to this approach is that much information is lost due to averaging heterogeneous voxels, and therefore, the functional relationship between a ROI-pair that evolve at a spatial scale much finer than the ROIs remain undetected. To address this shortcoming, we introduce a novel evolutionary computation (EC) based voxel-level procedure to examine functional plasticity between an investigator defined ROI-pair by simultaneously using subject-specific BOLD-fMRI data collected from two sessions seperated by finite duration of time. This data-driven procedure detects a sub-region composed of spatially connected voxels from each ROI (a so-called sub-regional-pair) such that the pair shows a significant gain/loss of functional relationship strength across the two time points. The procedure is recursive and iteratively finds all statistically significant sub-regional-pairs within the ROIs. Using this approach, we examine functional plasticity between the default mode network (DMN) and the executive control network (ECN) during recovery from traumatic brain injury (TBI); the study includes 14 TBI and 12 healthy control subjects. We demonstrate that the EC based procedure is able to detect functional plasticity where a traditional averaging based approach fails. The subject-specific plasticity estimates obtained using the EC-procedure are highly consistent across multiple runs. Group-level analyses using these plasticity estimates showed an increase in the strength of functional relationship between DMN and ECN for TBI subjects, which is consistent with prior findings in the TBI-literature. The EC-approach also allowed us to separate sub-regional-pairs contributing to positive and negative plasticity; the detected sub-regional-pairs significantly overlap across runs thus highlighting the reliability of the EC-approach. These sub-regional-pairs may be useful in performing nuanced analyses of brain-behavior relationships during recovery from TBI.

Introduction

The human brain is an ever-changing network of complex interactions between billions of neurons operating at multiple spatio-temporal scales (Hebb, 1947; Diamond et al., 1981; Greenough and Chang, 1988). Recent applications of blood oxygen level dependent (BOLD) functional magnetic resonance imaging (fMRI) have focused on examining the functional relationship between spatially discrete regions at a macroscopic level (Biswal et al., 1995; Raichle et al., 2001; Fox et al., 2005; Greicius et al., 2009). In fMRI literature, the functional relationship between a pair of regions of interest (ROIs) in the brain has been traditionally defined as the degree of temporal synchronization between the regions (Friston et al., 1993). A common research goal in the systems neuroscience is to examine changes in brain connectivity over time. For example, in a study examining the influence of exercise in cognitive aging, investigators examined the connectivity changes between two predefined temporal regions, bilateral parahippocampus, and bilateral middle temporal gyrus (Voss et al., 2013). Alternatively, in the context of learning, changes in the functional relationship between areas associated with motor and sensory cortices, such as A1 and the insula and M1 and thalamus, have been examined following musical training (Luo et al., 2012). The clinical neurosciences have a burgeoning literature examining the consequences of neurological insult and normal and abnormal aging on large-scale brain networks (for review of the literatures see Sheline and Raichle, 2013; Tijms et al., 2013; Hillary et al., 2015). For example, in traumatic brain injury (TBI), there has been increasing interest in documenting the change in pairwise functional relationships among the regions central to large scale networks, such as the default mode network (DMN), salience network, and the executive control network (ECN; Hillary et al., 2011a, 2014; Sharp et al., 2011; Bonnelle et al., 2012; Arenivas et al., 2014; Venkatesan et al., 2014). It is the primary aim in these literatures to understand how the functional relationship between predefined regions in the brain evolves over time.

Common approaches to examine functional brain connectivity include seed-based or anatomical-atlas (van den Heuvel and Hulshoff Pol, 2010) based analysis. In the former, the relationship of a targeted brain-region (a seed) is correlated with all other brain regions in a pairwise manner. Alternatively, examiners may use an atlas-based approach and ROIs are defined using standardized anatomical atlases (e.g., Automated Anatomical Labeling; Tzourio-Mazoyer et al., 2002) and the functional relationship between any pair of ROIs is defined as the correlation between the average signal representing each ROI (Luo et al., 2012; Hillary et al., 2014; Rajtmajer et al., 2015). In an ideal scenario, a ROI would be composed of voxels that represent homogeneous activity. However, the reliance upon structural (as opposed to functional) information (e.g., atlas-based approach) results in ROIs with significant signal heterogeneity. Within ROI heterogeneity may cause the functional relationship between two ROIs to change across time at finer spatial scales than the ROIs themselves, and these changes may not be easily detected using an averaging based approach. To minimize within ROI heterogeneity, a number of data-driven voxel-level methods have been developed to aid in ROI definition (Heller et al., 2006; Yan et al., 2011; Zalesky et al., 2012; Thirion et al., 2014; Yeo et al., 2014; Rajtmajer et al., 2015; Wang et al., 2015). However, traditional ROI parcellation approaches are not designed to detect functional connectivity changes between two ROIs that are occurring at finer spatial scales. Hence there is a need for a data driven approach that can automatically detect the correct spatial scales at which the functional connectivity between 2 ROIs are changing across time.

In order to better examine the fine scale functional connectivity changes between two brain regions across time all voxel-wise interactions between the regional-pairs must be examined. One way to handle this is by defining the strength of the functional relationship between a regional-pair as the number of statistically significant voxel-wise functional connections that exist between the two regions (where two voxels are functionally connected if the activity at each voxel is highly synchronized). However, by simply examining the change in the total voxel-wise connections between two brain regions across time, subtle effects may not be easily detectable. For example, one might imagine a scenario where there exists significant voxel-wise connectivity change over time between two large brain regions; however, due to both connectivity loss and gain in voxel pairings, the net change is close to zero. This would mislead an investigator to believe that no change has occurred. Therefore, to accurately characterize the change in the functional relationship over time between a pair of ROIs, connections between sub-regional-pairs that show significant change in their functional relationship strength must be detected, including both connection gain and loss. We define a sub-region as a cluster of spatially connected voxels within a ROI, where each voxel can have 6 neighbors in 3 dimensions. A sub-regional-pair consists of a sub-region from each of the two ROIs (see Figure 3B).

The goal of this study is to develop a novel approach for examining functional plasticity between two ROIs across time by detecting such significantly plastic sub-regional-pairs. In this article, the change in the strength of the functional relationship between two regions (either the entire ROI-pair or a constituent sub-regional-pair) across time will be referred to as plasticity. In particular, positive plasticity refers to situations where the strength of the functional relationship between two regions becomes stronger (i.e., two regions show an increase in synchronization over time); similarly, negative plasticity refers to cases where the strength of the functional relationship becomes weaker.

We introduce a novel voxel-wise evolutionary computation (EC) based procedure to document plasticity between an investigator defined ROI-pair across two fMRI-scan sessions that are separated by a finite duration of time; the procedure can be run for each subject separately, and it simultaneously utilizes the BOLD-fMRI data from both of the sessions. Of note, we choose a commonly used anatomical atlas for ROI selection (e.g., AAL atlas) and this space is refined over the course of the EC procedure based upon connectivity change. In this sense, the chosen ROI provides an initial constraint of the search space for the algorithm but the procedure is agnostic to the nature of the ROI. The procedure is data-driven, and detects sub-regional-pairs that exhibit significant positive and negative plasticity (see Section Objective Function). These sub-regional pairs are then used for developing a subject specific estimate of overall positive and negative plasticity. The challenge in detecting such plastic sub-regional-pairs is primarily due to the fact that many combinations of the connected subset (see Section Chromosome Encoding) of voxels from each ROI must be probed. Hence, for large ROI-pairs, examining different combinations of sub-regional-pairs becomes computationally very expensive. Therefore, we chose to use an evolutionary computation (EC) based approach. An advantage of using EC is that it is implicitly parallelizable (Rechenberg, 1971; Schwefel, 1974; De Jong, 2006; i.e., EC may be broken into multiple concurrent processes; this can drastically reduce total computation time) and has been widely used in handling problems of similar complexity in other research domains (Gruau and Whitley, 1993; Gruau, 1994; Kobayashi and Ohbayashi, 1999; Roy et al., 2013).

We illustrate (1) the insights gained by considering significantly plastic sub-regional-pairs and (2) the strength of this EC procedure in the study of brain plasticity. To do so, we examine the functional plasticity between four ROI-pairs in traumatic brain injury (TBI) and healthy control (HC) subjects with focus on large-scale networks that are now central to the cognitive neurosciences (e.g., DMN, ECN). Specifically, in a ROI-pair A-B, A is either DMN-frontal left or DMN-frontal right and B is either ECN-frontal left or ECN-frontal right, for a total of four ROI-pairs. The ROIs are defined using AAL atlas (Tzourio-Mazoyer et al., 2002); the brain regions constituting these ROIs are presented in Table S5. These ROIs have been chosen because there have been well documented functional changes in DMN and ECN post TBI (Bonnelle et al., 2011; Sharp et al., 2011; Venkatesan et al., 2014); these results may be used to validate the findings of the EC procedure. To examine the utility of this EC-based approach we conduct the following analyses:

1. Analysis-1: First, we use conventional ROI-averaging based procedure to define functional connectivity between the 4 ROI-pairs for each subject at two different time-points when BOLD-fMRI scans were conducted (see Section Subjects). Then we examine if the functional connectivity change can be detected across the two time-points at a group level for TBI and HC subjects. We hypothesize that this approach will detect little to no plasticity.

2. Analysis-2: Next, we define the strength of functional connectivity between a ROI-pair as the number of statistically significant voxel-wise functional connections that exist between the two regions. For a given voxel-pair, a significant connection indicates that the activities recorded at the two voxels are highly synchronized (have high correlation value). Thus, a great change in the number of statistically significant voxel-wise functional connections across time for a ROI-pair corresponds to a great change in the level of their synchronization (i.e., it represents significant plasticity). We then examine if this way of defining functional connectivity helps us to detect any plasticity effects in the TBI and the HC groups across two time-points when BOLD-fMRI scans were conducted (see Section Subjects). While we hypothesize that this approach will be able to detect some plasticity, we also hypothesize that a more sensitive procedure will detect greater plasticity.

3. Analysis-3: Next, we document the drawbacks of the simple voxel-based approach (defined in Analysis-2) to examine plasticity, and introduce the EC-based approach which is also built upon the basic idea of counting statistically valid voxel-wise connections between regions. However, the EC based approach makes use of an adaptive algorithm that probes for sub-regional-pairs that show significant decrease/increase in statistically valid voxel-wise connections across two time-points when BOLD-fMRI scans were conducted (see Section Subjects). Thus, the EC-based approach detects sub-regional-pairs within a ROI-pair that exhibit significant plasticity (the plasticity of each sub-regional-pair is either positive or negative). Using these sub-regional-pairs, an overall estimate of positive and negative plasticity will be derived for the given ROI-pair.

The EC approach involves some stochasticity, and therefore, it is very important to examine its performance across multiple independent runs. Since the primary goal of this article is to use the EC procedure to estimate the level of positive and negative plasticity between a ROI-pair, we document its performance by measuring the consistency of the plasticity estimates found across multiple independent runs. As a secondary measure of consistency, we also document the spatial overlap between the sub-regional-pairs found across independent runs.

Materials and Methods

Data

Subjects

The study included 26 subjects, 14 individuals with moderate and severe TBI and 12 age and education matched healthy control participants (HCs; see Table 1 for demographic information). Individuals with TBI completed two separate MRI sessions ~3 and 6 months following the resolution of post-traumatic amnesia. HCs completed two separate MRI sessions separated by ~3–4 months to provide context for natural variation in healthy brain networks (see Table 1). For the purposes of this article, the two time-point fMRI-scan sessions will be referred to as session-1 and session-2. For the TBI sample, each testing session included both MRI data collection and cognitive testing. Moderate and severe TBI were defined using Glasgow Coma Scale (GCS) at the time of injury, whereby a score of 3–8 was indicative of severe injury and a GCS of 9–12 was indicative of moderate injury (Teasdale and Jennett, 1974). Two subjects who received a GCS of 14 and one subject who received a 24-h GCS of 15 were included because of either mental status change after the initial GCS scoring or the presence of significant brain injury documented through medical records. Subjects were excluded if they were receiving treatment for concomitant injuries (e.g., orthopedic injuries or injury to the spinal cord) that would make it difficult for them to remain still and comfortable in the MRI environment. For the purposes of this study where we aim to examine functional plasticity, we use a well-documented clinical sample with known network changes from 3–6 months post injury (Hillary et al., 2014; Rajtmajer et al., 2015).

TABLE 1

Table 1. Subject demographic information.

Research was approved by the Institutional Review Board and the Office of Research Protections at the Pennsylvania State University. Individuals included in the study demonstrated some level of cognitive impairment. If an individual retained the ability to sign medical documents and/or function independently, then consent was accepted; if the individual was not functionally independent, then a caregiver's signature was required in addition to the subject's signature of assent.

Functional Data

Study participants were scanned using one of three MRI machines, including a Philips Achieva 3T scanner in the Department of Radiology at Hershey Medical Center, Hershey, PA and two identical Siemens Magnetom Trio 3T scanners (Social, Life, and Engineering Sciences Imaging Center at the Pennsylvania State University in University Park, PA; Department of Radiology at Hershey Medical Center in Hershey, PA). Of the subjects, 7 were scanned on the Achieva 3T scanner, 5 subjects were scanned on the Magnetom Trio 3T scanner at the Pennsylvania State University, and 14 subjects were scanned on the Magnetom Trio 3T scanner at Hershey Medical Center. For repeat scanning, all subjects were scanned on the same machine across time points. Prior to scanning subjects were made aware of the importance of minimizing motion within the scanner which at times included pre-scanning preparation in the mock-scanning environment.

Data Acquisition Parameters

Anatomical images with a spatial resolution of 1.0 × 1.0 × 1.0 mm were acquired using an MPRAGE sequence: 2000 ms/2.03 ms/9° (repetition time (TR)/echo time (TE)/flip angle (FA), 256 × 256 mm² field of view (FOV), and 256 × 256 acquisition matrix with 1 mm slices. Echo planar imaging (EPI) was used to examine the blood oxygen level dependent response for functioning imaging. Imaging parameters for EPI were 2000 ms/30 ms/90° (TR/TE/FA), 240 × 240 mm² FOV, and 80 × 80 acquisition matrix with 4 mm slices.

Data Preprocessing

Resting data were collected over the course of a 5 min period, resulting in 150 volumes of data. The first five volumes were removed from analyses to control for signal instability, resulting in a time series of 145 volumes. For all volumes bad slices were first repaired using the art-slice procedure, which is part of the ArtRepair toolbox (Mazaika et al., 2009). The volumes were then slice time corrected and realigned to the first volume using SPM8. Even after realignment, there remain important influences of motion. As recommended in the ArtRepair manual, this can be handled via Volume repair at the end of the preprocessing pipeline. Power et al. (2012) demonstrate that global signal regression is a powerful method for reducing the effects of motion. Given concerns published elsewhere (Murphy et al., 2009; Saad et al., 2012) we do not use this approach. Spike artifacts were eliminated using the despike filter available in the ArtRepair toolbox. Each subject's high-resolution (1 × 1 × 1 mm) TI image was co-registered to the mean functional image using SPM8. The co-registered T1 image was then segmented using SPM8, which produced a normalized (MNI space) gray matter image. The functional images were then normalized to the MNI space. The normalized functional volumes and the gray matter image were resliced, and the voxel dimensions were defined as 3 × 3 × 3 mm. The equivalence of voxel size allowed one to one mapping of voxels between a normalized functional and normalized T1 volumes, and provided a mechanism to extract gray matter signal (see the following section). To reduce the effect of ringing artifact (Lindquist and Wager, 2008) and to improve the signal to noise ratio, a 6 mm smoothing filter was applied. The CSF and the white matter nuisance signals were regressed out, and a bandpass filter of 0.01 to 0.12 Hz (Bassett et al., 2011; Rajtmajer et al., 2015) was applied using the CONN toolbox (Whitfield-Gabrieli and Nieto-Castanon, 2012). The removal of nuisance signal due to physiology and head motion is a developing issue receiving significant recent attention in the resting MRI literature (Power et al., 2012, 2014). Finally, the bad volumes (volumes with excessive motion related problems) were repaired using ArtRepair software (Mazaika et al., 2009).

Binary Gray Matter Mask

We developed an average gray matter image using the normalized/resliced gray matter images of the control subjects at both time-points (12 gray matter images × 2 scan sessions = 24 gray matter images). Since the segmentation procedure assigns a probabilistic value to each voxel of the gray matter segmented file (i.e., the probability that the voxel belongs to the gray matter), we used only the voxels with probability value greater than 0.6 (Zalesky et al., 2012) to define a group averaged binary gray matter mask.

Defining ROIs

To determine the gray matter voxels that are part of DMN (left/right) and ECN (left/right), we first co-registered Automated Anatomical Labeling (AAL) atlas (Tzourio-Mazoyer et al., 2002) to the group averaged binary gray matter mask (See Section Binary Gray Matter Mask). This allowed us to label each gray matter voxel with a unique AAL identification number. Since the binary gray matter mask and the functional images are composed of voxels of the same dimensions and are normalized to MNI-space (See Section Data Preprocessing), the voxels have one to one correspondence between these two image modalities. Hence, the ROIs are defined using the gray matter voxels with AAL identification numbers that represent the brain regions that are part of DMN (left/right) and ECN (left/right). The brain regions constituting these ROIs are presented in Table S5.

Developing Subject-Specific Connectivity Matrices

As discussed in the introduction, both analysis-2 (simple voxel-based approach) and analysis-3 (EC based approach) are based on counting statistically significant voxel-wise connections between the ROI-pairs. Thus given a ROI-pair A-B, we first create two voxel-level connectivity matrices for each subject: one for session-1 (M₁) and one for session-2 (M₂). For each session, the connectivity matrices are developed by performing a Pearson's correlation between the set of signals that are recorded from the voxels that are part of ROI-A to the set of signals that are recorded from the voxels that are part of ROI-B. For example, if ROI-A consists of 500 voxels and ROI-B consists of 1800 voxels, then the connectivity matrix will have a dimension of 500 × 1800. The correlation values that fail false discovery rate (FDR) test at a threshold of 0.05 are discarded; these cells in the connectivity matrices are set to 0.

For illustration purposes, for both analysis-2 and analysis-3 (see Introduction) we only consider positively correlated functional edges as we are interested in examining the change in positive synchronization between the four ROI-pairs. Therefore, we discard the cells that represent negative correlations (these cells are set to 0 as well). However, both the simple voxel-based approach and the EC process are not dependent on this choice, and an investigator may choose to use both positive and negative correlated edges. The remaining cells were set to 1. Since we binarized the voxel-wise connections, while computing the change in connectivity strength of a sub-regional-pair across two fMRI-scan sessions we only need to compute the total number of connections at session-1 and session-2 instead of storing two separate distributions of statistically significant (FDR corrected; See the paragraph above) voxel-wise correlation values for each session. This makes the process computationally very efficient (less memory required). Regardless of the memory requirement, if the session-1 and session-2 distributions of statistically significant voxel-wise correlation values are used for examining plasticity between a given sub-regional-pair across time, then additional statistical correction mechanisms may have to be implemented as the number of data-points in the distributions may not be the same. For example, in an extreme case, for one or both of the sessions the distribution may contain no data-points if there are no significant voxel-wise connections. Of course, it remains a possibility that the connection strengths may be used to examine plasticity, but for the sake of simplicity we chose to binarize the connectivity matrix. Even though we binarize and lose some information, the approach discussed in this article provides novel insights into connectivity changes which otherwise would be hard to detect with an averaging based approach.

Evolutionary Computation Based Approach to Detect Significantly Plastic Sub-Regional Pairs

Given an investigator defined ROI-pair, ROI-A and ROI-B (henceforth, we will refer to this ROI pair as A-B), the goal of the evolutionary search based procedure is to find significantly (see Section Objective Function) plastic sub-regional-pairs within A-B by simultaneously utilizing BOLD-fMRI data collected from two fMRI-scan sessions that are separated by a finite duration of time; in this article for both TBIs and HCs the sessions are separated by 3 months (See Section Subjects). Further, this procedure is designed such that it can be run for each subject separately.

The procedure is schematically shown in Figure 2 and can be summarized as follows. For a given ROI-pair, A-B, we search for significantly plastic sub-regional-pairs in a recursive manner; at each recursion level an EC based search is conducted for a significantly plastic sub-regional-pair. At recursion level-1, if a significantly plastic sub-regional-pair is found then all the edges between the detected sub-regional-pair are stored in a file, and then discarded from the region (from both session-1 and session-2 data) so that the EC does not converge to this solution again and searches for another significantly plastic sub-regional-pair at the next recursion level. At recursion level-2, the evolutionary search procedure is applied to the remaining edges; the process repeats (recursion level-3, recursion level-4,….) until no significantly plastic sub-regional-pair can be found. That is, we terminate the entire procedure at the recursion level at which the EC based search fails to find a significantly plastic sub-regional-pair. As an example, if the procedure is run such that session-1 and session-2 data are intentionally made identical, then the EC based search fails to find a significantly plastic sub-regional-pair (as there are none) at recursion level-1, and hence the entire procedure terminates at recursion level-1. Using the set of significantly plastic sub-regional-pairs detected by the EC based procedure, we then develop estimates for positive and negative plasticity (see Section Estimating Positive and Negative Plasticity) for the ROI-pair. In the following paragraph, we discuss the EC-based search that is conducted at each recursion level.

Evolutionary algorithms are designed to mimic biological evolutionary processes to solve multi-modal and multi-dimensional optimization problems (Rechenberg, 1971; Schwefel, 1974; De Jong, 2006). In evolutionary computation (EC), a candidate solution (which is analogous to biological phenotype) for a given optimization task is encoded as an array of numbers (which is analogous to biological genotype/chromosome). If the array of numbers (the chromosome) directly represents the candidate solution, then it is referred to as a direct mapping, otherwise it is referred to as an indirect mapping (Gruau and Whitley, 1993; Gruau, 1994; Kobayashi and Ohbayashi, 1999; Roy et al., 2013).

For our purposes, given an investigator defined ROI-pair, ROI-A and ROI-B, a candidate solution represents one sub-region from each ROI. Henceforth, we will refer to this ROI pair as A-B. Each sub-region is composed of a connected set of voxels. Two voxels are said to be connected (or neighbors) if they share a common physical boundary, i.e., there is a path of length 1 between these two voxels. A set of voxels are said to be connected if there is a path between every pair of voxels in the set. The chromosome will not encode these voxels directly; instead, a compression scheme (see section Chromosome Encoding below) will be used to define a sub-region within each ROI. Thus we use indirect mapping in this work.

At a given recursion level, an evolutionary search for a significantly plastic sub-regional-pair will begin with a randomly generated population of chromosomes, each representing a sub-regional-pair within A-B. Each chromosome will then produce an offspring chromosome using an asexual (single parent) reproduction operator (Fogel, 1999, 2003; Davoian et al., 2006) (See Section Reproduction Operator for details). The reproduction operator ensures that an offspring shall encode a sub-regional-pair that has some spatial similarity (many common voxels) to the sub-regional-pair encoded by the parent (see section Chromosome Encoding); this process ensures that the search is not random. The fitness value (or quality of each solution) of a chromosome (parent and offspring) will be evaluated based on the level of plasticity (see objective function below) that the sub-regional-pair exhibits across 2 fMRI-scan sessions. That is, if the number of voxel-wise functional connections between the sub-regional-pair changes significantly across two sessions (session-1 and session-2), then the sub-regional-pair is highly plastic (see Figure 4). To ensure that the best solutions are preserved, a population-elitist selection (Eshelman, 1991; Eshelman and Schaffer, 1991) procedure will be implemented in each generation to choose the fittest chromosomes (parents and offspring combined) for the next generation; the population size in every generation will remain constant (see Section Evolutionary Search Parameters and Figure 2). That is, in population-elitist selection all the offspring and parents are pooled together and top 50% best solutions are chosen. As a result of this, even the fitness landscape near poor solutions gets the chance to be explored across generations. Hence, this form of selection maintains good diversity in the population, and therefore, provides a mechanism to prevent the evolutionary search from premature convergence (Eshelman and Schaffer, 1991). The reproduction-process, fitness-evaluation, and population-elitist selection steps will be repeated until the stopping criteria (see Section Evolutionary Search Parameters) are met (see Figures 1, 2).

FIGURE 1

Figure 1. The sequence of steps involved in the evolutionary computation based procedure.

FIGURE 2

Figure 2. The figure illustrates the steps involved in the evolutionary computation (EC) based procedure to examine plasticity between an investigator defined ROI-pair (shown as the oval shaped white regions on the brain image), ROI-1 and ROI-2. (A) The procedure begins at recursion level-1 where the evolutionary search for a significantly plastic sub-regional-pair (SRP) is conducted. In this illustration, the population size (number of parental chromosomes) is 3, and each parental chromosome (P1, P2, and P3) produces 3 offspring (C1, C2, and C3). Based on population-elitist selection the top 3 chromosomes are chosen for the next generation, and relabeled as P1, P2, and P3 in a descending order of fitness. The color of each chromosome represents its fitness (|Z-score|); the higher the fitness of the chromosome, greater is the change in the connectivity strength between the sub-ROIs that the chromosome encodes. At the end of recursion level-1, a significantly plastic SRP is detected (SRP-1; shown in dark-red) and the chromosome encoding this SRP is stored in a file. SRP-1 is then added to the list of blocked SRPs, and therefore, at the beginning of recursion level-2 the list contains one blocked SRP; EC will not converge to this solution at recursion level-2. At recursion level-2 a new search begins, and SRP-2 (show in light-red) is detected. For recursion level-3, SRP-1 and SRP-2 are blocked. The recursive process stops at recursion level-K when no significantly plastic SRP is found. (B) To block a SRP all the voxel-wise edges between the SRP are deleted from connectivity matrices (See Section Developing Subject-Specific Connectivity Matrices), M1 and M2 (shown as colored rectangular patches). (C) We illustrate an example consisting of 3 significantly plastic SRPs to show how plasticity estimates are derived.

The basic idea behind this evolutionary search can be explained as follows. Since the fittest candidate solutions (i.e., fittest chromosomes) survive to reproduce, the space around the sub-regional-pairs encoded by those fittest candidates are searched thoroughly over generations as a result of the reproduction operation, thereby increasing the likelihood of discovering significantly plastic sub-regional-pairs. This entire procedure is referred to as “evolutionary” because it holds similarity to natural evolution where the phenotypes that are suited for survival propagate along the generations and eventually become prevalent in the population; likewise our evolutionary search eventually converges to sub-regional-pairs (the artificial phenotype) that represent strong plasticity. Below, we discuss the details of this procedure.

Chromosome Encoding

A chromosome indirectly encodes a sub-regional-pair within A-B (see Figure 3B). We will explain the chromosome encoding using a specific example shown in Figure 3A. In Figure 3A, the voxels 0–7 in region-A (the left grid) and voxels 0–7 in region-B (the right grid) represent a ROI pair A-B; for the purpose of clarity we have not shown the edges that exist between these regions. The thick black boundaries represent the smallest cuboids, cuboid-A and cuboid-B, that respectively encapsulate the regions A and B. For the ease of explanation, in Figure 3A we have considered a special case where the ROI pair A-B has no width in the Z-direction (the direction pointing out of the page), and therefore, the cuboid encapsulating them appear as a rectangle.

FIGURE 3

Figure 3. (A) The figure illustrates the chromosome structure used for the evolutionary-search. The voxels 0–7 in the left grid and the voxels 0–7 in the right grid represent a ROI pair. The chromosomal-variables X1, Y1, and Z1 encode a point (small black-dot) within the smallest cuboid (the thick rectangular boundary in left grid) encapsulating ROI-A, and X2, Y2, and Z2 encode a point (small black-dot) within the smallest cuboid (the thick rectangular boundary in right grid) encapsulating ROI-B. The closest gray matter voxels to these points are chosen as the root voxels (i.e., voxel 0 in both ROIs). The circles with numbers represent gray matter voxels. The circles are sequentially numbered using BFS strategy w.r.t the root voxel; the voxels with lowest z-axis values are numbered first, followed by the voxels with lowest y-axis values, and then lowest x-axis values. L1 and L2 define the number of voxels chosen around the root node from ROI-A and ROI-B using a BFS strategy to form the sub-regional-pair; the sub-regional-pair also includes the root node. (B) We illustrate an example of a sub-regional-pair where L1 = 3 and L2 = 4. The voxels within the orange boundaries represent the sub-regional-pair.

The objective of the evolutionary search is to find a pair of sub-regions (i.e., a sub-region within cuboid-A and a sub-region within cuboid-B) that qualifies as significantly (see Section Objective Function) plastic where each sub-region is composed of a set of connected voxels. Two voxels are said to be spatially connected (neighbors) if they share a common boundary: a cubic-voxel can have 6 possible neighbors. The first 3 locations of the chromosome (X1, Y1, and Z1) encode the coordinate of a point (shown as a small black-dot in Figure 3A) inside cuboid-A. The voxel that is closest to this point and part of the ROI is chosen as the root node (the black node in Region-A, i.e., voxel 0), and the fourth component of the chromosome (an integer L1) determines the number of neighbors that are selected for membership in the sub-region. The neighbors are chosen using a breadth first search (BFS) strategy. The variables X2, Y2, Z2, and L2 (see Figure 3A) encode a sub-region within cuboid-B in an analogous way. We note that while we use this BFS approach in this work, alternative approaches for choosing the neighboring voxels exist. The motivation for choosing the BFS strategy in this work is because it imposes the constraint that the sub-regional-pair must be composed of a connected set of voxels.

In some cases the BFS search may identify multiple candidate voxels for inclusion. In these cases we arbitrarily prioritize voxels with the smallest z-axis value, then similarly with the y-axis and then x-axis values. For example, in Figure 3B, if L1 = 3, then the voxels 1, 2, 3, and the root node (voxel 0) represent the sub-region within ROI-A. While this choice influences the shapes of candidate sub-regional-pairs, the shape of the solutions carries little meaning in the absence of structural information. Rather, the location and general size of the significantly plastic sub-regional-pairs is of particular interest specifically because this information may provide structural insight.

Reproduction Operator

In the context of the current work, we require a reproduction operator that allows the offspring chromosomes to encode sub-regional-pairs that have some spatial similarity to the sub-regional-pairs encoded by the parental-chromosomes. As many operators can be developed to accommodate this requirement, we decided to develop a simple mutation based reproduction operator (i.e., a single parent reproduction operator, See Figure 2) that allows an investigator to easily control the level of similarity between a parent and an offspring chromosome; this operator is similar to mutation based reproduction operators used in evolutionary programming (Fogel, 1999, 2003; Davoian et al., 2006). If offspring are made to be too similar to the parents, then the search is less exploratory. On the other hand, if the offspring are too different from the parents then the search becomes quite random. An appropriate level of exploration is often problem-specific, and it is generally determined by adjusting the exploration level till consistent good solutions are discovered across multiple independent runs. We have explained the reproduction operator below.

A parental-chromosome (of size 8 as discussed above) is given as follows:

Parent = {X1, Y1, Z1, L1, X2, Y2, Z2, L2}

An offspring will be generated by producing a random number for each of the eight positions within a bound defined by the offset-parameters {dx1, dy1, dx1, dl1, dx2, dy2, dz2, dl2}. For example, the value of X1 for a child, X1_c, will be determined by the value of X1 for its parent, X1_p, such that it falls in the range [X1_p−dx1, X1_p+dx1].

Here, we note that even though the evolutionary search begins with a population of randomly generated chromosomes encoding various combinations of sub-regional-pairs, the reproduction operation along with the population-elitist-selection procedure causes the search trajectory to tend toward sub-regional-pairs that are highly plastic.

Objective Function

Given a candidate solution (a sub-regional-pair encoded by a chromosome), the change in the functional relationship (plasticity) that occurs between 2 time points (i.e., 2 scan sessions) is defined using the number of connections that exist between the sub-regional-pair at session-1 (NC₁), the number of connections that exist between them at session-2 (NC₂), and the total possible connections (TC) that can exist between them. The values NC₁ and NC₂ can be easily derived by counting the total number of cells that have the value 1 in the connectivity matrices, M₁ and M₂ (see Section Developing Subject-Specific Connectivity Matrices), respectively.

Since the connections are binary (they exist or they do not), the objective function is defined as the Z-score of a binomial test performed using the fraction of connections between the candidate sub-regional-pairs existing at session-1 and session-2. The objective function is defined as follows.

$\begin{array}{rc}{P}_{\mathit{\text{1}}}=N{C}_{\mathit{\text{1}}}∕TC& (1)\end{array}$

Where P₁is the fraction of connections that exist between the sub-regional-pair at session-1.

$\begin{array}{rc}E{C}_{\mathit{\text{2}}}=P_{\mathit{\text{1}}}^{*}TC& (2)\end{array}$

Hence,

$\begin{array}{rc}E{C}_{\mathit{\text{2}}}=N{C}_{\mathit{\text{1}}}& (3)\end{array}$

where EC₂ is the expected number of connections at session-2 between the sub-regional-pair established using session-1 data.

The standard deviation for the binomial distribution established using session-1 data will be defined as:

$\begin{array}{rc}SD=Sqrt\left[T{C}^{*}{P}_{\mathit{\text{1}}}^{*}(1-P_{\mathit{\text{1}}})\right]& (4)\end{array}$

The fitness of the candidate solution will then be defined as the absolute (signs neglected) Z-score of the session-2 connection based on the binomial distribution established using session-1 data. That is,

$\begin{array}{rc}Fitness=Z-score=|(N{C}_{\mathit{\text{2}}}-E{C}_{\mathit{\text{2}}})∕SD|& (5)\end{array}$

Thus a candidate solution with higher absolute Z-score will be more fit than a candidate solution with a lower absolute Z-score. That is, a highly fit candidate solution would represent sub-regional-pair with high level of plasticity (i.e., such sub-regional-pairs are significantly plastic).

Note that the fitness function disregards the sign of the Z-score, as the purpose of the evolutionary search is to simply find significantly plastic sub-regional-pairs while disregarding the direction of plasticity. However, once the evolutionary process finishes running and locates the significantly plastic sub-regional-pairs, we run a simple post hoc analysis to determine the direction of the change according to the sign of the Z-score. This way we are able to detect both positively and negatively plastic sub-regional-pairs (Figure 4).

FIGURE 4

Figure 4. The figure illustrates an example of a positively and negatively plastic sub-regional-pair. In this figure, session-1 and session-2 refers to two fMRI-scan sessions that are separated by finite duration of time; chronologically, session-2 occurs after session-1. Edges indicate a voxel from ROI-A is significantly correlated with a voxel from ROI-B. The red sub-regional-pair gains functional edges from session-1 to session-2 (hence positive plastic), while the blue sub-regional-pair shows decline in the number of functional edges (hence negative plastic). The binomial test based objective functional detects such plastic sub-regional-pairs.

Evolutionary Search Parameters

A. Criterion for stopping the search for any new significantly plastic sub-regional-pair For a given investigator defined ROI-pair, A-B, the evolutionary search is run recursively until no new significantly plastic sub-regional pairs can be found. That is at a given recursion level, if an evolutionary search produces a solution with the objective function value less than 1 (i.e., Z-score less than 1), then the evolutionary search terminates, and no new search for a significantly plastic sub-regional-pair is initiated. Typically the EC based search starts by detecting plastic sub-regional-pairs with very high Z-scores (>>1). Therefore, for all practical purposes, when the evolutionary search starts producing solutions with a Z-score close to 1 (or less), most highly plastic sub-regional-pairs have already been detected. Of course, there is a possibility that the evolutionary search at a given recursion level may have converged to a poor solution. To confirm that no more good solutions (significantly plastic sub-regional-pairs) exist, either multiple runs should be executed or the search should be run with a sufficiently large population size. For the current work we chose the latter approach. The population size and other evolutionary parameters are discussed below.

B. Criteria for stopping the evolutionary search At a given recursion level, the evolutionary search is terminated when one of the following occurs:

1. Complete convergence (all chromosome encode the same solution).

2. The objective function value does not improve over a span of 100 generations. For a more rigorous search, one may choose this value to be much higher. In context of the present study, this value provided a good balance between runtime and ability of the procedure to detect statistically significant plastic sub-regional-pairs.

C. Evolutionary search parameters Below, we summarize the evolutionary search parameters:

1. Population size is set to 400.

2. Reproduction parameters, dx1, dy1, dx1, dx2, dy2, and dz2, were set to 6 mm; the parameters dl1 and dl2, were set to 4. This forced the mutation operator to choose a root-voxel for an offspring not more than two voxels away from the parent root-voxel (the voxel size of preprocessed functional image is 3 × 3 × 3 mm ^3). As above, these values were chosen to provide a balance between consistent good solutions and run time; automatic tuning of these parameters could in principle be used to further optimize this approach.

3. Given a ROI-pair A-B, in our experiments we allowed the chromosomal-parameters x1, y1, z1 to point to any spatial location within the smallest cube encapsulating ROI A, and x2, y2, z2 to point to any spatial location within the smallest cube encapsulating ROI B. The chromosomal-parameters, L1 and L2 were allowed to take any value at a step size of 5 between 64 (4 × 4 × 4 voxels) to the total number of voxels in ROI A and ROI B, respectively. Both the DMN regions consisted of ~500 voxels, and each ECN region consisted of about 1800 voxels.

Post hoc Correction for Multiple Tests

In reference to the objective function discussed above, a set of significantly plastic sub-regional-pairs detected using EC represents an array of statistically significant Binomial tests. Therefore, we perform Bonferroni correction to establish corrected p-values for the sub-regional-pairs; only the sub-regional-pairs with corrected p-values less than 0.05 are considered valid. In this article we chose Bonferroni correction because we wanted to be more conservative in our estimate of the plasticity levels of individuals with TBI and controls to examine robustly if there was heightened or diminished connectivity in TBI cases as compared to controls. Depending on context, other (e.g., data-driven) approaches may be appropriate; see for instance Derrac et al. (2011).

Estimating Positive and Negative Plasticity

Using the set of significantly plastic sub-regional-pairs that qualify post-hoc correction (see Section Post hoc Correction for Multiple Tests), we evaluate the percentage of total voxel-wise connections that contribute to positive and negative plasticity. For example, as shown in Figure 2 (Box C), let us consider an investigator defined ROI-pair, ROI-1, and ROI-2, where ROI-1 is composed of 100 voxels and ROI-2 is composed of 500 voxels. Therefore, at most there can be 50000 voxel-wise connections between ROI-1 and ROI-2. Let us also consider that using the EC based procedure (including post-hoc correction) we detected 3 significantly plastic sub-regional-pairs: SRP-1, SRP-2, and SRP-3. SRP-1 and SRP-3 were positively plastic and gained 10000 and 1000 voxel-wise connections from session-1 to session-2, respectively, while SRP-2 lost 5000 connections from session-1 to session-2. Thus, the percentage voxel-wise connections contributing to positive plasticity is [(10000+1000)*100]/50000% = 22% and percentage voxel-wise connections contributing to negative plasticity is [5000*100]/50000% = 10%.

Results

Averaging Based Approach

To examine the effectiveness of the averaging based approach to detect plasticity, for each subject for each session we developed an average signal for all the four ROIs: DMN-frontal left/right and ECN-frontal left/right. Then we established the functional connectivity strength between all the four ROI-pairs for each subject for each session using Pearson's correlation coefficient. We performed a paired t-test of Fisher z-transformed correlation values to examine if this approach can detect group level plasticity (change in functional connectivity strength between session-1 and session-2) for the TBIs and the HCs. As shown in Figures 5A,B this approach failed to detect plasticity. Only the functional relationship between DMN-frontal right and ECN-frontal right in HCs was found to be significant. This supports our hypothesis that the averaging based approach fails to detect plasticity.

FIGURE 5

Figure 5. DFL, DFR, EFL, and EFR, stand for default mode frontal-left, default mode frontal-right, executive frontal-left, and executive frontal-right, respectively. (A) The box plot illustrates the distribution of the connectivity strength (correlation between ROI-wise averaged signal) of the TBI-subjects at session-1 (red) and session-2 (blue) for all the four ROI-pairs. (B) The box plot illustrates the distribution of the connectivity strength (correlation between ROI-wise averaged signal) of the control cases at session-1 (red) and session-2 (blue) for all the four ROI-pairs. Using paired t-test (after Fisher z-transforming the correlation values), the difference between the group level means at 2 sessions was found to be insignificant for the TBIs. For the controls, except for the fourth ROI-pair, plasticity for rest all ROI-pairs was insignificant.

Simple Voxel-Level Approach

In order to demonstrate that a voxel-based approach will be able to detect plasticity better than an averaging based approach, we developed a simple edge count based procedure. That is, for each subject for each session, a connectivity matrix was developed for each ROI-pair (see Section Developing Subject-Specific Connectivity Matrices developing subject specific connectivity matrix). The total number of functional edges between a pair of ROIs for each subject at each session was evaluated by considering the connectivity matrices entry by entry and counting the cells that have the value 1. Using the total edge count we performed paired Wilcoxon signed-rank test for all the four ROI-pairs for HC and TBI group (see Figure 6); we chose a non-parametric test (Wilcoxon test) because based on Shapiro–Wilk tests we found the data was not consistently normal across all 4 ROI-pairs and across 2 time-points. As shown in Figure 6B, no plasticity was detected for the controls. For the TBI cases, the tests were insignificant at p = 0.05. However, at p = 0.1, three of the four ROI-pairs showed some signs of plasticity (see Figure 6A). This supports our hypothesis that a voxel-level procedure will be better able to capture plasticity than averaging based approaches. However, we anticipate that this simple voxel-level procedure may be underestimating the level of plasticity as it does not differentiate between positive and negative plastic sub-regional-pairs.

FIGURE 6

Figure 6. DFL, DFR, EFL, and EFR, stand for default mode frontal-left, default mode frontal-right, executive frontal-left, and executive frontal-right, respectively. Connection strength is defined as the number of voxel-wise functional edges between the ROI-pairs. (A) The box plot illustrates the distribution of connection strength of TBI-subjects at session-1 (red) and session-2 (blue) for all the four ROI-pairs. (B) The box plot illustrates the distribution of connection strength of control cases at session-1 (red) and session-2 (blue) for all the four ROI-pairs. Using paired Wilcoxon signed-rank test, the difference between the group level means at 2 sessions was found to be insignificant for TBIs and controls for all the 4 ROI-pairs. However, in contrast to averaging analysis (see Figure 5) there seem to be some signs of plasticity for TBI at a p-value of 0.1.

Evolutionary Computation Based Approach

Examining Functional Plasticity

In Figure 7, we have shown the percentage of total voxel-wise connections that contributed to positive and negative plasticity for each TBI-subject and the HCs between the ROI pair, DMN-frontal right and ECN frontal right. As illustrated in the figure the mean positive plasticity level of the TBI was found to be significantly higher than the HCs (p = 0.009). A similar trend was observed for the ROI-pairs, DMN-frontal left and ECN frontal left (p = 0.046), and DMN-frontal right and ECN frontal left (p = 0.06); the group-level comparison of the mean positive plasticity for each ROI-pair was conducted using Wilcoxon rank sum test. We chose a non-parametric test (Wilcoxon test) because based on Shapiro–Wilk tests we found the data was not consistently normal across all 4 ROI-pairs.

FIGURE 7

Figure 7. The percentage of total voxel-wise connections between DMN-frontal right and ECN-frontal right that contributed to statistically significant positive (red) and negative (blue) plasticity. A similar trend was observed for the ROI-pairs, DMN-frontal left and ECN frontal left, and DMN-frontal right and ECN frontal left. (A) The bar-graph illustrates the percentage of total voxel-wise edges that contributed to positive (red) and negative (blue) plasticity for each TBI subject. (B) The bar-graph illustrates the percentage of total voxel-wise edges that contributed to positive (red) and negative (blue) plasticity for each HC subject. The red and the blue horizontal lines indicate the mean positive and negative plasticity level across a group, respectively. The mean positive plasticity for the TBI was significantly higher than the controls.

This finding is consistent with previous findings that suggested such heightened functional connectivity in large scale networks in TBI subjects (Nakamura et al., 2009; Bonnelle et al., 2011; Hillary et al., 2011a; Caeyenberghs et al., 2012). Figure 7 thus suggests that while there is significant amount of positive and negative functional plasticity in the HC brain, the TBIs show more positive plasticity across time. Knowing the positive and negative plasticity level for HCs can be very useful for establishing thresholds for defining a healthy brain, and hence can be used for developing diagnostic systems for detecting neurological illnesses/insults.

To examine the consistency of the EC procedure in estimating plasticity across multiple independent runs, we ran 10 experiments per subject per ROI-pair. We then evaluated the standard deviation of the plasticity estimates (% voxel-wise connections contributing to positive and negative plasticity) for each subject. The standard deviation of the plasticity estimate obtained across multiple runs was very low (typically < 0.5%), thus highlighting the consistency of this approach (see Tables S1–S4).

The EC procedure was run on a serial mode (processor information: Intel(R) Xeon(R) CPU E5-2670 0 @ 2.60 GHz, 2 GB memory), and we developed an estimate of the runtime for each ROI-pair by averaging across 260 experiments (10 experiment × 26 subjects); the runtime results are presented in Table S6. The average runtime was ~15 min for four ROI-pairs. Because EC is implicitly parallelizable, the runtime can be significantly reduced using parallel programming. The EC procedure was performed in the MNI space. While the EC procedure can be run in the native space, for our purposes, normalization was required for 2 reasons. First for examining functional relationship between regions, it required extraction of all gray matter voxels in the brain and most segmentation algorithms perform some degree of normalization (e.g., SPM8) and hence implicitly require re-sampling. Second, in order to examine connectivity changes at separate time points simultaneously, identification of equivalent regions was necessary and this was achieved via normalization. Since the voxel size of the functional image before (3 × 3 × 4 mm) and after (3 × 3 × 3 mm) normalization was very similar we anticipate that the computational runtime presented in Table S6 should remain roughly the same.

In Figure 8, we have shown some significantly plastic sub-regional-pairs that were detected between DMN-frontal right to ECN-frontal left for TBI subject-11; for illustration purposes we have only shown the first 5 highly plastic (out of 41) sub-regional-pairs detected by the EC approach. The figure illustrates that sub-regional-pairs contributing to plasticity may be much smaller than the ROI-pair. The figure also illustrates that within a ROI-pair there may be both positively plastic sub-regional-pairs and negatively plastic sub-regional-pairs.

FIGURE 8

Figure 8. An illustration of the first 5 significantly plastic sub-regional-pairs that were detected by EC based procedure between DMN-frontal left to ECN-frontal right for TBI subject-11; the DMN and the ECN regions have been colored purple. The plastic sub-regional-pairs are colored yellow, green, pink, blue, and light blue (from left to right). The sub-regional-pairs are much smaller than the ROI-pair, and therefore, it may be hard to detect such subtle plasticity effects using a conventional averaging based approach. The +∕−ve signs indicate the plasticity direction.

Overlap Analysis

In order to test if the evolutionary process detects the same set of significantly plastic sub-regional pairs across multiple independent runs, per subject per ROI-pair the overlap analysis was conducted for 10 pair of independent runs.

The overlap between two sub-regional-pairs (a sub-regional-pair from run-1 and a sub-regional-pair from run-2) was defined using the number of common voxels they share; the overlap score was defined using the Sorenson-dice index (Sorenson, 1948). The Sorenson-dice index ranges from 0 to 1, with 0 indicating no overlap and 1 indicating complete overlap. For a given pair of independent runs, run-1 and run-2, the maximum spatial overlap score for each significantly plastic sub-regional-pair detected in run-1 was evaluated by comparing them with the significantly plastic sub-regional-pairs detected in run-2. In this way each sub-regional-pair from run-1 received an overlap score. Similarly, the maximum spatial overlap score for each sub-regional-pair detected in run-2 was evaluated by comparing them with run-1 sub-regional-pairs. The mean of the maximal overlap scores of all sub-regional-pairs (run-1 and run-2 included) was used for defining the overall overlap between the 2 runs. This procedure was repeated for all 10 pair of runs, per subject per ROI-pair and then the mean and standard deviation of the overlap score was established. As shown in Table 2A the overlap score was about 0.8. Also, we found no significant difference (t-test and Wilcoxon rank sum test, p > 0.05) between the Sorenson scores of the TBI and HC.

TABLE 2A

Table 2A. Sub-cluster-pair overlap analysis based on Sorenson score.

We also performed overlap analysis using the adjusted rand index (ARI; Vinh et al., 2009). In this context, the ARI considers the set of voxel-pairs between the ROI-pair and measures the overlap between 2 different sets of labelings (produced by 2 independent EC runs), which we denote T and U. The set T has (N_T + 1) labels representing N_T sub-regional-pairs and a label for voxel-pairs that are not part of any significantly plastic sub-regional pair. Similarly, the set U contains (N_U + 1) labels. Here, each voxel-pair that exists in one of the N_T significantly plastic sub-regional-pair is assigned a corresponding value from {1,2,…,N_T}; voxel-pairs that are not part of a significantly plastic sub-regional-pair are assigned a value of 0. This process is repeated for the N_U significantly plastic sub-regional-pairs identified in the second set of labelings. The ARI compares these labelings in such a way as to be invariant to label permutations. If the ARI is 1 then the overlap is perfect (i.e., the labelings are identical), and if the ARI is 0 then any observed overlap is due to chance. Thus the ARI metric examines how consistently the voxel-pairs are grouped into sub-regional-pairs across independent runs.

Per subject per ROI-pair, the ARI metric was calculated for 10 pair of independent runs. In Table 2B, we have summarized the mean and the standard deviation of the ARI score for the TBI and the HC group for all four ROI-pairs. The ARI score varied between 0.37 and 0.52. Also, based on Wilcoxon rank sum test we found no significant difference (p > 0.05) between the ARI scores of the TBI and HC for 3 of the 4 ROI-pairs; significant difference (p = 0.04) was only observed for the ROI-pair DMN-frontal-right and ECN-frontal-left. However, based on t-test no significant (p > 0.05) difference was found for all the 4 ROI-pairs.

TABLE 2B

Table 2B. Sub-cluster-pair overlap analysis based on Adjusted rand index.

Voxel-Pair Consistency Analysis

In this section we examine how consistently a voxel-pair is assigned to a significantly plastic sub-regional-pair across runs. In an ideal scenario, the EC approach will always (across multiple runs) identify a particular voxel-pair as either (a) existing in a significantly plastic sub-regional-pair or (b) not existing in a significantly plastic sub-regional-pair. To estimate the voxel-pair consistency we performed 5 independent runs per subject per ROI-pair. The results are summarized in Table 3; we found the voxel-pair consistency to be around 85% across all runs.

TABLE 3

Table 3. Link percentage consistency analysis.

Correlation Analysis of Plasticity Estimates

If we add the number of connections contributing to positive plasticity and the number of voxel-wise connections contributing to negative plasticity, then we obtain overall plasticity; i.e., edges gained/lost (gl₁).

In Section Simple Voxel-Level Approach we conducted a simple voxel-wise analysis by counting the number of voxel-wise functional connections that existed between a pair of ROI at session-1 and session-2. Thus, using the simple voxel-wise approach of section-Simple Voxel-Level Approach, we can derive another estimate for connections gained/lost (gl₂). Even though the plasticity estimate gl₂ is quite unrefined and may include functional connections that are not part of any statistically significant (see section Objective Function) sub-regional-pairs, it is reasonable to expect that the estimates gl₁and gl₂ should at least indicate similar plasticity effect in the brain. In order to examine this, for each subject for each ROI-pair we established the values: gl₁ and gl₂. Then we performed a simple correlation analysis. As we expected the correlation value was very high (between 0.97 and 1; see Figure S1). This helped us to further validate that the results obtained using evolutionary computation approach are reasonable.

Discussion

In this paper we introduced a novel voxel-level procedure based on evolutionary computation (EC) to examine plasticity between an investigator-defined ROI-pair. In addition to allowing the researcher to define the ROI-pair a priori, the procedure is to our knowledge unique in that it automatically defines the correct spatial scales at which the plasticity is occurring. The procedure reliably detects statistically significant positively and negatively plastic sub-regional-pairs by simultaneously analyzing subject-specific BOLD-fMRI data collected from two fMRI-scan sessions that are separated by a finite duration of time.

Compared to averaging based method of examining functional connectivity, we found that the EC-based procedure showed improved sensitivity to subtle changes in connectivity occurring over time which otherwise would be hard to detect (see Figure 7 and Tables S1–S4). The subject specific plasticity estimates (i.e., % voxel-wise connections contributing to positive and negative plasticity) obtained using the EC approach across multiple independent runs were highly consistent (see section Examining Functional Plasticity); the standard deviation of negative and positive plasticity estimates were less than 0.5% for almost all cases (and Tables S1–S4).

Even though the EC approach was developed to estimate plasticity, in the future these sub-regional-pairs may serve as functionally discrete regions for further analysis (i.e., plasticity-specific brain parcellation for network analysis). For such applications, it would be important that the sub-regional-pairs overlap significantly across independent runs. We used two measures to quantify the reliability of the EC approach in this secondary sense. First, we documented the spatial overlap between the sub-regional-pairs found across independent runs using the Sorenson-dice index (a score of 0 means no overlap and 1 means complete overlap). Second, we documented the extent to which the mappings of voxel-pairs to sub-regional-pairs are consistent across independent runs according to the ARI score (a score of 0 means minimal overlap due to chance and 1 means complete overlap). On average, the Sorenson overlap score was 0.8 and the ARI score was between 0.37 and 0.52, which are comparable to the spatial overlap across multiple runs obtained using other voxel-based approaches for brain parcellation (Thirion et al., 2014; Wang et al., 2015).

Since the overlap score is quite high (Sorenson score = 0.8), the solutions are very similar in terms of the spatial location they encode, and therefore, the sub-regional-pairs detected may be reliably used for further post hoc analysis. As one more step to document the overlap between the solutions found by EC across multiple runs, we analyzed the consistency by which voxel-pairs are detected as a part of a significantly plastic sub-regional-pair. We found that the voxel-pair consistency to be highly reliable (85%). This consistency is particularly impressive given that most voxelwise clustering approaches are designed primarily for data parcellation (clustering of homogeneous voxels) and do not consider pairs of clusters and do not utilize multiple fMRI-scan session data simultaneously, which are added dimensions of complexity in the current approach. Overall, the EC approach showed good consistency in terms of plasticity estimates and also spatial overlap across independent runs. It was a primary goal of this work to validate this EC-approach, so we did not include additional constraints on the algorithm, but we anticipate the spatial-overlap for identification of sub-regional pairs can be further improved by imposing physiologically or anatomically relevant constraints during the evolutionary search.

Using this EC procedure we can quantify plasticity on the population level, including the typical levels of positive and negative plasticity in the HC brain (see Figure 7). This information will provide a reference point for examining the positive/negative plasticity during any developmental process (e.g., learning, aging, recovery). For example, based upon this HC baseline variability, when examining the TBI cases, we see that 11/14 TBI subjects show heightened positive plasticity, and about 2/14 TBI subjects (subject-1 and subject-5) show heightened negative plasticity between DMN-frontal right and ECN-frontal right. These findings could be very important. A sub-regional-pair showing significant positive functional plasticity indicates that across time the two sub-regions have become more synchronized, and therefore, there is a possibility that the amount of metabolic energy consumed by the two regions to maintain synchronization may have increased as well (Nugent et al., 2015). In context of the current work, the positive plasticity could imply a heightened metabolism in the brain of a TBI associated with elevated synchronization between DMN and ECN. By contrast, negative plasticity after TBI has been linked to disruption in specific white matter tracts connecting the regions (Bonnelle et al., 2011; Fagerholm et al., 2015). The results found using the EC based approach is quite consistent with what has been demonstrated in the literature regarding the plasticity between DMN and ECN during recovery from TBI. The hyperconnectivity reported here has been observed in this and other datasets of individuals with TBI revealing increased synchrony in large-scale networks (Nakamura et al., 2009; Sharp et al., 2011; Hillary et al., 2011a,b; Bonnelle et al., 2012; Caeyenberghs et al., 2012).

There are several advantages afforded by the EC-approach used here. First, it offers previously unavailable sensitivity to the sub-regional response within an investigator defined ROI-pair. That is, the current approach permits the investigator to select theoretically relevant ROIs for hypothesis testing while maintaining sensitivity to the change in the functional relationship between the ROIs occurring at a much finer scale than the ROIs selected. As a secondary characteristic of this approach, the investigator is afforded the unique opportunity to examine functionally discrete subdivisions within the ROI-pair, which may be helpful to examine subtle functional responses between brain regions. For example, as shown in Figure 8, between DMN-frontal left and ECN-frontal right, we found several significantly plastic sub-regional-pairs which would have remained undetected by conventional averaging based approach. In the context of TBI literature, this approach may be very useful in examining the relationship of several coarsely defined hub regions to other sub-networks in the brain. For example, hub regions such as the PCC thought to be central to the DMN functioning is now also recognized to hold differential responses to goal-directed and internal-state networks depending upon where within the PCC the investigator is sampling data (Leech et al., 2011). Based on the results discussed in this article, we anticipate the EC based procedure to be well suited for such research problems where the investigator defined ROIs are coarsely defined and may contain within ROI heterogeneity.

Finally, a very important feature of the EC approach is that it uses two fMRI-scan session data simultaneously, thereby, allowing the sub-regional-pairs to be held constant across time. Simultaneously using multiple fMRI-scan session data is an important issue that is under-represented in the imaging literature. For any investigation of brain plasticity, if the regions are not held constant across time, then interpreting the results can become very challenging. The EC based procedure thus contributes to this important research problem.

It should be noted that, consistent with traditional pre-processing steps, we used a 6 mm smoothing kernel (2 voxels wide). While some degree of spatial smoothing is standard for maximizing signal-to-noise-ratio in fMRI datasets and helps to minimize ringing artifacts (Lindquist and Wager, 2008), spatial blurring can introduce spurious local connectivity. However, because we use an identical kernel for all analyses (e.g., averaging vs. EC comparison) and the physical distances between the ROIs chosen in these analyses are greater than the smoothing kernel, we do not anticipate a significant effect of smoothing on the primary findings. For future applications, if the EC approach is used to examine within-ROI plasticity (instead of between-ROI plasticity) then a large smoothing kernel could influence the local connectivity results. EC is a flexible paradigm and we hope to integrate adaptive smoothing (Lindquist and Wager, 2008) to this procedure for effective within-ROI plasticity analysis in the future.

In summary, we present a sensitive and reliable EC approach which has immediate importance for multiple areas of work in the cognitive and clinical neurosciences. We demonstrate sensitivity to the effects that would remain undetected using conventional approaches to examining between-ROI connectivity changes across time. The procedure will be very useful for any longitudinal design and, given recent emphasis on “precision medicine,” ideal for examining subtle connectivity changes over time in an individual while maintaining superior granularity with respect to the spatial dimensions of this change. The approach offers sensitivity to subtle between-ROI changes occurring during developmental, interventional, or clinical recovery processes, which may be very important for developing bio-markers for clinical applications. In the future, we hope to extend this work by integrating anatomically relevant information to constrain the evolutionary search and extend the EC framework to examine within-ROI plasticity.

Author Contributions

AR, CC, FH developed the methodology. AR, CC, FH designed the experiments and analysis. CC provided conceptual advice on computational aspect of this work. FH and RB provided conceptual advice on clinical utility of this work. AR developed the software programs and conducted the computational experiments. FH and RB collected the Traumatic brain injury data. AR and RB preprocessed the fMRI data. AR, CC, FH, and RB wrote the manuscript.

Funding

This project was supported by the Social Sciences Research Institute, in University Park, PA and the Penn State CTSI Grant (UL Tr000127) from the National Center for Advancing Translational Sciences, National Institutes of Health. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Dr. Reka Albert, Department of Physics, Pennsylvania State University, and her research group, Dr. Debashis Ghosh, Department of Biostatistics and Informatics, University of Colorado, Anschutz Medical Campus, and Dr. Peter Molenaar, Department of Human Development and Family Studies, Pennsylvania State University, for their valuable comments that allowed us to improve the article.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnins.2016.00146

References

Arenivas, A., Diaz-Arrastia, R., Spence, J., Cullum, M. C., Krishnan, K., Bosworth, C., et al. (2014). Three approaches to investigating functional compromise to the default mode network after traumatic axonal injury. Brain Imaging Behav. 8, 407–419. doi: 10.1007/s11682-012-9191-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassett, D. S., Wymbs, N. F., Porter, M. A., Mucha, P. J., Carlson, J. M., and Grafton, S. T. (2011). Dynamic reconfiguration of human brain networks during learning. Proc. Natl. Acad. Sci. U.S.A. 108, 7641–7646. doi: 10.1073/pnas.1018985108

PubMed Abstract | CrossRef Full Text | Google Scholar

Biswal, B., Yetkin, F. Z., Haughton, V. M., and Hyde, J. S. (1995). Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn. Reson. Med. 34, 537–541. doi: 10.1002/mrm.1910340409

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonnelle, V., Ham, T. E., Leech, R., Kinnunen, K. M., Mehta, M. A., Greenwood, R. J., et al. (2012). Salience network integrity predicts default mode network function after traumatic brain injury. Proc. Natl. Acad. Sci. U.S.A. 109, 4690–4695. doi: 10.1073/pnas.1113455109

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonnelle, V., Leech, R., Kinnunen, K. M., Ham, T. E., Beckmann, C. F., De Boissezon, X., et al. (2011). Default mode network connectivity predicts sustained attention deficits after traumatic brain injury. J. Neurosci. 31, 13442–13451. doi: 10.1523/jneurosci.1163-11.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Caeyenberghs, K., Leemans, A., Heitger, M. H., Leunissen, I., Dhollander, T., Sunaert, S., et al. (2012). Graph analysis of functional brain networks for cognitive control of action in traumatic brain injury. Brain 135(Pt. 4), 1293–1307. doi: 10.1093/brain/aws048

PubMed Abstract | CrossRef Full Text | Google Scholar

Davoian, K., Reichel, A., and Lippe, W. M. (2006). Comparison and analysis of mutation-based evolutionary algorithms for ANN parameters optimization. Training 6, 13.

Google Scholar

De Jong, K. A. (2006). Evolutionary Computation: A Unified Approach. Optimization. Available online at: http://mitpress.mit.edu/0262041944.

Derrac, J., García, S., Molina, D., and Herrera, F. (2011). A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms. Swarm Evol. Comput. 1, 3–18. doi: 10.1016/j.swevo.2011.02.002

CrossRef Full Text | Google Scholar

Diamond, M. C., Dowling, G. A., and Johnson, R. E. (1981). Morphologic cerebral cortical asymmetry in male and female rats. Exp. Neurol. 71, 261–268. doi: 10.1016/0014-4886(81)90087-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Eshelman, L. J. (1991). “The {CHC} adaptive search algorithm: how to safe search when engaging in nontraditional genetic recombination,” in Foundations of Genetic Algorithms, ed G. J. E. Rawlins (San Mateo, CA: Morgan Kaufmann Publishers), 265–283.

Eshelman, L. J., and Schaffer, J. D. (1991). “Spurious correlations and premature convergence in genetic algorithms,” in Foundations of Genetic Algorithms, ed G. J. E. Rawlins (San Mateo, CA: Morgan Kaufmann Publishers), 102–112.

Google Scholar

Fagerholm, E. D., Hellyer, P. J., Scott, G., Leech, R., and Sharp, D. J. (2015). Disconnection of network hubs and cognitive impairment after traumatic brain injury. Brain 138, 1696–1709. doi: 10.1093/brain/awv075

PubMed Abstract | CrossRef Full Text | Google Scholar

Fogel, D. B. (2003). “Evolving neural networks: selected medical applications and the effects of variation operators,” in Modeling and Simulation: Theory and Practice, eds B. A. Bekey and B. Y. Kogan (La Jolla, CA: Springer), 217–248.

Fogel, L. J. (1999). Intelligence through Simulated Evolution: Forty Years of Evolutionary Programming. New York, NY: John Wiley.

Google Scholar

Fox, M. D., Snyder, A. Z., Vincent, J. L., Corbetta, M., Van Essen, D. C., and Raichle, M. E. (2005). The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proc. Natl. Acad. Sci. U.S.A. 102, 9673–9678. doi: 10.1073/pnas.0504136102

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J., Frith, C. D., Liddle, P. F., and Frackowiak, R. S. (1993). Functional connectivity: the principal-component analysis of large (PET) data sets. J. Cereb. Blood Flow Metab. 13, 5–14. doi: 10.1038/jcbfm.1993.4

PubMed Abstract | CrossRef Full Text | Google Scholar

Greenough, W. T., and Chang, F-. L. F. (1988). “Plasticity of synapse structure and pattern in the cerebral cortex,” in Cerebral Cortex, eds A. Peters and E. G. Jones (New York, NY: Plenum Press), 391–440. doi: 10.1007/978-1-4615-6619-9_11

CrossRef Full Text | Google Scholar

Greicius, M. D., Supekar, K., Menon, V., and Dougherty, R. F. (2009). Resting-state functional connectivity reflects structural connectivity in the default mode network. Cereb. Cortex 19, 72–78. doi: 10.1093/cercor/bhn059

PubMed Abstract | CrossRef Full Text | Google Scholar

Gruau, F. (1994). Automatic definition of modular neural networks. Adapt. Behav. 3, 151–183. doi: 10.1177/105971239400300202

CrossRef Full Text | Google Scholar

Gruau, F., and Whitley, D. (1993). Adding learning to the cellular development of neural networks: evolution and the baldwin effect. Evol. Comput. 1, 213–233. doi: 10.1162/evco.1993.1.3.213

CrossRef Full Text | Google Scholar

Hebb, D. O. (1947). The effects of early experience on problem solving at maturity. Am. Psychol. 2, 306–307.

Heller, R., Stanley, D., Yekutieli, D., Rubin, N., and Benjamini, Y. (2006). Cluster-based analysis of FMRI data. Neuroimage 33, 599–608. doi: 10.1016/j.neuroimage.2006.04.233

PubMed Abstract | CrossRef Full Text | Google Scholar

Hillary, F. G., Medaglia, J. D., Gates, K., Molenaar, P. C., Slocomb, J., Peechatka, A., et al. (2011b). Examining working memory task acquisition in a disrupted neural network. Brain 134, 1555–1570. doi: 10.1093/brain/awr043

PubMed Abstract | CrossRef Full Text | Google Scholar

Hillary, F. G., Rajtmajer, S. M., Roman, C. A., Medaglia, J. D., Slocomb-Dluzen, J. E., Calhoun, V. D., et al. (2014). The rich get richer: brain injury elicits hyperconnectivity in core subnetworks. PLoS ONE 9:e104021. doi: 10.1371/journal.pone.0104021

PubMed Abstract | CrossRef Full Text | Google Scholar

Hillary, F. G., Roman, C. A., Venkatesan, U., Rajtmajer, S. M., Bajo, R., and Castellanos, N. D. (2015). Hyperconnectivity is a fundamental response to neurological disruption. Neuropsychology 29, 59–75. doi: 10.1037/neu0000110

PubMed Abstract | CrossRef Full Text | Google Scholar

Hillary, F. G., Slocomb, J., Hills, E. C., Fitzpatrick, N. M., Medaglia, J. D., Wang, J., et al. (2011a). Changes in resting connectivity during recovery from severe traumatic brain injury. Int. J. Psychophysiol. 82, 115–123. doi: 10.1016/j.ijpsycho.2011.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobayashi, K., and Ohbayashi, M. (1999). “A new indirect encoding method with variable length gene code to optimize neural network structures,” in IJCNN'99. International Joint Conference on Neural Networks. Proceedings (Cat. No.99CH36339) (Washington, DC), 4409–4412.

Leech, R., Kamourieh, S., Beckmann, C. F., and Sharp, D. J. (2011). Fractionating the default mode network: distinct contributions of the ventral and dorsal posterior cingulate cortex to cognitive control. J. Neurosci. 31, 3217–3224. doi: 10.1523/JNEUROSCI.5626-10.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindquist, M. A., and Wager, T. D. (2008). Spatial smoothing in fMRI using prolate spheroidal wave functions. Hum. Brain Mapp. 29, 1276–1287. doi: 10.1002/hbm.20475

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, C., Guo, Z. W., Lai, Y. X., Liao, W., Liu, Q., Kendrick, K. M., et al. (2012). Musical training induces functional plasticity in perceptual and motor networks: insights from resting-state fMRI. PLoS ONE 7:e36568. doi: 10.1371/journal.pone.0036568

PubMed Abstract | CrossRef Full Text | Google Scholar

Mazaika, P. K., Hoeft, F., Glover, G. H., and Reiss, A. L. (2009). Methods and software for fMRI analysis for clinical subjects. Neuroimage 47:S58. doi: 10.1016/s1053-8119(09)70238-1

CrossRef Full Text

Murphy, K., Birn, R. M., Handwerker, D. A., Jones, T. B., and Bandettini, P. A. (2009). The impact of global signal regression on resting state correlations: are anti-correlated networks introduced? Neuroimage 44, 893–905. doi: 10.1016/j.neuroimage.2008.09.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakamura, T., Hillary, F. G., and Biswal, B. B. (2009). Resting network plasticity following brain injury. PLoS ONE 4:e8220. doi: 10.1371/journal.pone.0008220

PubMed Abstract | CrossRef Full Text | Google Scholar

Nugent, A. C., Martinez, A., D'Alfonso, A., Zarate, C. A., and Theodore, W. H. (2015). The relationship between glucose metabolism, resting-state fMRI BOLD signal, and GABAA-binding potential: a preliminary study in healthy subjects and those with temporal lobe epilepsy 7. J. Cereb. Blood Flow Metab. 35 (1559–7016 (Electronic)), 583–591. doi: 10.1038/jcbfm.2014.228

PubMed Abstract | CrossRef Full Text | Google Scholar

Power, J. D., Barnes, K. A., Snyder, A. Z., Schlaggar, B. L., and Petersen, S. E. (2012). Spurious but systematic correlations in functional connectivity mri networks arise from subject motion. Neuroimage 59, 2142–2154. doi: 10.1016/j.neuroimage.2011.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Power, J. D., Mitra, A., Laumann, T. O., Snyder, A. Z., Schlaggar, B. L., and Petersen, S. E. (2014). Methods to detect, characterize, and remove motion artifact in resting state fMRI. Neuroimage 84, 320–341. doi: 10.1016/j.neuroimage.2013.08.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Raichle, M. E., MacLeod, A. M., Snyder, A. Z., Powers, W. J., Gusnard, D. A., and Shulman, G. L. (2001). A default mode of brain function. Proc. Natl. Acad. Sci. U.S.A. 98, 676–682. doi: 10.1073/pnas.98.2.676

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajtmajer, S. M., Roy, A., Albert, R., Molenaar, P. C., and Hillary, F. G. (2015). A voxelwise approach to determine consensus regions-of-interest for the study of brain network plasticity. Front. Neuroanat. 9:97. doi: 10.3389/fnana.2015.00097

PubMed Abstract | CrossRef Full Text

Rechenberg, I. (1971). Optimierung Technischer Systeme nach Prinzipien der Biologischen Evolution. Dissertation of Dr. Rechenberg.

Roy, A., David Schaffer, J., and Laramee, C. B. (2013). Evolving spike neural network sensors to characterize the alcoholic brain using visually evoked response potential. Procedia Comput. Sci. 20, 27–32. doi: 10.1016/j.procs.2013.09.234

CrossRef Full Text | Google Scholar

Saad, Z. S., Gotts, S. J., Murphy, K., Chen, G., Jo, H. J., Martin, A., et al. (2012). Trouble at rest: how correlation patterns and group differences become distorted after global signal regression. Brain Connect. 2, 25–32. doi: 10.1089/brain.2012.0080

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwefel, H. P. (1974). Adaptive Mechanismen in der biologischen Evolution und ihr Einfluss auf die Evolutionsgeschwindigkeit. Interner Bericht der Arbeitsgruppe Bionik und Evolutionstechnik am Institut für Mess-und Regelungstechnik Re.

Sharp, D. J., Beckmann, C. F., Greenwood, R., Kinnunen, K. M., Bonnelle, V., De Boissezon, X., et al. (2011). Default mode network functional and structural connectivity after traumatic brain injury. Brain 134(Pt. 8), 2233–2247. doi: 10.1093/brain/awr175

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheline, Y. I., and Raichle, M. E. (2013). Resting state functional connectivity in preclinical Alzheimer disease. Biol. Psychiatry 74, 340–347. doi: 10.1016/j.biopsych.2012.11.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Sorenson, T. J. (1948). A method of establishing groups of equal amplitude in plant sociology based on similarity of species content. Kongelige Danske Videnskabernes Selskab 5, 4–7.

Teasdale, G., and Jennett, B. (1974). Assessment of coma and impaired consciousness. A practical scale. Lancet 2, 81–84. doi: 10.1016/S0140-6736(74)91639-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Thirion, B., Varoquaux, G., Dohmatob, E., and Poline, J. B. (2014). Which fMRI clustering gives good brain parcellations?. Front. Neurosci. 8:167. doi: 10.3389/fnins.2014.00167

PubMed Abstract | CrossRef Full Text | Google Scholar

Tijms, B. M., Wink, A. M., de Haan, W., van der Flier, W. M., Stam, C. J., Scheltens, P., et al. (2013). Alzheimer Disease: connecting findings from graph theoretical studies of brain networks. Neurobiol. Aging 34, 2023–2036. doi: 10.1016/j.neurobiolaging.2013.02.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix, N., et al. (2002). Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. Neuroimage 15, 273–289. doi: 10.1006/nimg.2001.0978

PubMed Abstract | CrossRef Full Text | Google Scholar

van den Heuvel, M. P., and Hulshoff Pol, H. E. (2010). Exploring the brain network: a review on resting-state fMRI functional connectivity. Eur. Neuropsychopharmacol. 20, 519–534. doi: 10.1016/j.euroneuro.2010.03.008

PubMed Abstract | CrossRef Full Text

Venkatesan, U. M., Dennis, N. A., and Hillary, F. G. (2014). Chronology and chronicity of altered resting-state functional connectivity after traumatic brain injury. J. Neurotrauma 32, 252–264. doi: 10.1089/neu.2013.3318

PubMed Abstract | CrossRef Full Text | Google Scholar

Vinh, N. X., Epps, J., and Bailey, J. (2009). “Information theoretic measures for clusterings comparison: is a correction for chance necessary,” in Proceedings of the 26th Annual International Conference on Machine Learning (New York, NY: ACM), 1073–1080.

Google Scholar

Voss, M. W., Erickson, K. I., Prakash, R. S., Chaddock, L., Kim, J. S., Alves, H., et al. (2013). Neurobiological markers of exercise-related brain plasticity in older adults. Brain Behav. Immun. 28, 90–99. doi: 10.1016/j.bbi.2012.10.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z. I., Jones, S. E., Jaisani, Z., Najm, I. M., Prayson, R. A., Burgess, R. C., et al. (2015). Voxel-based morphometric magnetic resonance imaging (mri) postprocessing in mri-negative epilepsies. Ann. Neurol. 77, 1060–1075. doi: 10.1002/ana.24407

PubMed Abstract | CrossRef Full Text | Google Scholar

Whitfield-Gabrieli, S., and Nieto-Castanon, A. (2012). A functional connectivity toolbox for correlated and anticorrelated brain networks. Brain Connectivity 2, 125–141. doi: 10.1089/brain.2012.0073

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, X., Kelley, S., Goldberg, M., and Biswal, B. B. (2011). Detecting overlapped functional clusters in resting state fMRI with connected iterative scan: a graph theory based clustering algorithm. J. Neurosci. Methods 199, 108–118. doi: 10.1016/j.jneumeth.2011.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Yeo, B. T., Krienen, F. M., Chee, M. W., and Buckner, R. L. (2014). Estimates of segregation and overlap of functional connectivity networks in the human cerebral cortex. Neuroimage 88, 212–227. doi: 10.1016/j.neuroimage.2013.10.046

CrossRef Full Text | Google Scholar

Zalesky, A., Fornito, A., Egan, G. F., Pantelis, C., and Bullmore, E. T. (2012). The relationship between regional and inter-regional functional connectivity deficits in schizophrenia. Hum. Brain Mapp. 33, 2535–2549. doi: 10.1002/hbm.21379

PubMed Abstract | CrossRef Full Text | Google Scholar