Tracking Whole-Brain Connectivity Dynamics in the Resting State

Abstract

Spontaneous fluctuations are a hallmark of recordings of neural signals, emergent over time scales spanning milliseconds and tens of minutes. However, investigations of intrinsic brain organization based on resting-state functional magnetic resonance imaging have largely not taken into account the presence and potential of temporal variability, as most current approaches to examine functional connectivity (FC) implicitly assume that relationships are constant throughout the length of the recording. In this work, we describe an approach to assess whole-brain FC dynamics based on spatial independent component analysis, sliding time window correlation, and k-means clustering of windowed correlation matrices. The method is applied to resting-state data from a large sample (n = 405) of young adults. Our analysis of FC variability highlights particularly flexible connections between regions in lateral parietal and cingulate cortex, and argues against a labeling scheme where such regions are treated as separate and antagonistic entities. Additionally, clustering analysis reveals unanticipated FC states that in part diverge strongly from stationary connectivity patterns and challenge current descriptions of interactions between large-scale networks. Temporal trends in the occurrence of different FC states motivate theories regarding their functional roles and relationships with vigilance/arousal. Overall, we suggest that the study of time-varying aspects of FC can unveil flexibility in the functional coordination between different neural systems, and that the exploitation of these dynamics in further investigations may improve our understanding of behavioral shifts and adaptive processes.

dynamics, fMRI, functional connectivity, independent component analysis, intrinsic activity, resting state, state variability

Introduction

Assessment of functional connectivity (FC) from functional magnetic resonance imaging (fMRI) time series, particularly during resting-state/task-free periods, has revealed a great deal of knowledge about the macro-scale spatiotemporal organization of the brain. Based on the correlations between intrinsic low-frequency oscillations (Biswal et al. 1995; Cordes et al. 2001) mediated by underlying structural connectivity (Honey et al. 2009), connectivity analysis has shifted focus away from merely localizing activations and deactivations and toward characterization of coactivation patterns, that is, network identification. A number of intrinsic connectivity networks (ICNs) are now widely recognized, in particular the default-mode network (Raichle et al. 2001; Buckner et al. 2008), ventral and dorsal attention networks (Corbetta and Shulman 2002; Fox, Corbetta et al. 2006; Vincent et al. 2008), and salience network (Seeley et al. 2007), and the relationships between them have been intensely studied in basic and clinical cognitive neuroscience. Recent work suggests a more refined and fine-grained parcellation of these large-scale networks into a multitude of smaller constituents (Kiviniemi et al. 2009; Abou-Elseoud et al. 2010; Allen et al. 2011), and also shows that these networks are not conditional upon a task-free resting state but are equally involved in task performance (Calhoun et al. 2008; Smith et al. 2009). Such substructure reveals the modular organization of different systems, with communication “hubs,” in graph theoretical terms (Hagmann et al. 2008; Buckner et al. 2009). This dramatically different view on aspects of brain function may in turn help improve diagnostic relevance for neuropsychiatric disorders, in particular where activation differences are subtle (Fornito and Bullmore 2012).

Despite such progress, we argue that the assessment of FC has been limited, in large part, by an implicit assumption of spatial and temporal stationarity throughout the measurement period. While this assumption is convenient, in that it keeps whole-brain connectivity analysis from becoming vastly more complex, it also unfortunately represents a gross oversimplification. Spontaneous/intrinsic fluctuations of activity and connectivity have long been appreciated in electrophysiological recordings of single cells, local fields, and surface electroencephalograms (EEGs). These fluctuations have been exploited in studies where high temporal resolution allows trial-by-trial exploration of the dynamics and adaptability of cognitive processes (Arieli et al. 1996; Makeig et al. 2004; Onton et al. 2006), and are also increasingly employed in single-trial analyses of task-based fMRI studies (e.g., Debener et al. 2006; Fox, Snyder et al. 2006; Eichele et al. 2008; Sadaghiani et al. 2009; Coste et al. 2011). Dynamics are potentially even more prominent in the resting state, during which mental activity is unconstrained. It is well established that individuals freely engage in several types of mental activity during resting periods (e.g., Delamillieure et al. 2010), and that the predominance of activity (e.g., imagery or inner language) affects FC and modular organization throughout the brain (Doucet et al. 2012). Relatively subtle modulations in cognitive load, for example, by instructing participants to keep eyes closed, open, or fixated, alter the spectral content of spontaneous activity and patterns of FC throughout subcortical nuclei, sensorimotor cortex, and default-mode regions (McAvoy et al. 2008; Yan et al. 2009; Wu et al. 2010). More straightforward attempts to modulate internal activity, such as requesting that individuals remember the events of their day or silently recall song lyrics, also result in pronounced changes in whole-brain FC, with differences sufficiently large and robust to permit highly accurate classification of cognitive states (Shirer et al. 2012). Finally, explicit investigations of resting-state FC dynamics have unambiguously demonstrated the time-varying nature of both connectivity strength and directionality (i.e., positive or negative) (e.g., Chang and Glover 2010; Kiviniemi et al. 2011; Hutchison et al. 2012), with indications that capturing this variability may engender new understanding of the FC differences found in neuropsychiatric diseases such as Alzheimer's dementia (Jones et al. 2012), autism (Starck et al. 2012), and schizophrenia (Sakog˘lu et al. 2010).

Considering the strong evidence for fluctuations in FC, how then should we best investigate spontaneous variations in the framework of a large group study? In contrast to investigations with task designs or paced experimental manipulations, the variety of mental states experienced during rest and the transitions between them cannot be specified a priori. Stable connectivity patterns and change points must be learned from the data directly, and formal models for this process are just being developed (Cribben et al. 2012).

