Abstract
Current theory proposes that healthy neural dynamics operate in a metastable regime, where brain regions interact to simultaneously maximize integration and segregation. Metastability may confer important behavioral properties, such as cognitive flexibility. It is increasingly recognized that neural dynamics are constrained by the underlying structural connections between brain regions. An important challenge is, therefore, to relate structural connectivity, neural dynamics, and behavior. Traumatic brain injury (TBI) is a pre-eminent structural disconnection disorder whereby traumatic axonal injury damages large-scale connectivity, producing characteristic cognitive impairments, including slowed information processing speed and reduced cognitive flexibility, that may be a result of disrupted metastable dynamics. Therefore, TBI provides an experimental and theoretical model to examine how metastable dynamics relate to structural connectivity and cognition. Here, we use complementary empirical and computational approaches to investigate how metastability arises from the healthy structural connectome and relates to cognitive performance. We found reduced metastability in large-scale neural dynamics after TBI, measured with resting-state functional MRI. This reduction in metastability was associated with damage to the connectome, measured using diffusion MRI. Furthermore, decreased metastability was associated with reduced cognitive flexibility and information processing. A computational model, defined by empirically derived connectivity data, demonstrates how behaviorally relevant changes in neural dynamics result from structural disconnection. Our findings suggest how metastable dynamics are important for normal brain function and contingent on the structure of the human connectome.
Introduction
To understand how cognitive–emotional functions emerge from the brain and are affected by disease requires an account of how neural ensembles act in concert to generate behavior (Uhlhaas and Singer, 2006; Deco et al., 2008; Chialvo, 2010; Tognoli and Kelso, 2014). One approach is to consider the brain as a complex system (Friston, 1997; Chialvo, 2010). In this framework, current theory suggests an essential property of neural activity is metastability, a dynamical regime in which neural ensembles are able to coordinate rapidly, flexibly engaging and disengaging without becoming locked into fixed interactions (Friston, 1997; Shanahan, 2010; Tognoli and Kelso, 2014). Metastability is thought to confer optimal information processing capabilities, flexible behavior, and memory (Werner, 2007; Deco et al., 2009a; Shanahan, 2010). However, there is little empirical evidence to support this hypothesis. We have previously shown, with both empirical and computational approaches, that metastability at rest is higher than during a focused cognitive task (Hellyer et al., 2014). During a task, high metastability may be undesirable, where instead a specific configuration of brain systems is maintained over time (e.g., in coordinating specific visual and motor systems to perform a visually cued motor task).
Theoretical studies demonstrate that the emergence of metastable dynamics is contingent on the coupling between modules of a dynamical system (Friston, 1997; Strogatz, 2001; Shanahan, 2010; Cabral et al., 2011). In particular, dynamic patterns of functional connectivity, consistent with metastable dynamics, emerge when coupling has “small-world” topology with short average path lengths and high clustering (Wildie and Shanahan, 2012) of modules. Recently, networks of anatomical connections have been incorporated within computational simulations of large-scale neural dynamics, suggesting metastable dynamics provide a link between structural and functional connectivity (Deco et al., 2009b; Honey et al., 2009; Cabral et al., 2011; Hellyer et al., 2014).
The disruption of neural dynamics is thought to be important in brain disorders (Uhlhaas and Singer, 2006), likely caused by abnormal structural connectivity (Friston, 2002; Bassett and Bullmore, 2006; Cabral et al., 2012; Sharp et al., 2014). Indeed, through the examination of brain disorders, we are able to explore the importance of structural connectivity for the organization of functional connectivity. A large body of work has examined the link between focal damage to the brain, such as in ischemic stroke, and cognition (Bird et al., 2004; Sharp et al., 2010b; Gratton et al., 2012; Warren et al., 2014). However, structural disconnection is often intermingled with gray matter damage (Sharp et al., 2011; Bonnelle et al., 2012). In contrast, traumatic brain injury (TBI) frequently results in diffuse axonal injury (DAI), which disrupts long-distance white matter tracts connecting brain regions (Sharp et al., 2011; Johnson et al., 2013a,b), but with neuronal bodies relatively spared. As such, it is a pre-eminent example of a white matter disconnection disorder (Sharp et al., 2014). Damage to white matter connectivity in TBI alters the spatiotemporal properties of functional brain networks (Kinnunen et al., 2011; Hellyer et al., 2013; Caeyenberghs et al., 2014; Jilka et al., 2014), resulting in long-term cognitive problems, including impairments in cognitive flexibility, memory, and information processing speed (Bonnelle et al., 2011; Kinnunen et al., 2011; Jilka et al., 2014). Cognitive inflexibility after TBI may be observed as poor performance on tests of task switching (Kinnunen et al., 2011; Hellyer et al., 2013; Caeyenberghs et al., 2014; Jilka et al., 2014). Extreme inflexibility may manifest as perseveration, the repetition of a particular response, such as a phrase or gesture, despite the cessation of a stimulus. Therefore, TBI provides an ideal paradigm to examine how three important levels of description (structure, functional dynamics, and behavior) converge.
Here, using empirical and computational approaches, we investigate how metastability, defined as the SD of the Kuramoto order parameter (Shanahan, 2010; Cabral et al., 2011), arises from the structural connectome and relates to behavior. We test whether (1) structural disconnection after TBI (measured using diffusion tensor imaging) is associated with reduced metastability [measured using resting-state functional MRI (fMRI)] and (2) metastability is associated with behavioral measures of cognitive flexibility, memory, and information processing. Furthermore, we use computational simulations to investigate the consequences of structural disconnection on large-scale neural dynamics, to demonstrate how disconnection after TBI results in altered metastability.
Materials and Methods
Overview.
Our experimental approach is split into three different sections, described in detail below. An overview of the methods is provided in Figures 1 and 2. In brief, first we used resting-state fMRI to estimate empirical measures of metastability in healthy control subjects and in the presence of structural disconnection in TBI patients (Fig. 1A,C). Second, we used a diffusion tensor imaging (DTI) approach to estimate both voxelwise [fractional anisotropy (FA)] and region of interest (ROI)-level connectivity for each subject (Fig. 2). Finally, we used a range of both computational (using the Kuramoto model of phase oscillators) and empirical imaging-based approaches to describe the relationship between structural connectivity, metastability, and neuropsychological performance (Fig. 1C).
Overview of experimental design. A, fMRI was used to estimate global measures of network dynamics during rest in 63 patients and 26 controls. B, We used a computational model to simulate neural dynamics using dynamic systems framework constrained by structural connectivity. C, We used a Hilbert transformation of the fMRI data and the phase output of the computational model to compare the global dynamics of empirical data and the dynamics of a computational model constrained by white matter structural connectivity.
Estimation of structural disconnection in traumatic brain injury patients versus controls. A, To define a standardized connectivity matrix to perform computational simulation, probabilistic tractography was performed in 10 independent age-matched healthy control subjects, resulting in a binary connectivity graph and spatial estimates of probable tract location for each connected edge (see Materials and Methods). B, Measures of tract integrity (FA) were estimated in each of the 63 patients and 26 healthy controls by generating a “skeleton” for each subject using the preprocessing steps of TBSS. C, For each edge of the reference connectivity matrix, each subject's FA skeleton was projected though the spatial mask for each edge, resulting in a 164-region “white matter integrity” matrix for each subject. C, For each subject, the reference binary connectivity map was weighted according to relative reduction in tract integrity in each individual subject (see Materials and Methods), resulting in a individual weighted connectivity matrix for each of the 63 TBI patients and 26 healthy controls.
Image acquisition.
Standard protocols were used to acquire functional, structural, and diffusion tensor MRI data using a Phillips Intera 3.0 tesla MRI scanner, with an eight-array head coil, and sensitivity encoding (SENSE) with an under sampling factor of 2. For each participant, diffusion-weighted volumes with gradients applied in 64 noncollinear directions were collected. The following parameters were used: 73 contiguous slices; slice thickness, 2 mm; field of view, 224 mm; matrix, 128 × 128 (voxel size, 1.75 × 1.75 × 2 mm3); b = 1000; four images with no diffusion weighting (b = 0 s/mm2). Earplugs and padded headphones were used to protect participants' hearing during the scanning procedure. We additionally collected a standard high-resolution T1 image for segmentation and image coregistration. During the resting-state fMRI scan, subjects were asked to keep their eyes closed and to try not to fall asleep. Functional volumes were collected using a T2*-weighted gradient echoplanar imaging sequence with whole-brain coverage (repetition time/echo time, 2000/30 ms; 31 ascending slices with thickness of 3.25 mm; gap, 0.75 mm; voxel size, 2.5 × 2.5 × 5 mm; flip angle, 90°; field of view, 280 × 220 × 123 mm; matrix, 112 × 87). Quadratic shim gradients were used to correct for magnetic field inhomogeneities within the brain.
Participants.
Sixty-three TBI patients (16 female; mean age ± SD, 37.4 ± 12.37 years) and 26 healthy control subjects (12 female; mean age ± SD, 35.96 ± 17.61 years) were scanned using standard functional and structural MRI protocols (see below). The patients and controls overlapped with those studied by Fagerholm et al. (2015), where group demographics are reported in more detail. At the group level, patients and controls were matched for age (t(87) = −0.47, p = 0.64) and gender. We did not hypothesize that there would be any effects of handedness, and so we did not select subjects according to handedness. TBI patients were scanned in the chronic phase [5.48 ± 3.33 (months ± SD)] after injury. Injury severity of TBI patients was classified according to the Mayo system (Malec et al., 2007): moderate-severe, 55 patients; mild (probable), 5 patients; symptomatic (possible), 3 patients. Forty-nine patients had a clinically relevant episode of post-traumatic amnesia after TBI. Where the mechanism of injury was known, the mechanism was road traffic accident (21 patients), assault (17), fall/syncope (15), and sports injury/concussion (2). For eight patients, the mechanism of injury was unknown. All participants gave written consent, were checked for contraindications to MRI scanning, and had no history of significant neurological or psychiatric illness before TBI. The Hammersmith, Queen Charlotte's, and Chelsea research ethics committee awarded ethical approval for the study.
Neuropsychological assessment.
All but one patient (62; mean age ± SD, 37.53 ± 12.45 years) performed a paper-and-pencil, neuropsychological test battery. Our analysis focused on cognitive measures shown previously to be sensitive to impairments after TBI: associative memory, using the immediate recall and retention measure of the People Test from the Doors and People Test (Baddeley, 1986, 1992; Baddeley et al., 1994), and executive function, using the Trail Making Test alternating switch cost index (SCI; Reitan, 1958). In addition, a subset of 49 patients also completed the computerized Choice Reaction Time task that assesses speed of processing (Rabbitt, 1966; Logan et al., 1984). It was for technical reasons, e.g., relating to equipment error, that not all subjects completed all three neuropsychological tests. Neuropsychological assessment was performed immediately before the MRI scanning session by a trained experimenter; scoring for each test was performed according to the protocols provided by the original publisher of each test, and no further selection based on neuropsychological outcome was performed. We present all collected data for each of these groups (i.e., no statistical subsampling on the basis of performance on any individual test has been performed).
Analysis of functional imaging data.
Preprocessing of functional data was performed according to standard analysis approaches: briefly, this included realignment of echoplanar images to remove the coarse effects of motion between scans using FMRIB's motion correction tool MCFLIRT (Smith et al., 2004; Fig. 1A). T1 images for each subject were segmented into 164 cortical and subcortical regions using the Destreux Freesurfer atlas (Fischl et al., 2004). The segmented T1 images were registered to the motion-corrected data using boundary-based registration (Greve and Fischl, 2009). Subsequently, mean BOLD time series for each ROI were extracted from the resting-state scans. We bandpass filtered the data between 0.01 and 0.2 Hz to remove sources of non-neural noise and focus on slow modulations in BOLD that have previously been associated with intrinsic connectivity networks (ICNs; Niazy et al., 2011). Such bandpass filtering is an important step in transforming raw time series into phase space (see below). To account for variance related to head motion or non-neural physiological noise, during preprocessing we regressed out from the time course for each of the 164 ROIs, the six motion parameter time courses estimated by MCFLIRT (Smith et al., 2004), the motion parameters squared, as well as time series sampled from regions of white matter and CSF. To further reduce the possibility that effects are driven by head motion, we controlled for motion at the group level (i.e., across subjects) by including the estimate of mean framewise displacement in higher-level analyses as a regressor of no interest.
To facilitate comparison of measures of metastability and synchrony between the computational model and empirical data, we transformed the empirical data into a complex phase representation, so the same analysis can be applied to both the empirical and computational simulation data. The transformation of functional neuroimaging data into phase representation for analysis has been previously performed using a variety of different approaches such as wavelet analysis (Kitzbichler et al., 2009; Chang and Glover, 2010), as well as the computationally simpler Hilbert transform on bandpass filtered data (Glerean et al., 2012). For simplicity, we perform the latter on each of the bandpass filtered 164 ROI time series from the empirical data, resulting in 164 phase time series. Measures of network dynamics were either calculated on all regions simultaneously (global) or within specific predefined intrinsic connectivity networks (local).
Definition of intrinsic connectivity networks from functional imaging data.
We defined a set of ICNs from the resting-state fMRI data to allow us to optimize the computational model and to define networks from which to sample neural dynamics. To do this, we performed temporal concatenation independent component analysis (ICA) on each of the 164 ROI mean BOLD time series for an independent group of 10 healthy control subjects, using FSL MELODIC (Beckmann et al., 2005). The optimal decomposition estimated during fast ICA resulted in the identification of 15 independent networks. These ICNs were used to tune the computational model (see below). Resulting components were thresholded at z > 2.3, and surviving brain regions were included in that network. Each of these networks were then labeled by eye based on their resemblance to the canonical ICNs produced by Smith et al. (2009). This resulted in the identification of seven canonical ICNs that were used in subsequent analyses (visual, auditory, default mode, dorsal attention, salience, and both left and right frontoparietal control networks; Fig. 3). These networks were used to sample metastability for both the empirical and simulated data.
ICA of fMRI data in macroscopic ROIs. ICA was used to decompose ROI data from 164 cortical and subcortical regions of the brain in 10 healthy control subjects into 15 components. Here, seven of the resulting components that best resembled canonical “resting-state networks” have been reconstructed into MNI 152 space according to the 164-region model with a Gaussian blur of 3 mm (isotropic) for ease of visualization. L, Left.
Estimation of healthy structural connectivity network.
The mean location and probability of structural connections was estimated in an additional group of 10 independent healthy control subjects (Fig. 2A). Structural T1 images were segmented into white matter and the same 164 cortical and subcortical gray matter ROIs as used to sample the fMRI data, using Freesurfer (Greve and Fischl, 2009). This produced a mask for each region in each participant's T1 native space. Diffusion imaging data were reconstructed using the FSL diffusion toolkit using standard protocols (Behrens et al., 2003b). We also modeled the probability distribution of fiber direction within each voxel to account for crossing fibers (Behrens et al., 2003a). Nonlinear registration was used to calculate a warp field between conformed Freesurfer space and the DTI b0 image, using the FSL nonlinear image registration tool (Smith et al., 2004). The warp field was then applied to masks for white matter and each of the 164 cortical and subcortical ROIs using nearest-neighbor interpolation. Individual gray matter masks were dilated by a single voxel and multiplied by the white matter mask to generate ROIs to be used as seeds and targets for tractography at the boundary between white and gray matter surfaces (Gong et al., 2009).
Probabilistic tractography, using 5000 random streamline samples per voxel, was used to estimate the connectivity matrix C between each of the 164 other regions alongside a spatial distribution of connective fibers between each region. The probability of connections between two regions C(i,j) was defined as the proportion of all fibers sent from region i, which successfully reached region j. Because probabilistic tractography cannot determine directionality of connections between cortical regions and the size of seed and target ROIs may differ for each connection, we define C(i,j) as the mean of the forward and reverse connections between regions, i.e., C(i,j) =
Estimation of individual structural connectivity.
Global and focal reductions of FA in TBI patients have been shown to bias tractography estimation in TBI patients (Squarcina et al., 2012), potentially resulting in spurious differences including false increases in structural connections after injury (Fig. 2B). This previous work suggested that sampling FA along a known tract distribution estimated from an independent group of healthy control subjects is preferable to estimating tracts from TBI patients themselves. Therefore, we sampled FA projected through a set of tracts linking the 164 regions, defined on the independent group of 10 healthy controls. First, we used tract-based spatial statistics (TBSS; Smith et al., 2006) to align the FA map of each patient and control subject to a common template. To reduce partial volume effects, these tracts were then skeletonized, resulting in a voxelwise map of the white matter skeleton for each subject (Fig. 4).
We used random permutation testing (Nichols and Holmes, 2002; Smith et al., 2004; Winkler et al., 2014) to assess whether there was reduced FA in patients compared with controls (Kinnunen et al., 2011; Hellyer et al., 2013). We fitted a general linear model for each voxel within the skeleton. Although this approach gives a good overview of spatial distribution of damage, we also wanted to estimate the integrity of individual predefined “tracts” for each patient and control to define a structural connectivity matrix. Therefore, mean FA values were calculated from masks generated at the intersection of voxels of the skeletonized map and each tract region of interest (as defined using tractography on independent controls; see the previous section). This resulted in a “tract integrity” FA matrix for each patient and control (Fig. 2C). To define “damage” to tracts, rather than use an arbitrary threshold of these matrices (which is likely to remove valid tracts in both the patient and control population with naturally low FA because of factors such as crossing fibers), we determined the tracts within which damage is likely by estimating the normal distribution (mean and SD) of FA values for each connection within the independent group of healthy controls (used to define the tractography). This information was then used to z-transform the FA values for each tract within each individual subject in the patient and control groups. To simulate damage to the connectivity matrix C, if any edge from the scaled FA matrix for an individual fell below a certain threshold (−1.6 SD, which represents the position within the normal Gaussian distribution where an individual event is distinct from noise with a nominal probability of approximately p < 0.05), it was “lesioned.” Rather than removing the tract (which is overly destructive given the nature of the traumatic axonal injury, where the tract typically remains but shows evidence of damage), we instead reduced connectivity by a fixed amount (Fig. 2C). This has the effect that each connectivity matrix has the same number of connections, while capturing any pattern of damage. Results reported are for a reduction of 50%, but the results were robust to a range of different damage values.
Graph theoretic metrics for structural connectivity.
To define changes in structural connectivity in relation to changes in neural dynamics, we assessed large-scale structural connectivity using measures from graph theory, calculated using the Brain Connectivity Toolbox (Rubinov and Sporns, 2010). These measures are described briefly below.
Degree (D) is defined as follows:
The degree of each node Di within a weighted graph C is defined as the sum of all directly connected edges to the node within the network. The mean degree defines how strongly interconnected all nodes within the network are (Freeman, 1978).
Characteristic path length (L̄) is defined as follows:
Path length (Li) is the average distance of an individual vertex to all its connected neighbors [v ∈ V(C)] in a network C, weighted by the inverse of the weight of connectivity (D), i.e., higher weight connections are interpreted as a shorter connection length. The mean of this value across all nodes (L̄) is the measure of characteristic path length within a network (Watts and Strogatz, 1998).
Clustering coefficient (K̄) is defined as follows:
The weighted clustering coefficient of a node (Ki) is the average connectivity strength of all “triangles,” i.e., all neighbors (ei), that also directly connect to each other as pairs (ki) around a specific node (i). The mean across nodes (K̄) is used as a measure of network clustering (Watts and Strogatz, 1998).
Small-worldness (σ) is defined as follows:
Small-world networks have a low-characteristic path length and a high-clustering coefficient. An often applied metric of small-worldness is the small-world index (σ − SWI); Sporns, 2006; Humphries and Gurney, 2008). This compares the path length (L) and clustering coefficient (K) to equivalent measures of a suitable Erdös-Rényi random network (Humphries and Gurney, 2008; K̂ and L̂, respectively). If σ > 1, a network is considered small-world.
Differences in σ and K̄ may be artifactually driven by a reduction in the mean degree of individual networks. To correct for this potential source of estimation error, we controlled for variation of mean degree across each subject's structural connectivity graph by normalizing the mean degree of each graph (i.e., dividing C by the mean of connected vertices in D) before calculating σ and K̄. In this way, changes in each of these measures will be driven by altered network topology, rather than simply global changes in connectivity strength.
Computational simulation of neural dynamics.
To explore the effect of structural disconnection on neural dynamics, we used the simple Kuramoto model of coupled phase oscillators (Acebrón et al., 2005; Fig. 1B). We chose this model, partly because it has been shown to be able to simulate macroscopic neural dynamics related to underlying structural connectivity (Shanahan, 2010; Cabral et al., 2011; Hellyer et al., 2014; Messé et al., 2014). It has also been shown to capture the same essential aspects of macroscopic dynamics as far more complex models based on tens of thousands of simulated neurons (Bhowmik and Shanahan, 2013). Compared with such complex and computationally far more intensive models, the Kuramoto model provides a good trade-off between complexity and plausibility, modeling a few key parameters of structural and functional relationships between nodes. By using the less computationally intensive Kuramoto model, we were able to explore large parameter spaces and simulate many individual subjects' dynamics. The activity of each of the 164 brain regions (which we define here as a node) is represented in our model by the phase of a Kuramoto oscillator. The phase of each node over time is described by the following Kuramoto equation (Kuramoto, 1984; Acebrón et al., 2005):
The natural frequency ω defines the phase change of an uncoupled node per time step. In our simulations, as in previous work (Cabral et al., 2011), we fixed the natural frequency to match known oscillations within the gamma frequency range (ω = 60 Hz). The connectivity matrix C is a binary connectivity matrix determined by the empirical strength of white matter connections or lesioned using individual tract integrity data (see above). The distance matrix D, determined by the empirical length of connections between regions, imposes time delay on phase interactions between nodes. This is analogous to the simulation of a delay caused by neural conduction between regions of the brain.
We introduce two control parameters to the coupling and delay of the network: the global coupling parameter (k) and the mean global velocity (v), such that C = kC and D = D/v. The behavior of the Kuramoto model in terms of global metastability and synchrony, by modulation of the parameters k and v, has been explored previously (Shanahan, 2010; Cabral et al., 2011). For completeness, we also present a traversal of this parameter space (see Fig. 8) using the baseline binary connectivity dataset defined using 10 healthy controls (see above, Estimation of healthy structural connectivity network). The presence of multiple local maxima of both metastability and synchrony within the k, v plane makes it challenging to optimize through a gradient-descent approach. Therefore, we randomly selected 6000 pairs of model parameters within the k, v plane and executed the model for each pair. We then used nearest-neighbor interpolation to create the parameter space k, v presented in Figure 8.
Validation of computational simulation.
To validate the computational simulation against empirical functional connectivity derived from fMRI BOLD data, we followed the approach to simulating BOLD activity previously demonstrated in the literature (Cabral et al., 2011), using the sine of the high-frequency activity of the Kuramoto model as the neural input to the Balloon-Windkessel hemodynamic model (Friston et al., 2000), low-pass filtered the resulting time courses at <0.25 Hz, and downsampled to a 2 s sampling rate.
Unlike previous approaches to model validation that have assumed that the spatiotemporal organization of correlations within empirical fMRI BOLD time courses are univariate, we used an ICA approach to compare empirical BOLD activity (defined using the ICA decomposition described above) with the output of our computational models. Simulated BOLD time courses from each of the computational models were decomposed into 15 spatially independent time courses. Spatial components from the empirical and modeled ICAs were then “matched” using spatial correlation of their maps. Since the order of components extracted by MELODIC varies, we determined maximal correspondence between empirical and model components by calculating the pairwise spatial correlation between functional connectivity maps for all pairs of components and reordering the resulting correlation matrix so as to maximize entries along the diagonal. An evaluation function was defined as the mean correlation between the empirical spatial maps and the modeled spatial maps for the top n (here n = 5) matched components in the computational model compared with the empirical data, providing an objective measure by which the fit of individual regions of the global parameter space to empirical data may be compared against one another (see Fig. 8). We present the results for ICA using 15 components in both empirical and modeled data, and n = 5; however, varying each of these parameters produced qualitatively similar results.
Synchrony and metastability measures of neural dynamics.
The order parameters R(t) and φ(t) can be jointly defined by the following:
where N is the total number of regions within the network and the level of synchrony between phase time courses is described by R(t), in terms of how coherently phase changes over time (Shanahan, 2010; Cabral et al., 2011; Fig. 1C). During fully synchronous behavior, R(t) = 1, whereas R(t) = 0 where phase across all phase time series is fully asynchronous. The phase of all the input phase time series is described by φ(t), but is not used in the present work. For both empirical and simulated time series, we measured neural dynamics in terms of mean global synchrony (R̄) and global metastability as the SD σR of global synchrony across the same period (Shanahan, 2010; Cabral et al., 2011). In addition to global measures of dynamics, to evaluate measures of network dynamics in both the empirical and simulated data, we calculated separate mean synchrony and metastability measures for the phase time series of regions within the seven empirical ICNs defined previously (see above, Definition of intrinsic connectivity networks from functional imaging data).
Results
Widespread disruption to the structural connectome after traumatic brain injury
To demonstrate that white matter connectivity is disrupted after TBI in the group of patients studied, we performed standard TBSS to compare white matter integrity, measured by FA, between patients and healthy controls using age and total gray matter volume as covariates of no interest. In the between-group contrast of patients < controls, there was widespread reduction in FA across the white matter skeleton (Fig. 4A). FA reduction was particularly pronounced in the interhemispheric fibers of the corpus callosum (where damage was widespread, but most extensive in the body and genu), as well as in tracts within the superior longitudinal fasiculus, corticospinal tract, and the anterior and posterior limbs of the internal capsule. Additionally, a strong reduction in FA was observed within the fornix and corona radiata. This distribution of widespread changes to white matter integrity is typical of injury after TBI and is consistent with our previous TBSS-based findings (Kinnunen et al., 2011).
Structural brain damage after traumatic brain injury. A, Widespread white matter disruption after traumatic brain injury measured by TBSS of FA. Contrasts between traumatic brain injury < healthy control subjects are shown in red-yellow. Contrasts are overlaid on a standard MNI 152 T1 1 mm brain and the mean FA skeleton (in green), thresholded at p ≤ 0.05, and corrected for multiple comparisons using Threshold Free Cluster Enhancement (Smith et al., 2006). B, Lesion probability maps of cortical contusions across 63 TBI patients, estimated by a neuroradiologist on the T1 structural images. The color bar indicates the number of patients who had lesions at each site, overlaid on a standard MNI 152 T1 2 mm brain. L, Left.
To explore how network-level measures of structural connectivity are changed after TBI, we tested for group-wise differences between graph theoretical measures of large-scale structural connectivity between healthy controls and TBI patients using age as a covariate of no interest. There was a significant reduction in patients compared with controls in small-worldness (t(86) = −3.10, p < 0.01), clustering coefficient (t(86) = −2.04, p < 0.05), and mean degree (t(86) = −3.42, p < 0.001). In addition, patients had a significantly higher characteristic path length compared with healthy controls (t(86) = 3.07, p < 0.01).
Empirical measures of metastability of large-scale neural dynamics are reduced after traumatic brain injury
We assessed the metastability of large-scale neural dynamics after TBI, which was measured using 164 regional phase time courses derived from resting-state fMRI BOLD data in both patients and controls (Fig. 1A). We compared metastability across the whole brain and within ICNs between patients and controls, including age, gray matter volume, and mean absolute movement as covariates of no interest. Global metastability was significantly reduced in TBI patients compared with controls (t(84) = −2.63; p < 0.05, one-tailed; Fig. 5). Patients also showed lower metastability in the salience network (t(84) = −3.68, p < 0.001), a left frontoparietal network (t(84) = −2.41, p < 0.02), and a dorsal attention network (t(84) = −2.27, p < 0.05); these survive multiple-comparison correction by false discovery rate (FDR; q < 0.1). A potential confound for measuring widespread neural dynamics in the TBI patient population is the presence of cortical contusions. Thirty-two patients in the TBI group were found to have focal gray matter lesions, suggestive of cortical contusions, on T1-weighted structural imaging. It is possible that these lesions affected the BOLD time courses extracted and the resulting metastability calculations. We therefore repeated the calculation of global metastability after removing brain regions whose anatomical segmentation overlapped with focal lesions in any patients (Fig. 4B; i.e., time courses from affected regions were not analyzed in any patients or controls). Reduction in global metastability in TBI patients compared with controls after this adjustment was similar to the original analysis (t(84) = −2.63, p < 0.01).
Empirical metastability at rest is significantly reduced in patients compared with controls. Mean measures of metastability (±1 SEM) estimated using a phase-transformed functional time course extracted from 63 patients and 26 control subjects suggest that global measures of metastability are reduced after TBI (*p < 0.05, **p < 0.01, ***p < 0.001, FDR corrected).
To establish whether the differences we report were attributable to changes in metastability (i.e., temporal variability in how synchronous the brain is) rather than simple temporal variability in the BOLD signal, we performed a follow-up analysis. We calculated the SD of (1) the mean BOLD signal across all 164 regions and (2) the mean BOLD signal of each region. In a GLM using age, gray matter volume, and mean absolute displacement as covariates of no interest, we found no group-wise effect in the SD of the mean global BOLD signal (t(84) = 0.22, p = 0.82). For each of the 164 regions, no single region showed a significant effect of group (FDR correcting for multiple comparisons; even without correcting for multiple comparisons, only four regions had a group-wise effect with p < 0.05, with the lowest value p = 0.04).
Empirical measures of metastability after traumatic brain injury predict cognitive performance
To investigate whether empirical measures of global metastability relate to the cognitive impairments seen in the TBI population, we regressed measures of metastability against measures of cognitive flexibility, associative memory, and information processing speed, including age and mean absolute movement during the fMRI run as covariates of no interest. In patients, there was a significant negative relationship between global metastability and SCI (where higher SCI suggests poorer flexibility; t(58) = −2.21, p < 0.05; Fig. 6, top) and median reaction time (t(44) = −3.46, p < 0.01; Fig. 6, middle). Global metastability was positively related to immediate memory recall, i.e., improved performance (t(58) = 2.49, p < 0.05; Fig. 6, bottom), and retention (t(58) = 2.780, p < 0.01; Fig. 6, bottom). These results were from a multiple regression model containing age, motion, and total gray matter volume as covariates of no interest. The results survive multiple-comparison correction with FDR (q < 0.1).
Global metastability predicts behavioral outcome. Measures of global metastability during rest significantly correlate with scores of task-switching (top, red, n = 62), information processing speed (middle, purple, n = 49), and logical memory (bottom, blue and green, n = 62). CRT, Choice reaction time.
Frontal disconnection after TBI predicts global measures of empirical metastability
We explored how regional structural disconnection may relate to metastability. We calculated the mean FA across the whole white matter skeleton and compared it with global metastability in patients and controls, with age as a covariate of no interest. We observed a significant main effect of group (t(84) = −2.12, p < 0.05) and interaction between group and FA (t(84) = 2.04, p < 0.05). We also explored this relationship at a regional level, performing a voxelwise regression of metastability on FA within both groups separately. In healthy control subjects, no region of the white matter skeleton was significantly associated with global metastability. Within the patient group, there was a significant association between global metastability and integrity of the white matter skeleton, predominantly within the white matter linking the frontal lobe with the thalamus through the anterior thalamic radiation (Fig. 7). However, when the two groups were compared directly, there were no voxels surviving multiple-comparison correction for a group by FA interaction, so these group differences should be treated with some caution.
Empirical metastability in patients is associated with cortical connectivity within the frontal lobes. Multiple regression of empirical global metastability (red-yellow) with skeletonized FA values in 63 TBI patients (age included as a covariate of no interest) is shown. Results are overlaid on a standard MNI 152 T1 1 mm brain and the mean fractional anisotropy skeleton (in green), thresholded at p ≤ 0.05, and corrected for multiple comparisons using Threshold Free Cluster Enhancement (Smith et al., 2006).
A computational simulation of macroscopic neural dynamics resembles empirically defined intrinsic connectivity networks
To explore how measures of metastable neural dynamics respond to structural disconnection (i.e., after TBI), we examined a computational simulation of the brain at rest using the Kuramoto model, constrained by the white matter structural connectivity between 164 brain regions, defined using white matter tractography. We validated the model by comparing it with resting-state functional connectivity derived from fMRI BOLD.
Previous work has explored the dynamics of the Kuramoto model in relation to the strength and structure of coupling between nodes (Shanahan, 2010; Cabral et al., 2011, 2012; Wildie and Shanahan, 2012). This work suggests that the model behavior is highly sensitive to two constants, the global coupling parameter (K) and mean global velocity (V), which is determined by a distance matrix, D. To understand the effects of these constants on our 164 node model, using coupling defined by the reference binary reference connectivity matrix (see Materials and Methods, Estimation of healthy structural connectivity network), we performed a parameter space search using 6000 randomly generated pairs of parameters k, v within the plane k(1…6) and V(1…16). The behavior of global metastability as well as mean global synchrony of the system as a function of K is shown in Fig. 8. We observed that for increasing values of K, the system tends toward maximum global synchrony after passing through an intermediate phase where metastability is maximal. To reduce the complexity of further computations, we selected a point in the V dimension based on plausible physiology, such that V = 11 ms−1, following Cabral et al. (2011).
A large-scale neural network model predicts functional connectivity where metastability is maximized. Plots of a parameter space exploration of the Kuramoto model using 6000 randomly generated pairs of the coupling and delay parameters k, v, within the plane k(1…6) and V(1…16), here expressed as mean velocity (ms−1), are shown. Top, Global synchrony; middle, global metastability; bottom, correlations of the five best-matching simulated and empirical ICNs. For ease of further computation, we set velocity at a biologically plausible value (Cabral et al., 2011), which equates to ∼11 ms−1 (dashed line).
For the model output for each k, v pair, we used ICA to decompose the 164 node time courses into a set of simulated ICNs. These simulated ICNs were then correlated with a set of resting-state networks derived from the empirical BOLD fMRI resting-state data in the same 10 independent healthy control subjects used in the tractography step (see Materials and Methods). We found that the correlation between the simulated and empirically defined networks was highest near the point of maximal metastability (Fig. 8). The emergence of functional networks around the region of maximal metastability suggests that such a rich dynamical regime is an important organizing principle of how structural connectivity may allow functional networks to form in the brain.
Empirically defined macroscopic structural disconnection leads to reduced metastability in a simulation of macroscopic neural dynamics
To examine the effect of macroscopic structural disconnection after TBI on simulated neural dynamics, we used individualized structural connectivity matrices in patients and controls (see Materials and Methods; Fig. 2C) to define coupling within the Kuramoto model (Fig. 1B). Global metastability was calculated for simulations of the model executed separately for each subject's connectivity matrix. Runs were repeated for a range of values of the coupling constant, K (Fig. 9, left). In the region of ideal maximum metastability identified in the parameter search (K = 3.5; see above), global metastability was significantly reduced in patients compared with controls (t84= −4.90, p < 0.0001; Fig. 9, right).
Simulated global metastability is decreased in traumatic brain injury patients compared with controls. Left, mean metastability (±1 SEM) in patients (n = 63) and controls (n = 26) in the computational model for a range of different coupling strengths (K). v = 11 ms−1. Right, Mean measures of metastability estimated within the model (±1 SEM) for 63 patients and 26 healthy control subjects suggest that global measures of metastability simulated from structural connectivity in TBI patients are reduced compared with healthy controls (*p < 0.05, **p < 0.01, ***p < 0.001, FDR corrected).
To further explore the effect of structural disconnection on simulated dynamics, we applied the same analysis to compute metastability within subsets of regions involved in canonical ICNs. Simulations using structural connectivity from individual TBI patients had significantly lower metastability within the dorsal attention network (t(84) = 4.15, p < 0.001), a right frontoparietal network (t(84) = −1.99, p < 0.05), a default mode network (t(84) = −3.75, p < 0.001), a salience network (t(84) = −3.62, p < 0.001), the primary auditory (t(84) = −4.06, p < 0.001), and low-level visual networks (t(84) = −2.45, p < 0.02). Simulated results were obtained from multiple regression using age and total gray matter volume as covariates of no interest and are FDR corrected for multiple comparisons (q < 0.1).
We assessed how well the simulated metastability predicts empirical metastability. Global empirical metastability was entered as the dependent variable of a linear regression, with global simulated metastability as a predictor variable and age, total gray matter volume, and the motion estimate for the empirical data as covariates of no interest. The overall model was able to significantly predict empirical measures of metastability (F(83) = 5.95, R2 = 0.24, p < 0.0001), with a significant effect of simulated metastability (t(83) = 3.30, p < 0.001)
To evaluate the extent to which global simulated metastability is determined by changes to macroscopic structural connectivity, we used linear regression with graph theoretical measures and group (patients or controls) as a covariate. This analysis showed that small-world index (t(84) = 2.91, p < 0.001), clustering coefficient (t(84) = 2.36, p < 0.05), and mean degree (t(84) = 5.96, p < 0.001) were significant positive predictors of simulated global metastability. An increase in characteristic path length was associated with reduced metastability (t(84) = −5.62, p < 0.001).
Discussion
We considered two interrelated questions: How does high-level behavior arise from the structural connectivity of the brain, and how does disruption of network structure alter behavior? Metastability has been suggested as a fundamental property of neural dynamics, serving as a conceptual bridge between brain structure and behavior (Tognoli and Kelso, 2014). Here, we used TBI as a model to interrogate the relationship between metastability (here defined as the SD of the Kuramoto order parameter), structural connectivity, and behavior. After TBI, metastability measured using fMRI is reduced compared with age-matched healthy control subjects. The level of metastability relates to behavioral impairment on a range of cognitive tasks. Importantly, using both empirical and computational modeling, we show that reduction in metastability after TBI is associated with damage to structural network topology, providing a demonstration of how metastable dynamics relate to behavior through structural connectivity.
DAI is a common pathology in TBI, accounting for much of the morbidity and mortality after injury, preferentially damaging long-distance tracts (Adams et al., 1989; Geddes et al., 1997). We demonstrate a significant reduction in the small-worldness of the structural connectome in TBI patients compared with controls, alongside a reduction in metastability. Previously, the relationship between network topology and metastability has been shown in computational simulations (Shanahan, 2010; Cabral et al., 2012; Wildie and Shanahan, 2012). Our computational findings, alongside empirical observations, provide additional support for a relationship between altered topology and metastable dynamics. We show that this relationship depends on the amount of structural damage (e.g., mean FA and average node strength), but it also depends on higher-order metrics such as clustering coefficient and small-worldness. The relationship between small-worldness of the structural connectome and metastability makes intuitive sense. Short overall path lengths facilitate increased global synchronization, while local modular architecture may provide some reservoir of different states, preventing the system from getting “stuck” in a synchronized state. However, the relationship between network topology and metastability may not be straightforward, involving heterogeneous time delays between nodes. Alternative descriptions of the network structure may be more effective at relating damage to altered metastability (e.g., the importance of a scale-free or rich-club structure; Senden et al., 2014); more computational and theoretical work is needed.
Our modeling results support the empirical findings, demonstrating that alterations in structural topology from TBI reduce simulated metastability. This suggests a mechanistic link between reduction in small-worldness and neural dynamics. The results of the simulations are consistent with computational models based on abstract network architectures (Friston, 1997; Shanahan, 2010) and those defined by anatomical connectivity (Deco et al., 2009b; Cabral et al., 2011). This work suggests how network topology allows the emergence of metastability, implicating sparseness (Friston, 1997) and small-worldness (Shanahan, 2010; Cabral et al., 2011, 2012; Wildie and Shanahan, 2012; Messé et al., 2014). More recently, networks with “rich-club” organization have been shown to support a broad repertoire of dynamic states (Senden et al., 2014). Such rich dynamics are reminiscent of the emergence of metastability; however, such a link is speculative. Future work could explore in more detail whether other graph theoretical properties, such as the presence of a strong rich club, better explain the alterations in metastability that we observed, leading to a more refined explanation of how complex neural dynamics emerge from the network topology of the brain.
Our fMRI empirical results were acquired at rest, in the absence of any explicit behavioral requirements. The rest state is when metastability is likely to be best suited to efficient cognitive flexibility. Rest, which must be distinguished from low arousal states such as sleep or sedation, can be thought of as a “jack-of-all-trades” state when the brain is in a broad exploratory regime. The dynamical regime during rest may constitute an upper limit for flexibility of the neural dynamics. We have previously shown, with both empirical and computational approaches, that metastability at rest is higher than during a focused cognitive task (Hellyer et al., 2014). During a task, high metastability is less desirable, since a specific configuration of brain systems is recruited (e.g., in coordinating specific visual and motor systems to perform a visually cued motor task). However, the dynamical regime during this task will still reflect the level of metastability at rest in terms of how the system can transition from rest to a task state efficiently and reliably. With low metastability at rest (e.g., after TBI), the system is likely to take longer and be less reliable, transitioning between cognitive states, showing a reduced repertoire of brain configurations required to facilitate task performance.
Metastability may be consistent with other descriptions of the brain as a dynamical system, such as self-organized criticality. Critical systems balance the competing demands of information propagation around a system with the need to maintain stable functional long- and short-scale functional relationships (Beggs and Plenz, 2003; Beggs, 2008). Such behaviors maximize information flow and capacity (Shew and Plenz, 2013), which is likely important for efficient cognitive function. Previous empirical work has provided evidence of metastable dynamics in systems demonstrating signatures of self-organized criticality (Haldeman and Beggs, 2005; Kitzbichler et al., 2009). A closely related question is whether metastable dynamics underlie cognitive function at finer spatiotemporal scales, in common with other measurements of brain activity and structure that show “scale-free” properties (Werner, 2007). Within the framework of self-organized criticality, an alteration in metastability may accompany a shift away from the critical state, associated with a decrease in efficiency of information storage, or processing capacity of the brain. Such a description has an intuitive link to cognitive deficits after TBI. The extent to which this change in metastable activity is invariant of scale (i.e., is present at the level of microscopic neural circuits as well as the macroscopic scale) is unclear. If metastability is indicative of a self-organized critical system, then it may be expected that the macroscopic dynamics changes described here may cascade across all spatial and temporal scales of the brain, from local neuronal circuits to the system-level description explored here. Examination of neural dynamics at a range of spatial and temporal scales using electrophysiological or optogenetic approaches in behaving animal models (Scott et al., 2014) may, therefore, be helpful in exploring these mechanisms in more detail.
We found reduction in metastability related to cognitive impairments on three tasks assessed: cognitive flexibility, speed of information processing, and associative memory. The switch cost index of the Trail Making Test, which involves rapidly and accurately switching between competing task demands, assesses cognitive flexibility and intuitively maps onto reduced metastability (which reflects reduced dynamical flexibility). This relationship may help explain perseveration after TBI: structural damage to white matter tracts limits the metastability of the brain, which limits cognitive flexibility. However, the relationship between metastability and behavior was not specific to cognitive flexibility, being present for the other two measures tested. All three tasks involve the integration of information across large-scale brain networks (Sharp et al., 2010a; Spreng and Grady, 2010; Bonnelle et al., 2011; Erika-Florence et al., 2014; Jilka et al., 2014). These tasks require communication between sensory, motor, and cognitive control regions, so altered global metastability may affect them all. This suggests large-scale metastability may be an important dynamical mechanism underlying general cognitive function. However, there may be different ways metastability can break down, evident as altered dynamics within specific brain networks, resulting in different profiles of impairments. The prominence of frontal network damage relating to metastability supports the clinically established association between cognitive inflexibility and frontal lesions. Future work with a larger group of patients, performing a broader range of behavioral tasks in the scanner, and imaging with higher temporal resolution would better characterize both across-subject variance in network metastability and associated patterns of behavioral impairment beyond the domains studied here. Information about altered neural dynamics could provide a sensitive biomarker to stratify patients and be used to design individualized treatments, involving electrical stimulation, pharmacological intervention, or neurofeedback (Sharp et al., 2014).
There are a number of limitations to the work. Our results show a relationship between metastability, individual cognitive performance, and underlying structure in the TBI patients rather than in both the patients and controls. This may in part be because we had limited behavioral data on the healthy control subjects; also, patients tend to be much more variable (in terms of behavioral impairment, neural dynamics, and structure), and therefore relationships may be easier to detect. However, we also acknowledge that the lack of relationships in the healthy controls mean we cannot claim that metastability is important for understanding individual differences in cognitive performance or structural network topology, beyond the patient group. A second limitation is based on the constraints inherent in tractography measured with diffusion MR. The structural connectivity matrices and graphs generated were undirected, in so much as feedforward and feedback connectivity of individual regions had a uniform effect on node–node functional interaction, which is unlikely to be the case in vivo. In addition, long-distance connections, e.g., interhemispheric pathways, may be difficult to resolve accurately since uncertainty in streamline location increases with the length of the tract (Jones, 2010a,b). The computational model, a system of coupled oscillators, is obviously a dramatic simplification of brain function. For example, the simulation is built on a relatively low-dimensional connectivity matrix of 164 regions. However, despite these limitations, the simulation provides important insights into the relationship between brain structure and function, broadly consistent with empirical findings. Such models, at least at the level of global network dynamics, replicate the broad changes in BOLD seen with fMRI, even though the model is based only on network topology (Deco et al., 2008, 2009b; Cabral et al., 2011, 2012; Messé et al., 2014), with no modeling of the functional specialization of individual nodes. These limitations mean that precise, quantitative comparisons between the simulations and the brain were not expected, although perform surprisingly well. Future work could consider whether our findings generalize to other computational oscillator models at comparable scales, such as the Wilson-Cowan model. Difficulties with the measurement of the BOLD fMRI signal, such as partial volume effects, regional differences in vascular reactivity or susceptibility artifacts would also make quantitative comparisons challenging. Finally, the data we present here were collected in the eyes-closed state. Recent work has suggested subjects may not consistently remain conscious throughout an extended resting scan (Tagliazucchi and Laufs, 2014). This opens up the possibility that subjects varied in whether they drifted in and out of sleep, and this could have affected empirical measurements of metastability. Therefore, future work should exclude this possibility, with the use of concurrent EEG to directly detect sleep states (Tagliazucchi and Laufs, 2014).
In summary, we found that large-scale structural disconnection is associated with reduced metastability, linked to impaired cognitive flexibility and other behavioral impairments. The link between damaged structural connectivity after TBI and reduced metastability (demonstrated both in empirical and simulated data) provides evidence that metastability is contingent on the integrity of the underlying structural network topology. This suggests a mechanistic link between structure, neural dynamics, and behavior. The results indicate a compelling link between brain structure and function and suggest that the framework of metastable dynamics offers an account for understanding the brain in health and disease.
Footnotes
P.J.H. was supported by a Medical Research Council Doctoral Training Award. G.S. was supported by a clinical research fellowship awarded in the Wellcome Trust-GlaxoSmithKline Translational Medicine Training Programme. We thank Frantisek Vasa for useful conversations on comparing simulated and empirical networks. We also thank the participants who were involved in this research.
The authors declare no competing financial interests.
- Correspondence should be addressed to Robert Leech, Computational, Cognitive, and Clinical Neuroimaging Laboratory, 3rd Floor, Burlington Danes, Hammersmith Hospital, Du Cane Road, London W12 0NN, UK. r.leech{at}imperial.ac.uk