In this article, we describe a simple data-driven approach to assess FC dynamics based on established techniques, including spatial independent component analysis (ICA), sliding time-window correlation, and k-means clustering of windowed correlation matrices. Group ICA (Calhoun et al. 2001) is used to decompose multisubject resting-state data into functionally homogeneous regions (Kiviniemi et al. 2009; Abou-Elseoud et al. 2010), enabling a whole-brain analysis without resorting to atlas-based regions of interest that may merge distinct areas (e.g., see Shirer et al. 2012) or fail to capture intersubject spatial variability (Allen et al. 2012). Time-varying FC is estimated by computing correlations between component time courses (TCs; Jafri et al. 2008) using a series of sliding windows (Sakog˘lu et al. 2010). We then evaluate the degree of variability in resulting FC time series to identify brain regions with particularly variable (or flexible) connections. Lastly, we use k-means clustering to identify patterns of FC that reoccur in time and across subjects (Lloyd 1982). We describe these clusters as “FC states” in a conceptual analogy to EEG microstates, short periods during which scalp topography remains quasi-stable (Lehmann 1990; Pascual-Marqui et al. 1995). The FC states observed here are highly replicable and in part diverge strongly from stationary connectivity patterns, challenging current descriptions of interactions between large-scale networks. Moreover, the differential occurrence of specific FC states over time motivates theories regarding their functional roles and relationships with arousal. We conclude that the study of time-varying aspects of FC can unveil underappreciated flexibility in the functional coordination between different neural systems, and that the further investigation and exploitation of these fluctuations may improve our understanding of cognitive and behavioral dynamics.

Materials and Methods

Data Acquisition and Preprocessing

Data used in this work comprise resting-state scans from 405 healthy participants (mean age: 21.0 years, range: 12–35 years, 200 females) collected on the same scanner and combined across 34 studies and 18 principal investigators at the Mind Research Network. Informed consent was obtained from all subjects according to institutional guidelines at the University of New Mexico and data were anonymized prior to group analysis. The current dataset represents a subset of the 603 subjects used in Allen et al. 2011, where more stringent inclusion criteria have been incorporated to limit the influence of motion (subject data with a maximum translation of >1.5 mm were excluded) and poor spatial normalization (subject data with spatial correlation to EPI template <0.93 were excluded), and to improve sample homogeneity (subjects older than 35 were excluded). Further details on sample demographics can be found in Allen et al. 2011.

Imaging was performed on a 3-T Siemens Trio scanner with a 12-channel radio frequency coil. T₂*-weighted functional images were acquired using a gradient-echo EPI sequence with TE = 29 ms, TR = 2 s, flip angle = 75°, slice thickness = 3.5 mm, slice gap = 1.05 mm, field of view = 240 mm, matrix size = 64 × 64, voxel size = 3.75 mm × 3.75 mm × 4.55 mm. Resting-state scans were a minimum of 5 min, 4 s in duration (152 volumes); any additional volumes were discarded to match data quantity across participants. Subjects were instructed to keep their eyes open during the scan and fixate on a foveally presented cross.

Functional images were preprocessed using an automated pipeline based around SPM 5 (http://www.fil.ion.ucl.ac.uk/spm/software/spm5). Preprocessing included the removal of the first 4 image volumes to avoid T1 equilibration effects, realignment using INRIalign, slice-timing correction using the middle slice as the reference frame, spatial normalization into Montreal Neurological Institute space, reslicing to 3 mm × 3 mm × 3 mm voxels, and smoothing with a Gaussian kernel (FWHM = 5 mm). Voxel time series were z-scored to normalize variance across space, minimizing possible bias in subsequent variance-based data reduction steps. Note that variance normalization may not be optimal for all investigations as it removes amplitude information that may be relevant for differences between subjects; however, our focus here is temporal modulation rather than magnitude, thus variance normalization is preferred.

Group ICA and Postprocessing

Data were decomposed into functional networks using a group-level spatial ICA as implemented in the GIFT toolbox (http://mialab.mrn.org/software/gift/; Fig. 1A). We used a relatively high model order (number of components, C = 100) to achieve a “functional parcellation” of refined cortical and subcortical components corresponding to known anatomical and functional segmentations (Kiviniemi et al. 2009; Smith et al. 2009; Abou-Elseoud et al. 2010). Subject-specific data reduction via principal components analysis (PCA) retained 120 principal components using a standard economy-size decomposition and group data reduction retained C = 100 PCs using the expectation–maximization (EM) algorithm to avoid otherwise prohibitive memory requirements (Roweis 1998). The Infomax ICA algorithm (Bell and Sejnowski 1995) was repeated 10 times in Icasso (http://www.cis.hut.fi/projects/ica/icasso) and aggregate spatial maps (SMs) were estimated as the modes of the component clusters. Subject-specific SMs (S_i) and time courses (TCs, R_i) were estimated using the GICA1 back-reconstruction method based on PCA compression and projection (Calhoun et al. 2001; Erhardt et al. 2011), which is equivalent to least-squares based dual regression (Filippini et al. 2009) in the subject-specific PC-reduced space (see Appendix of Allen et al. 2012). As in Allen et al. 2011, we characterized a subset of C₁ = 50 components as ICNs, as opposed to physiological, movement related, or imaging artifacts (ARTs). Components were evaluated based on expectations that ICNs should exhibit peak activations in grey matter, low spatial overlap with known vascular, ventricular, motion, and susceptibility artifacts, and should have TCs dominated by low-frequency fluctuations (Cordes et al. 2000).

Figure 1.

Illustration of analysis steps. (A) Group ICA decomposes resting-state data from M = 405 subjects into C = 100 components, C₁ = 50 of which are identified as intrinsic connectivity networks (ICNs). GICA1 back reconstruction is used to estimate the TCs (R_i) and SMs (S_i) for each subject. (B) Stationary FC between components (left, ∑_i) is estimated as the covariance of R_i. Dynamic FC (right, ∑_i^L1(w)) is estimated as the series of regularized covariance matrices from windowed portions of R_i.

Component TCs underwent additional postprocessing to remove remaining noise sources. These include low-frequency trends related to scanner drift, motion-related variance which may not be wholly captured in distinct components given the spatial nonstationarity inherent to movement, and other nonspecific “spikes” or noise artifacts that are not decomposed well by a linear mixing model. Postprocessing included 1) detrending linear, quadratic, and cubic trends, 2) multiple regression of the 6 realignment parameters and their temporal derivatives, 3) removal of detected outliers, and 4) low-pass filtering with a high-frequency cutoff of 0.15 Hz. The outlier removal approach used here is similar to the “scrubbing” method proposed by (Power et al. 2012), but rather than remove affected time points from data (which would compromise the subsequent sliding window approach), we replaced outliers with the best estimate using a third-order spline fit to the clean portions of the TCs. Outliers were detected based on the median absolute deviation, as implemented in 3Ddespike (http://afni.nimh.nih.gov/afni). Improvement in the root-mean-square of the temporal derivative over component TCs, referred to as “DVARS” in Power et al. 2012 (see SupplementarySupplementary DataA) and removal of anticipated motion-related biases from FC estimates (see Supplementary Fig. S1B), suggest satisfactory correction of motion artifacts. As a final step in postprocessing, we normalized the variance of each TC, thus covariance matrices (below) correspond to correlation matrices. In exploratory work, we repeated all analyses on component TCs that underwent minimal postprocessing (only detrending and low-pass filtering) and found nearly identical results with regard to FC temporal variability (Fig. 4) and connectivity states (Fig. 5), suggesting that the somewhat aggressive postprocessing applied here did not fundamentally alter dynamic structures.

FC Estimation and Temporal Variability

For each subject i = 1 … M, stationary FC was estimated from the TC matrix R_i as the C × C sample covariance matrix ∑_i (Fig. 1B, left). Dynamic FC was estimated with a sliding window approach, wherein we computed covariance matrices ∑_i(w), w = 1 … W, from windowed segments of R_i (Fig. 1B, right). We used a tapered window (see Fig. 1B, right), created by convolving a rectangle (width = 22 TRs = 44 s) with a Gaussian (σ = 3 TRs) and slid in steps of 1 TR, resulting in W = 126 windows. Because relatively short time segments may have insufficient information to characterize the full covariance matrix, we estimated covariance from the regularized precision matrix (inverse covariance matrix, ∑_i⁻¹(w)) (Varoquaux et al. 2010; Smith et al. 2011). Following the graphical LASSO method of Friedman et al. 2008, we placed a penalty on the L1 norm of the precision matrix to promote sparsity. The regularization parameter lambda (λ) was optimized separately for each subject by evaluating the log-likelihood of unseen data (windowed covariance matrices from the same subject) in a cross-validation framework. Final dynamic FC estimates for each window, ∑_i^L1(w), were concatenated to form ∑_i^L1, a C × C × W array representing the changes in covariance (correlation) between components as a function of time. Both stationary and dynamic FC estimates were Fisher transformed to stabilize variance prior to further analysis.

FC estimates between some ICNs exhibited greater temporal variability than others (see Figs. 3 and 4A). We used a simple algorithm to separate ICNs into groups with more variable FC (Partition 1, P₁), referred to as the “zone of instability” (ZOI), and less variable FC (P₂). The algorithm proceeded with 3 steps: 1) ICNs were randomly assigned to P₁ or P₂ with equal probabilities, 2) in repeated iterations, membership for a single component was changed in such a way to maximize the Separation Index

$$\hbox{SI} = \left( {\displaystyle{1 \over {h_1 }}\mathop \sum \limits_{u,v\; \in \; P_1 } a_{u,v} - \displaystyle{1 \over {h_2 }}\mathop \sum \limits_{u,v\; \in \; P_2 } a_{u,v} } \right)/{\sigma}_2 $$

where a_u,v is the average low-frequency (<0.025 Hz) amplitude of FC oscillations between components u, v = 1, ... C₁, σ₂ is the standard deviation over a_u,v,|$u,v \in \hbox{}P_{\rm 2} $|⁠, and h₁, h₂, are the number of components in each respective partition, and 3) stopping criteria were reached when any change in membership would result in a decrease in SI. To obtain a robust partitioning solution that incorporated data variability and was independent of initial conditions, we repeated the algorithm on b = 1000 bootstrap resamples of the data, that is, M subjects were drawn with replacement and a_u,v was recomputed as the average over that sample. ZOI scores for each ICN (see Fig. 4B) were then calculated as the fraction of repetitions in which the component was assigned to P₁.

Clustering Analysis

To assess the frequency and structure of reoccurring FC patterns we applied the k-means clustering algorithm (Lloyd 1982) to windowed covariance matrices. We used the L1 distance function (Manhattan distance), guided by work suggesting the L1-norm may be a more effective similarity measure than the L2 (Euclidean) distance for high-dimensional data (Aggarwal et al. 2001). Only covariances between the C₁ = 50 ICNs were used in the clustering analysis, resulting in (50 × (50 − 1))/2 = 1225 features. Prior to clustering, subject arrays ∑_i^L1 were subsampled along the time dimension. Subsampling was performed both to reduce redundancy between windows (the chosen time step of 1 TR induces high autocorrelation in FC time series) and to reduce computational demands. Similar to EEG microstate analysis (Pascual-Marqui et al. 1995), subject exemplars were chosen as those windows with local maxima in FC variance, resulting in 7.5 ± 1.5 (mean ± SD) windows per subject (range: 4–12) for 3026 instances. The clustering algorithm was applied to the set of all subject exemplars and was repeated 500 times to increase chances of escaping local minima, with random initialization of centroid positions. The resulting centroids (cluster medians) were then used to initialize a clustering of all data (405 subjects × 126 windows = 51 030 instances). To ensure that exemplar selection did not bias group clusters, we repeated the analysis using alternative methods: 1) we first applied clustering at the subject level and propagated subject centroids to the group level, and 2) we selected 6 windows from each subject at random. Both methods produced clusters almost identical to those observed using windows at local maxima. Likewise, we repeated the clustering using different distance functions (correlation, Euclidean, and cosine, rather than the L1-norm) and also found extremely similar results. For group clustering (and subject-level clustering where applicable), the number of clusters (k) was determined using the elbow criterion of the cluster validity index, computed as the ratio between within-cluster distance to between-cluster distance, though additional exploratory analyses using hierarchical clustering or explicitly varying k (see Supplementary Fig. S4) demonstrated consistent results over a large range of k.

Reproducibility of clusters was established via replication on bootstrap resamples and nonoverlapping split-half samples of subjects (see SupplementarySupplementary Data). We additionally validated that clusters were not due to nonspecific differences across subjects or time by applying k-means to surrogate data. Surrogate datasets were created via phase randomization in the Fourier domain (Prichard and Theiler 1994) (see SupplementarySupplementary DataA). To create surrogate dataset 1 (SR1), each subject array ∑_i^L1 was Fourier transformed and the same random sequence of phases was added to all FC phase spectra, maintaining the covariance structure across all windows. For SR2, an identical process was applied, but a different random sequence of phases was added to each FC phase spectrum, disturbing the covariance structure. The mean, variance, and temporal autocorrelation of FC time series in SR1 and SR2 were identical to the original data (see SupplementarySupplementary DataB). As seen in Supplementary Figure S6C, clusters were found in SR1, but not SR2, demonstrating that it is the co-occurrence of FC between specific ICNs that drives the clustering, rather than distinctions in FC mean or variance across subjects.

As a validation of the clustering approach, we used the SimTB framework (Erhardt et al. 2012) (http://mialab.mrn.org/software/simtb/) to simulate fMRI time series under a model of dynamic neural connectivity, then applied k-means clustering to estimate FC states from windowed covariance matrices in a manner identical to the real data. As shown in Supplementary Figure S7, the clustering provided excellent estimates of both the discrete neural states and the transitions between states, suggesting that the clusters derived from real data faithfully reflect the structure and temporal properties of dynamic connectivity.

Results

Figure 2A displays the ICNs identified with group ICA. Based on their anatomical and presumed functional properties, ICNs are arranged into groups of subcortical (SC), auditory (AUD), somatomotor (SM), visual (VIS), cognitive control (CC; referring loosely to the planning, monitoring, and adapting one's behavior), default-mode (DM), and cerebellar (CB) components. The manual arrangement of ICNs is very similar to various orderings provided by empirical methods, including spectral clustering and algorithms based on the optimization of modularity and diagonal structure as implemented in Brain Connectivity Toolbox (http://www.brain-connectivity-toolbox.net/). Detailed images of each ICN are shown in Supplementary Figure S2 and coordinates of peak activations are listed in Supplementary Table S1. ICNs are similar to those observed in previous high model order ICA decompositions (Kiviniemi et al. 2009; Smith et al. 2009; Abou-Elseoud et al. 2010; Allen et al. 2011) and cover the majority of subcortical and cortical gray matter. Figure 2B displays the FC between ICNs, computed over the entire scan length and averaged over subjects. Patterns of FC are consistent with prior literature, showing modular organization within sensory systems and default-mode regions, as well as anticorrelation between these domains (e.g. Fox et al. 2005; Chang and Glover 2010; Shirer et al. 2012). We note that, based on average connectivity, language areas (L MTF + IFG and R cerebellum) cluster with default-mode components, presumably because of a high proportion of time spent in self-narrative and inner speech (Delamillieure et al. 2010).

Figure 2.

ICN SMs (A) and the stationary FC between them (B). ICNs are divided into groups and arranged based on their anatomical and functional properties. Within each group, the color of the component in (A) corresponds to the colored flag shown along the axes of (B). FC was averaged over all subjects and inverse Fisher transformed (r = tanh(z)) for display, facilitating comparisons with previous studies. ICN labels in (B) denote the brain region with peak amplitude and refer to bilateral activations unless specified as left (L) or right (R). See Supplementary Figure S2 and Table S1 for more detailed information on each component. STG, superior temporal gyrus; PreCG, precentral gyrs; PoCG, postcentral gyrus; SMA, supplementary motor area; ParaCL, paracentral lobule; SPL, superior parietal lobule; MTG, middle temporal gyrus; FFG, fusiform grys; MOG, middle occipital gyrus; SOG, superior occipital gyrus; IPL, inferior parietal lobule; ITG, inferior temporal gryus; MCC, middle cingulate cortex; pInsula, posterior insula; MiFG, middle frontal gyrus; IFG, inferior frontal gyrus; aInsula, anterior insula; PHG, parahippocampal gyrus; PCC, posterior cingulate cortex; AG, angular gyrus; ACC, anterior cingulate cortex; SFG, superior frontal gyrus; CB, cerebellum.

Results from the sliding window analysis are shown in Figure 3A–C, which displays FC dynamics for 3 representative subjects. As seen from the FC time series (Fig. 3A2–C2) and corresponding videos (Supplementary Movies 1–3), FC between some ICNs is highly nonstationary, in some cases exhibiting both strongly positive and strongly negative correlations within the 5-min scan. Fourier analysis of the time series (Fig. 3A3–C3) shows low-frequency FC oscillations peaking between 0.005 and 0.015 Hz, corresponding to a period on the order of 100 s. We note that FC oscillations between ICN components are significantly larger than those between ART components (see Supplementary Fig. S3A, spectral area under the curve: paired t(404) = 39.4, P ≈ 0), and are also weighted toward lower frequencies (see SupplementarySupplementary DataB, spectral center of mass: paired t(404) = −27.7, P ≈ 0). These distinctions suggest that dynamics between ICNs are related to changes in connectivity, rather than nonspecific physiological changes (e.g., cardiac or respiratory shifts) or subject movement that would be expected to affect all components similarly.

Figure 3.

Examples of FC dynamics for subject 124 (A), subject 267 (B) and subject 360 (C). (A1–C1) FC for each subject, averaged over all windows. (A2–C2) FC time series for connections between select pairs of ICNs. Correlation coefficients are plotted at the time point corresponding to the center of the window. Top panels show ∑_i^L1(w) for select windows. Highlighted connections are PreCG [2] to Thalamus [15] (light blue), L MOG [89] to R PoCG [10] (red), L IPL [76] to MOG [80] (orange), ACC [26] to R IPL [67] (dark blue), and MiFG + SFG [48] to L AG [75] (green). Highlighted windows are a subsample of the exemplars used in the clustering analysis (see Fig. 5A). (A3–C3) FC spectra for the time series in (A2–C2). Filled colored arrows marking the FC element locations in (A1–C1) correspond to the colored lines in (A2–C2) and (A3–C3).

Focusing on ICNs, we observe greater FC variability between some pairs of components than others. For example, for the subjects shown in Figure 3, FC between DM components 48 and 75 (green line) is relatively stable throughout the scan, while FC between CC component 76 and VIS component 80 (orange line) undergoes larger fluctuations. Comparing the amplitude of low-frequency FC oscillations between all pairs of ICNs (Fig. 4A) suggests an interesting pattern: a distinct set of ICNs that have more variable FC with one another than with other components. ICNs centered in superior occipital cortex exhibit more variable FC with default-mode regions than other areas, and both these domains show variable connections with components of inferior parietal cortex. We use an iterative partitioning algorithm on repeated bootstrap resamples of the data to distinguish the set of components that show more variable connections with one another (Fig. 4B, see Materials and Methods section for details). As expected from our initial examination of oscillation amplitude (Fig. 4A), dorsal attention areas, default-mode regions, and superior occipital networks are consistently assigned to the partition with more variable connectivity (Fig. 4C). We refer to this set of regions as comprising a ZOI, potentially reflecting a large-scale network capable of flexible behavior and heterogeneous function.

Figure 4.

Assessment of FC variability. (A) Amplitude of low-frequency (<0.025 Hz) FC oscillations between ICNs, averaged over subjects. Bins with greater amplitude indicate more variable FC. (B) Bootstrap partitioning procedure used to identify zone of instability (ZOI) scores for each component (see Materials and Methods section) (C) Surface rendering of ICNs with a ZOI score of >0.5.

Connectivity States

As seen in Figure 3A2–C2 and Supplementary Movies 1–3, the variability in individual FC time series is hardly random. Fluctuations give rise to highly structured patterns of FC that emerge and dissolve over tens of seconds. To explore the possibility that certain connectivity patterns may be quasi-stable, that is, they reoccur over time and are present in numerous subjects, we applied k-means clustering to the windowed FC matrices (see Materials and Methods section). Figure 5 displays the clustering results with k = 7. Each matrix represents the centroid of a cluster and putatively reflects a connectivity state stably present within the data. We note that these clusters are fully reproducible in bootstrap resamples of subjects and split-half analyses (see SupplementarySupplementary DataA), are decomposed over a large range in the number of clusters (k) (see SupplementarySupplementary Data), and are not found when using “null” data (i.e., disrupted covariance structure) that is matched to the real data in terms of mean, variance, and temporal autocorrelation (see SupplementarySupplementary DataA).

Figure 5.

Clustering approach (A) and result (B) for k = 7. Each cluster (States 1–7) is summarized with its centroid (left), modularity partition obtained with the Louvain algorithm for finding community structure (top right), and number of occurrences as a function of time (bottom right). The total number and percentage of occurrences is listed above each centroid and the number of modules (n) and modularity index (Q*, as defined in Rubinov and Sporns 2011) are adjacent to module depictions. Where possible, module colors (blue, red, green, and yellow) were matched across states such that similar partitions have the same color. As modularity partitions vary slightly from run-to-run, the Louvain algorithm was repeated on 100 bootstrap resamples (resampling ∑_i^L1(w) within each cluster) and consistency in modular assignment was mapped to color opacity (completely opaque = assigned to same module on all resamples; completely transparent = assigned to same module on 1/n resamples). beta Values indicate the slope (in units of percentage per minute) of the best linear fit (red) to the occurrence trend (blue). Light gray lines show occurrence profiles for 100 bootstrap resamples (resampling subjects).

FC states are arranged in order of emergence (see Supplementary Fig. S4). State 1, which accounts for >30% of all windows, resembles the stationary FC (compare to Fig. 2) and, based on our observations from higher order and hierarchical clustering, signifies the average of a large number of additional states that are not sufficiently distinct or frequent to be separated. FC patterns in States 2–7 are observed less frequently (ranging between 7 and 15%), but represent connectivity diverging substantially from the mean. Differences are apparent in terms of both the magnitude and the sign of connectivity between ICNs and the network modularity, that is, the partitioning of ICNs into subgroups. Modular partitions (color-coded in the axial, coronal, and sagittal slices in Fig. 5) are found using the Louvain algorithm with the definition for optimal modularity suggested by (Rubinov and Sporns 2011) and implemented in the Brain Connectivity Toolbox, which operates on the unthresholded FC matrix.

We describe 3 notable features that differ between FC states, though many additional distinctions are also present. First, states are differentiated by connectivity between DM regions. In States 2 and 7, a large DM module (blue) is present, comprising ICNs dedicated to the bilaterial hippocampi, precuneus, posterior cingulate cortex (PCC), anterior cingulate cortex, angular gyri, lateral temporal cortex, and inferior frontal cortex. In contrast, States 5 and 6 reveal a functional segregation between posterior DM nodes (precuneus and PCC) and anterior and lateral parietal regions. ICNs covering posteromedial cortex act in synchrony with some VIS and CC components and show weak asynchrony with other DM regions, leading to modified module affiliations. States are also differentiated by the FC between DM components and other regions throughout the brain. In most states, SM components show negative correlations with the DM system and positive correlations with sensory regions, often forming a large sensorimotor module (red). However in States 6 and 7, several SM ICNs exhibit “positive” correlations with DM components (and negative correlations with visual, auditory, and subcortical regions), resulting in the inclusion of these motor components in the DM module. Similarly, note that several frontal CC components that are typically are anticorrelated or uncorrelated with the DM system (e.g., States 2, 6, 5, and 7) exhibit different behavior in State 4 and show positive FC with DM regions (and negative FC with sensorimotor areas). The DM module expands to include a larger number of frontal regions that form their own module in other states (e.g., yellow in States 1 and 5; green in States 2 and 6). A third discriminating feature is the FC between cortical and subcortical components, which largely relates to State 3. With the exception of State 3, the thalamus and parts of the cerebellum consistently show synchronous activation with cortical components dedicated to the sensorimotor system. In striking contrast, State 3 shows strongly asynchronous activation between subcortical regions (bilateral thalamus, basal ganglia, and cerebellum) and sensorimotor cortex. This dissociation is accompanied by a substantial reduction in connectivity between DM regions, such that a prototypical DM module is no longer present. These large deviations in FC and modularity suggest that State 3 represents a cognitive state quite distinct from those represented in other clusters, a topic we explore further in the analysis below and in the Discussion section.

State Transitions

In addition to the describing the connectivity differences that distinguish FC states, we can also examine their occurrence as a function of time (Fig. 5, lower right panels) and the transitions between them. Figure 6A shows the state assignments as a function of time for the 3 example subjects. As we would expect from the very slow connectivity dynamics (Fig. 3), FC tends to be assigned to single states for long periods, though shorter periods are commonly seen when transitioning between states (e.g., for Subject 267, FC is briefly assigned to State 5 as it changes from States 1–6). Similar behavior is seen in iterations of dynamic FC simulations (see SupplementarySupplementary Data), with misclassification most likely to occur in transition periods. We can characterize transition behavior by considering FC as a Markov chain (MC), a system that undergoes transitions between a discrete number of states. In Figure 6B, we show the average transition matrix (TM) for this system, which represents the probability of changing from one state to another. White squares along the diagonal signify a very high probability of staying in the same state. For the off-diagonal elements, hotter colors in the State 1 column indicate a higher probability of entering State 1 from the other states, and cooler colors in the State 3 row indicate a lower probability of exiting State 3. Because the average TM is regular (all transition probabilities are >0), we can also approximate the stationary probability vector (π) as the principal eigenvector of the TM (Meyer 2000). The vector π, displayed in Figure 6C, represents the equilibrium probability of finding the MC in a particular state, that is, the expected behavior of the system in the long run; note that this is a distinct quantity from the number or fraction of occurrences presented in Figure 5. The stationary probability for State 3 exceeds 0.4, far greater than the probabilities for other states, which range from roughly 0.05 to 0.15, meaning that in the long run, the system is most likely to be found in State 3. These results are congruent with the temporal trends in state assignment displayed in Figure 5. There is a steady increase in the number of windows clustered into State 3 over time (increasing 3-fold), and a corresponding decrease in the frequency of windows assigned to States 2 and 7. Based on transition dynamics, we hypothesize that State 3 represents a state of drowsiness or light sleep, which subjects are more likely to enter into as time continues. Supporting this hypothesis, both the thalamocortical disconnection and weakening of DM connectivity that distinguish State 3 are consistent with changes to FC that occur as subjects move from wakefulness to sleep (Spoormaker et al. 2010; Larson-Prior et al. 2011; Sämann et al. 2011)

Figure 6.

Transitions between FC states. (A) State vectors for the 3 example subjects shown in Figure 3. Assigned states are plotted at the time point corresponding to the center of the sliding window. (B) The state transition matrix (TM), averaged over subjects. High values along the diagonal indicate a high probability of staying in a state. Note that transition probability is color-mapped on a log-scale. (C) The stationary probability vector (π, principal eigenvector of the TM) shows the steady state, or “long run” behavior. Error bars indicate the nonparametric 95% confidence intervals obtained from 1000 bootstrap resamples of the average TM (resampling subjects).

Discussion

Here, we explore dynamic patterns of FC with ICA, sliding windows, and clustering. Adding to the growing literature on connectivity dynamics during rest (Chang and Glover 2010; Kiviniemi et al. 2011; Hutchison et al. 2012; Jones et al. 2012) and internally driven states (Harrison et al. 2008; Shirer et al. 2012), our analysis of connectivity dynamics in a large sample (n = 405) provides, to our knowledge, the first whole-brain characterization of regional differences in FC variability and distinction of discrete FC states. These results and their implications are discussed in turn.

A Zone of Instability

We use a partitioning algorithm to identify a group of ICNs with more variable connections between them, conceived as a ZOI. This approach distinguishes heteromodal occipital cortex, lateral parietal cortex associated with the dorsal attention system (Corbetta and Shulman 2002), and lateral and medial aspects of the default network (Buckner et al. 2008). Perhaps unsurprisingly, ZOI regions are of the most globally connected (Cole et al. 2010) and some consistently emerge as structural and functional hubs (Hagmann et al. 2008; Buckner et al. 2009; Honey et al. 2009) or regions of high network centrality (Zuo et al. 2012) in topological descriptions of brain connectivity, collectively suggesting heterogeneous and integrative function. The characterization of a ZOI is somewhat at odds previous work delineating these same regions into 2 (Fox et al. 2005) or 3 (Vincent et al. 2008) different systems with functionally distinct, and potentially competing roles. While ZOI components certainly participate in systems with different primary roles, our analyses highlights that the relationships between them are flexible (more so than between other regions) with functional connections that emerge and dissolve, arguing against a labeling scheme with fixed segregation. Furthermore, our results may motivate analyses focused on a sort of “FC variability mapping” for both resting state and task datasets. Rather than focus on regions that show relatively constant patterns of connectivity (implicit to conventional connectivity analyses), there may be great benefit to identify those regions whose FC is notably more variable, suggesting membership in multiple large-scale systems and roles in adaptive processes.

Connectivity States

We use k-means clustering to identify reoccurring short-term connectivity patterns, which we describe as FC states. FC states are well predicted by large-scale models of neuronal connectivity that consider the repertoire of functional motifs generated by a given structural architecture. As described by Deco et al. (2011) in their recent review of connectivity models and dynamics, “[m]easurements over longer time windows recapitulate the anatomical connectivity, reflecting the RSNs that have been characterized in the literature … Shorter time windows emphasize the small departure from the RSN pattern, in which different nodes form new functional networks for a short period of time and then return to the RSN pattern. […T]he deterministic structure provided by the anatomy allows for certain functional networks to be repeated frequently in time, but that at any given point the precise configurations depend on the part of the dynamic repertoire that is being explored … ”. In excellent agreement with these models, we identify stable “departures” that in some cases are strikingly different from FC characterized over long time scales. Focusing specifically on DM connectivity, the FC states shown in Figure 5 challenge the notion of a singular and stable default-mode network. Rather, the constituents of this network and their covariation with nondefault-mode regions, such as motor and premotor cortex, are variable over time. Of particular note, precuneus ICNs exhibit affiliation with the DM module in only some states (1, 2, and 7), reflecting and potentially resolving the confusion and ongoing discussion regarding the inclusion of the precuneus in the default network (see Buckner et al. 2008 and Margulies et al. 2009 for excellent discussions). Furthermore, established “core” DM regions such as posterior cingulate and lateral parietal cortex may temporarily fail to exhibit synchronous activity, as in State 3. These findings, along with previous investigations of dynamics (Chang and Glover 2010; Kiviniemi et al. 2011; Hutchison et al. 2012) highlight that characterizations of connectivity and functional networks are strongly shaped by the time scale over which connectivity is measured. As studies increasingly explore the relationship between behavior and connectivity patterns (e.g., Kelly et al. 2008; Hamilton et al. 2011; Thompson et al. 2012), it will become more important to understand the relevant time scales over which FC (and FC changes) should be assessed.

Though our study is the first to characterize discrete FC states from resting-state imaging data, the observed states relate well to previous findings. In their study of connectivity dynamics in macaques and humans, Hutchison et al. (2012) describe periods of “hypersynchronization” found in both species during which oculomotor and motor networks showed extremely high intranetwork connectivity, that is, strong and positive correlations between all nodes. Hypersynchronization periods persisted for tens of seconds and reoccurred throughout the scan, suggesting analogy to the FC states described here. Specifically, States 2 and 6 represent periods with synchronous activation between oculomotor/dorsal attention areas, while State 4 signifies time windows with high correlations throughout the motor system. Though not investigated by (Hutchison et al. 2012), our results predict that periods of hypersynchronization between motor or oculomotor nodes would be accompanied by synchronization of default-mode regions and strong anticorrelation between the systems, highlighting the benefit of whole-brain connectivity studies.

As mentioned in the Results section, we posit that State 3 represents connectivity related to drowsiness/light sleep based on the temporal properties of occurrence and whole-brain FC. Because the resting-state scan is unconstrained and subjects are free to “mind-wander” in any fashion they choose, it is unlikely that FC patterns representing specific cognitive states would follow similar timelines across subjects. In contrast, the only feature expected to be common across subjects is the increased likelihood of drowsiness or sleep. Intriguingly, State 3 exhibits the expected temporal profile for drowsiness, becoming increasingly frequent as the scan progresses (Fig. 5). Corresponding results from a MC analysis suggest that subjects are most likely to be found in this state in the long run (Fig. 6). Along with increased prevalence in time, the FC pattern for State 3 is consistent with a transition from wakeful rest to sleep. The descent to sleep is marked by reduced thalamocortical connectivity (Spoormaker et al. 2010), increased subcortical connectivity (Larson-Prior et al. 2011), and a breakdown of default-mode connectivity (Spoormaker et al. 2010; Larson-Prior et al. 2011), all of which are found in State 3. Also in agreement with this theory, FC states showing the largest decreases in occurrence over time (States 2 and 7) are those with extensive intradefault network connectivity and greater antagonism between DM and CC components, signatures previously associated with greater task performance and presumably awareness (Kelly et al. 2008; Thompson et al. 2012). Though the hypothesized representation of State 3 should be tested with concurrent EEG-fMRI recordings, we believe that its identification in the current dataset speaks to the strength of the dynamic estimation/clustering approach and supports accessibility of spontaneous state-transitions from imaging data alone.

Finally, it is interesting to consider the relationship between FC states and EEG microstates, both of which are proposed to reflect the coordination of large-scale neural assemblies supporting different cognitive processes (e.g., Lehmann et al. 1998; Shirer et al. 2012). Although microstates persist for only hundreds of milliseconds, their dynamics exhibit fractal (scale free) organization that spans time scales as long as tens of seconds (Van De Ville et al. 2010). Thus, it is possible that FC states and EEG microstates capture very similar physiological phenomena, albeit as seen through very different spatial and temporal filters. Indeed, several recent studies have pursued microstate-based fusion of EEG-fMRI resting-state datasets, with promising results suggesting some correspondence between microstates and a limited number of fMR-based ICNs (Britz et al. 2010; Musso et al. 2010; Yuan et al. 2012). Large-scale ICNs have also been linked to covarying amplitude fluctuations in alpha and beta neural oscillations between distant regions (Brookes et al. 2011), with emerging work focused on the time-varying aspect of this connectivity (de Pasquale et al. 2011; Baker et al. 2012). Inherently greater temporal resolution in EEG and MEG make these methods naturally more suited for investigations of dynamics, though the unambiguous source localization and full coverage of deep cortex and subcortical nuclei afforded by fMRI advocates the complementary use of both techniques in studies of FC and connectivity changes. We believe further investigations into the consistency and specificity between different fMRI-based FC states and various electrophysiological features (including microstates and band-limited power fluctuations) will be fruitful and important lines of work to elucidate spatiotemporal dynamics associated with spontaneous cognition and behavioral transitions.

Limitations and Future Directions

The results presented here must be considered in the context of several experimental and methodological limitations. First and foremost, our ability to make inferences from FC dynamics and states is limited. Owing to the unconstrained nature of the resting-state scan, we have few tools with which to interrogate changes in FC; the functional roles of dynamics and their relationships to subjects’ cognitive states (if any) remain unknown. As demonstrated by (Hutchison et al. 2012), fluctuations in FC are readily observed in the anesthetized brain and may simply represent the repertoire of spontaneous patterns that one would expect based on underlying anatomical connectivity (Deco et al. 2011). Critical tests for a functional role of connectivity dynamics would be their influence on perception, cognition, or behavior, as demonstrated previously for a number of electrophysiological and hemodynamic signatures (e.g., Arieli et al. 1996; Boly et al. 2007; Eichele et al. 2008). Such evidence is emerging, with demonstrations that trial-to-trial variations in a large-scale connectivity affect response speed in a psychomotor vigilance task (Thompson et al. 2012), and that FC dynamics between particular regions may be related to fluctuations in autonomic system activity and overall awareness (Chang et al. 2012).

One should also consider that observed dynamics may be driven by time-varying noise (e.g., subject motion and variable respiratory and cardiac rhythms) despite our attempts to minimize these influences via aggressive outlier removal (Power et al. 2012) and a high model order ICA decomposition (Birn et al. 2008; Starck et al. 2010). The observed distinctions between ICN and ART components (see Supplementary Fig. S3) offer some confirmation that the observed effects are not solely due to nonspecific physiological changes. Future work should consider multimodal approaches such as concurrent EEG-fMRI, to determine electrophysiological differences between FC states, as well as subtle experimental manipulations of subjects’ internal states, to permit the mapping (and decoding) of cognitive states from connectivity data (e.g., Richiardi et al. 2011; Gonzalez-Castillo et al. 2012; Shirer et al. 2012).

A second experimental limitation of this study is the quantity of data available for each participant. Though we were able to compile data from a very large number of subjects, which benefits several aspects of the analyses, each subject was only scanned for approximately 5 min, precluding robust characterization of connectivity dynamics and state transitions at the level of the individual. Longer scanning times (ideally tens of minutes) will improve estimates of FC variability and permit patterns of connectivity to reoccur several times, which may be critical for future investigations that examine relationships between FC dynamics and behavioral variability within and between subjects.

Methodological limitations relate to the assumptions and constraints inherent to specific analyses. 1) In using spatial ICA, we impose a model of spatial stationarity on the data and assume that the structure of ICNs remains relatively constant over time. This may appear at odds with the work of Kiviniemi et al. (2011) who demonstrate substantial spatial dynamics when employing sliding time window ICA. However, these authors used low model order ICA decompositions (on average 15 components) to iteratively identify a single, large DM network undergoing structural changes over time. In this work, we capture prototypical DM regions in a number of components that show large changes in FC with each other, providing excellent agreement with (Kiviniemi et al. 2011). Future work could consider higher model order ICA or clustering to delineate regions that are more functionally homogenous (Craddock et al. 2012); however, as the number of nodes grows, it will become increasingly difficult to estimate the covariance matrix from short windows. Faster imaging methods, such as MR-encephalography (Zahneisen et al. 2011) or multiband excitation (Moeller et al. 2010), may help to overcome such limitations by increasing the number of samples from which to compute covariance and potentially enhancing signal quality by removing confounding physiological signals. It is also possible that increased sampling might reveal additional, faster dynamics in FC, though current data suggest intrinsic fluctuations are primarily low-frequency phenomena, adequately sampled with typical TRs.

Notably, the assumption of spatial stationarity made here is in contrast to the assumption of temporal stationarity made by Smith et al. in their “temporal” ICA analysis of resting-state data (Smith et al. 2012). While both studies pursue the same general goal of a more detailed analysis of intrinsic connectivity, they do so with different aims and approaches. We focus on the variability of temporal interactions between subnetworks. In doing so, we assume that these subnetworks are spatially fixed over time. Smith et al. focus on the spatial complexity and overlap between networks, and implicitly assume that temporal relationships are consistent over time. Thus, their approach is less suited for explicit investigations of temporal dynamics. Importantly, the results of the respective approaches need not be contradictory; in fact, both approaches highlight subdivisions of the DM system and their complex interactions with other brain regions (e.g., see Fig. 2 in Smith et al. 2012). A more detailed discussion of spatial and temporal ICA and other possible models for decomposing intrinsic activity can be found in our recent commentary (Calhoun et al. 2012).

Additional methodological limitations follow: 2) we characterized FC as the covariance between ICN TCs, rather than use metrics based on higher order statistics, such as mutual information, or lag-insensitive measures such as cross-correlation or coherence. While the use of covariance restricts the detection of nonlinear dependencies and the resolution phase of and spectral relationships, it is preferred for its straightforward interpretation and tractability (e.g., the use of coherence as in Chang and Glover 2010 would estimate connectivity in terms of magnitude and phase as a function of time and frequency for “each pair” of regions; applying this method to the 50 regions studied here would present immense challenges with regard to data analysis and visualization). Future work utilizing both real and simulated data should explore the suitability of different connectivity metrics as applied to studies of dynamics in large-scale networks. 3) Dynamics were estimated using a sliding window size of 22 TRs (44 s). In initial work, we varied window size from 30 s to 2min and found relatively little impact on dynamics beyond the expected result that larger windows reduce variability (Chang and Glover 2010; Hutchison et al. 2012). Comparisons between window sizes suggested 44 s provided a good tradeoff between the ability to resolve dynamics and the quality of covariance matrix estimation, in agreement with demonstrations that cognitive states may be correctly identified from covariance matrices estimated on as little as 30−60 s of data (Shirer et al. 2012), and that topological assessments of brain networks begin to stabilize at window lengths of roughly 30 s (Jones et al. 2012). 4) To identify patterns of FC, we used k-means clustering. Though k-means is an efficient and robust partitioning algorithm, it has several known limitations, in particular difficulty separating clusters of different sizes and densities as well as a high susceptibility to outliers. Alternative clustering models (e.g., density-based or fuzzy-clustering methods) that are not subject to the same limitations may be better suited to FC distributions. A variety of other approaches to identify FC states are also possible, including using topological descriptions of brain connectivity as features (e.g., modularity or community membership [Bassett et al. 2011; Jones et al. 2012; Kinnison et al. 2012]) rather than the connectivity values themselves, or using formal models for detecting change points in connectivity, as recently introduced by Cribben et al. (2012), rather than clustering. We are hopeful that further work will develop and improve methods for the identification of FC states and state transitions.

Funding

This research was supported by NIH 1R01-EB006841, NIH 1R01-EB005846, NIH 2R01-EB000840, NIH 1 P20 RR021938-01, and DOE DEFG02-08ER64581 to V.D.C. and a Nevronor grant from the Norwegian Research Council to T.E. E.A.A. is funded by the K.G. Jebsen Center for Research on Neuropsychiatric Disorders.

Notes

We thank Martin Havlicek, Arvind Caprihan, and J. Bruce Morton for helpful discussions and comments that improved this work. Conflict of Interest: The authors declare no competing financial interests in relation to the work presented.