# Time-resolved resting-state brain networks

^{a}Melbourne Neuropsychiatry Centre, The University of Melbourne and Melbourne Health, Melbourne, VIC 3010, Australia;^{b}Melbourne School of Engineering, The University of Melbourne, Melbourne, VIC 3010, Australia;^{c}Monash Clinical and Imaging Neuroscience, School of Psychological Sciences and Monash Biomedical Imaging, Monash University, Melbourne, VIC 3168, Australia;^{d}Queensland Brain Institute, The University of Queensland, Brisbane, QLD 4072, Australia; and^{e}QIMR Berghofer Medical Research Institute,^{f}The Royal Brisbane and Women's Hospital, and^{g}Metro North Mental Health, Brisbane, QLD 4029, Australia

See allHide authors and affiliations

Edited by Michael S. Gazzaniga, University of California, Santa Barbara, CA, and approved May 19, 2014 (received for review January 13, 2014)

## Significance

Large-scale organizational properties of brain networks mapped with functional magnetic resonance imaging have been studied in a time-averaged sense. This is an oversimplification. We demonstrate that brain activity between multiple pairs of spatially distributed regions spontaneously fluctuates in and out of correlation over time in a globally coordinated manner, giving rise to sporadic intervals during which information can be efficiently exchanged between neuronal populations. We argue that dynamic fluctuations in the brain’s organizational properties may minimize metabolic requirements while maintaining the brain in a responsive state.

## Abstract

Neuronal dynamics display a complex spatiotemporal structure involving the precise, context-dependent coordination of activation patterns across a large number of spatially distributed regions. Functional magnetic resonance imaging (fMRI) has played a central role in demonstrating the nontrivial spatial and topological structure of these interactions, but thus far has been limited in its capacity to study their temporal evolution. Here, using high-resolution resting-state fMRI data obtained from the Human Connectome Project, we mapped time-resolved functional connectivity across the entire brain at a subsecond resolution with the aim of understanding how nonstationary fluctuations in pairwise interactions between regions relate to large-scale topological properties of the human brain. We report evidence for a consistent set of functional connections that show pronounced fluctuations in their strength over time. The most dynamic connections are intermodular, linking elements from topologically separable subsystems, and localize to known hubs of default mode and fronto-parietal systems. We found that spatially distributed regions spontaneously increased, for brief intervals, the efficiency with which they can transfer information, producing temporary, globally efficient network states. Our findings suggest that brain dynamics give rise to variations in complex network properties over time, possibly achieving a balance between efficient information-processing and metabolic expenditure.

The coordination of brain activity between disparate neural populations is a dynamic and context-dependent process (1⇓–3). Although dynamic patterns of neural synchronization may be evident in time-dependent measures of functional connectivity (4, 5), the temporal stability of high-level topological properties is unknown. The topology of large-scale cortical activity—such as its efficient network layout (6), community structure (7), network hubs (8), rich-club organization (9, 10), and small worldness (11, 12)—may reflect fundamental aspects of cortical computation. Temporal fluctuations in these graph-theoretic measures may hence speak to adaptive properties of neuronal information processing.

With international connectome mapping consortia such as the Human Connectome Project (HCP) (13) and the developing Human Connectome Project in full swing, resting-state functional magnetic resonance imaging (rsfMRI) data of unprecedented temporal resolution are now available to map the time-resolved properties of functional brain networks. Imaging the brain at rest reveals spontaneous low-frequency fluctuations in brain activity that are temporally correlated between functionally related regions (14⇓⇓–17). Interregional correlations are referred to as functional connections, and they collectively form a complex network (18).

Functional brain networks are typically mapped in a time-averaged sense, based on the assumption that functional connections remain relatively static (stationary) in the resting brain. However, recent investigations have furnished compelling evidence challenging the “static” conceptualization of resting-state functional connectivity (5). In particular, the application of time-resolved methodologies for analyzing time series data has consistently revealed fluctuations in resting-state functional connectivity at timescales ranging from tens of seconds to a few minutes (19⇓⇓⇓⇓–24). Furthermore, the modular organization of functional brain networks appears to be time-dependent in the resting state (25, 26) and modulated by learning (27) and cognitive effort (28, 29). It is therefore apparent that reducing fluctuations in functional connectivity to time averages has led to a very useful but static and possibly oversimplified characterization of the brain’s functional networks. For example, connections that toggle between correlated and anticorrelated states are reduced to zero in a time-averaged sense, assuming equal dwell times in each state.

Conventionally, rsfMRI data are sampled at a resolution of 2 s or slower. Using multiband accelerated echo planar imaging, the HCP has acquired high-quality rsfMRI data at a subsecond resolution (30). This order of magnitude improvement in temporal resolution is highly advantageous to the feasibility of time-resolved functional connectomics. Faster sampling rates enable a richer temporal characterization of resting-state fluctuations, denser sampling of physiological confounds, and greater degrees of freedom (30).

Using a sliding-window approach applied to HCP rsfMRI data, we mapped the evolution of functional brain networks over a continuous 15-min interval at a subsecond resolution. For each of 10 individuals, this yielded a time series of correlation matrices (regions × regions × time), where matrix elements quantified the functional connectivity at a given time instant between cortical and subcortical regions comprising established brain parcellation atlases. We developed a statistic to test the time-resolved connectome for evidence of nonstationary temporal dynamics and applied it to the 10 individuals as well as a replication data set and simulated rsfMRI data.

Our main aim was to investigate the consequences of nonstationary fluctuations on the topological organization of functional brain networks. We hypothesized that dynamic behavior is coordinated across the brain so that transitions between distinct states are marked by reorganization of the brain’s functional topology. Evidence for this hypothesis is provided by the coordinated fluctuations in network measures, such as hub centrality (31), that have been observed in simulated rsfMRI data (32, 33).

## Results

Time-resolved functional brain connectivity was mapped using a sliding-window approach applied to high-resolution rsfMRI data acquired in 10 healthy, young adults participating in the HCP (*Materials and Methods*). Connectivity was estimated using pairwise linear correlation in regionally averaged rsfMRI time series data falling within fixed-length time windows (19, 21, 22, 26). We used a tapered window of length 60 s (83 time points per window). Sliding the window in time yielded a continuous series of snapshots characterizing the evolution of each individual’s functional brain network at a temporal resolution of 720 ms over a 15-min interval.

We developed a statistic to test for time-varying (nonstationary) connectivity. To estimate the statistic’s distribution under the null hypothesis of stationarity, 250 null data sets were generated using stable vector autoregressive (VAR) models (34) that approximately preserved the power and cross-spectrum of the actual rsfMRI data (*Materials and Methods* and *SI Appendix*, Fig. S1). The null hypothesis was rejected pairwise when the time-resolved correlation coefficients fluctuated on a timescale longer than the window length and/or between more extreme correlation levels than expected by chance [*P* < 0.01, familywise error rate (FWER) corrected across all connections]. This enabled partitioning of each individual’s functional brain network into dynamic (nonstationary) and static (stationary) connections.

We first tested all 6,670 connections defined by the 116-region Automated Anatomical Labeling (AAL) (35) atlas for evidence of time-varying connectivity. Across the 10 individuals, the number of connections for which the null hypothesis of stationarity was rejected varied between 118 and 570 [mean: 293 (∼4%); SE: 54]. The top-100 most dynamic connections for each individual, according to test statistic magnitude, could therefore be declared nonstationary for *all* individuals.

An index of consistency was then used to establish whether the top-100 most dynamic connections were consistent (i.e., overlapped) across the 10 individuals. For each individual, a binary graph comprising only the top-100 most dynamic connections was constructed. The degree of each region in these 100-edge graphs was determined and then summed across the 10 individuals to give a region-specific index of consistency/overlap. Regionally sorted from lowest-to-highest overlap, this index of consistency is shown in Fig. 1*A* for the actual data (blue line) and for the 250 VAR null data sets described above (black lines). The 19 regions residing to the right of the *P* = 0.01 cutoff (vertical red line) were therefore associated with dynamic behavior more consistently than what would be expected by chance alone. Cortical renderings of consistency are shown in Fig. 1*B*. The 19 regions are as follows: angular gyrus (l/r); supramarginal gyrus (l/r); rectus gyrus (l/r); medial orbitofrontal cortex (l/r); inferior parietal lobule (l/r); inferior frontal operculum (l/r); middle frontal gyrus (l/r); amygdala (r); superior temporal pole (r); olfactory cortex (l); and postcentral gyrus (l) and precentral gyrus (l) with left (l), right (r), and both (l/r) hemispheres. Several of these regions are frontal and parietal association areas that are known hubs of the structural human connectome (36), comprise some parietal areas of the default mode network (37), and include rich-club nodes (9, 10). The cerebellum, vermis, and temporal regions were consistently the least dynamic.

Time-resolved correlation coefficients for the top-100 most dynamic connections are shown in Fig. 2*A* for two representative individuals and a VAR null data set. It can be seen that fluctuations in the time-resolved correlation coefficients are synchronized across the brain, occurring at distinct moments in time, where multiple connections transition en masse between different correlation levels. This is in contrast to the simplest null hypothesis where transitions between correlation levels occur independently among connections.

To test this null hypothesis, the total number of connections (from the top-100 most dynamic connections) that transitioned at each time point was enumerated. This yielded a time series of counts for each individual ranging between 0 and 100, which we referred to as the “transition count.” A connection underwent a “transition” at times when its series of time-resolved correlation coefficients crossed their median value. The median correlation value was used here to represent the crossing point between two correlation levels. Using this simple two-state characterization, we tested the null hypothesis of uniformly distributed transitions across time. Null data sets were generated by randomly redistributing groups of related transition events in time, thus preserving synchrony owing to inherent properties of correlation networks (38), but randomizing synchrony resulting from specific aspects of dynamic connectivity (*SI Appendix*). The null hypothesis was rejected within each individual because at least one time point was always found where the transition count exceeded the minimum cutoff, satisfying a FWER of 0.01. Transition counts (blue lines) and the 0.01 FWER cutoff value (horizontal red lines) are provided in Fig. 2*B* for the same two individuals. The null hypothesis can be rejected at time points where the transition count exceeds the cutoff value. This result verifies that nonstationary fluctuations are synchronized across functional brain networks in such a way that multiple connections transition en masse at distinct moments in time between different correlation levels. Finally, we tested for power law scaling in the distribution of transition counts to establish a possible correspondence with the neuronal avalanche phenomenon (39). However, exponentials and stretched exponentials (40) provided better fits than a power law.

The results displayed in Figs. 1 and 2 were replicated using two different parcellation atlases: Craddock-200 (41) and random-90 (42). Random-90 comprises parcels that are equal in volume. The null hypothesis of uniformly distributed transitions was also rejected when using the top-500 or top-1,000 most dynamic connections.

To further verify this effect, a principal component analysis was performed on the top-100 most dynamic connections to identify any prevailing temporal patterns. For each individual, at least one principal component explaining at least 20% of the variance was present. In comparison, the probability of a principal component explaining at least 20% of the variance in the VAR null data was not significant for all individuals (*P* < 0.01).

### Dynamics of Efficiency.

We next performed a time-resolved analysis of network efficiency to test whether the dynamic fluctuations that we identified at the pairwise level extend to a complex network property.

Time-resolved connectivity estimates were thresholded to yield a continuous series of weighted networks, each with a fixed connection density of 20%. Regional efficiencies (43) were calculated for these weighted networks, resulting in a time series of efficiency values for each region (*SI Appendix*). These time-resolved regional efficiencies are shown in Fig. 3*A* for two representative individuals and a VAR null data set. It can be observed that several regions show globally coordinated transitions between low- and high-efficiency states, whereas others show less pronounced fluctuations. No such patterns were observed in the null data. To quantify this observation, the variance in regional efficiency over time was computed in the actual data and in the VAR null data. It was found that the largest variance in efficiency across all 250 null data sets never exceeded the smallest variance in the actual data (*SI Appendix*, Fig. S2), confirming that all regions displayed fluctuations in efficiency that were more variable than the stationary null data. This effect was replicated using several alternative connection densities (*SI Appendix*, Fig. S3) and also observed at the level of global efficiencies (*SI Appendix*, Fig. S4). It can also be observed that transitions from low-to-high efficiencies are sudden, whereas high-to-low transitions are gradual. To quantify this observation, the skewness in the forward difference of time-resolved efficiency was computed in the actual data and in the VAR null data. Forward differences were significantly positively skewed for all individuals (*P* < 0.01; skewness range: 0.17–1.3; *SI Appendix*, Fig. S5).

Cortical renderings of time-resolved regional efficiencies were compiled into movies for two of the individuals featured in Fig. 3*A* (Movies S1 and S2 and *SI Appendix*). Snapshots from these movies are shown in Fig. 3*B* at time points residing within an interval spanning a low-efficiency state and at another soon after the transition to a high-efficiency state.

We next tested whether high-efficiency states were more “costly” in terms of the anatomical distance between interconnected regions. Interregional distances were computed for pairs of suprathreshold connections, with the distance between pairs of regions approximated by the Euclidean distance between regional centers of mass. Mean connection length in high-efficiency states was then estimated by averaging these interregional distances over all time points where the global efficiency value exceeded its median. This was repeated for all individuals to yield an estimate of the mean connection length in low- and high-efficiency states. When pooling data across the original and replication data sets (see below), connections were found to be significantly longer in the high-efficiency states (*P* < 0.01; low: 67.1 ± 2.1 mm; high: 73 ± 1.9 mm).

Network efficiency is thought to reflect a network’s capacity for information transfer (6). High-efficiency networks may be energy demanding (43, 44), as suggested by increased cerebral blood flow to the strong, long-range functional connections facilitating integration across disparate network elements (45). Sporadic emergence of metabolically costly, high-efficiency states lasting for brief intervals may have evolved to minimize energy demands while maintaining the connectome in a globally integrated, responsive state.

### Dynamics of Modular Organization.

Modules refer to communities of regions that are more strongly interconnected with each other than with regions outside their community. Time-resolved modularity analysis (46) suggests that modular organization of functional brain networks is dynamic (25) and shaped by learning-dependent plasticity (27). We performed a conventional, time-averaged modular decomposition of functional brain networks (7) to understand the spatiotemporal dynamics of inter- and intramodular connections.

An algorithm to determine a modular decomposition representing consensus among the 10 individuals yielded four modules (*SI Appendix*). Modules rendered onto the cortical surface are shown in Fig. 4*A* and broadly overlap established resting-state networks, although the correspondence is not precise because modules were constrained to conform to coarse AAL regional boundaries. The values of the statistic used to test whether a connection was static/dynamic are shown in Fig. 4*B* in matrix form and averaged across the 10 individuals. The four blocks positioned along the matrix diagonal correspond to modules and encapsulate intramodular connections. The visual, default mode and somatomotor modules can be seen to comprise disproportionately many static connections. This suggests that static connections predominate between regions composing the same module, whereas intermodular connections are dynamic. To test this observation in each individual, the proportion of the top-100 most dynamic connections that were intermodular was divided by the overall proportion of *all* intermodular connections. This ratio was significantly greater than unity (ratio: 1.13 ± 0.02; *P* = 0.0001; 83 ± 1.5% of top-100 intermodular, 69% of all connections intermodular).

Dynamic connections were significantly more prevalent between spatially distant regions. The mean interregional distance for *all* connections was 77 mm, whereas the top-100 most dynamic connections were on average 83 ± 1.7 mm (*P* = 0.004). For all individuals, the time-averaged correlation coefficients were significantly anticorrelated with the statistic developed to test for nonstationarity (*r* < −0.4, *P* < 0.01; *SI Appendix*, Fig. S6). Hence, the most dynamic connections were also the connections with the weakest time-averaged correlations. Intermodular connections likely comprised the most dynamic connections because dynamic connections are near zero in the time averages used by the modular decomposition algorithm.

### Noise Confounds.

Time-resolved analysis of rsfMRI data are particularly susceptible to noise confounds. Confounds include scanner drift, head motion (47), and physiological noise due to variation in respiratory depth/rate (34) and cardiac rate (48). To exclude these confounds as possible causes of the results, for each individual, the principal component explaining the most temporal variation in the top-100 most dynamic connections was regressed against estimates of physiological noise as well as estimates of instantaneous head motion, namely, displacement and rotation in *x*, *y*, and *z* directions and all associated first-order derivatives. After controlling the false discovery rate at a relatively liberal threshold of 10% across all regressors, no measure of head motion or physiological noise was a significant predictor of connectivity fluctuations for any individual. Noise confounds and principal components are shown in *SI Appendix*, Fig. S7.

### Simulated rsfMRI Data.

To further exclude noise confounds as a potential explanation for the dynamic behavior that we identified, our findings were replicated using simulated rsfMRI data, which was necessarily free of any head motion, scanner drift, and physiological noise. Neuronal population dynamics were simulated for 47 neural masses interconnected according to the axonal connectivity of the macaque neocortex (31). The Balloon–Windkessel model was then applied to the simulated neuronal dynamics to generate realistic rsfMRI data matched in length and temporal resolution to the HCP data (*SI Appendix*). Using this noise-free, simulated rsfMRI data, we replicated our findings, with regional efficiencies displaying coordinated fluctuations akin to those seen in the HCP data (Fig. 3 and *SI Appendix*, Fig. S8).

### Replication Data Set.

Our main findings were replicated in an independent data set comprising an unrelated group of 10 individuals (*Material and Methods*). The replication data set was acquired with a different phase encoding and preprocessed using an alternative method for noise removal (30). The null hypothesis of uniformly distributed transitions across time was rejected within all individuals at 0.05 FWER and within 9 of the 10 individuals at 0.01 FWER. Highly synchronized fluctuations in regional efficiency were once again evident in all individuals, as shown for three individuals in *SI Appendix*, Fig. S9.

## Discussion

By mapping time-resolved functional brain networks at a subsecond resolution, we here report evidence for dynamic (nonstationary) behavior in the brain’s resting state from the scale of simple pairwise temporal correlations to a complex network property. We found that dynamic behavior was coordinated across the cortex, with hemodynamic activity between multiple pairs of spatially distributed regions spontaneously transitioning in and out of correlation over time in a globally coordinated manner.

Our results accord with electroencephalographic and theoretical investigations of dynamic connectivity (1⇓–3) that suggest a “natural partitioning” (1) of functional dynamics into synchronous epochs. Alternating patterns of correlation and anticorrelation may constitute fundamental dynamics of information processing by allowing the formation and dissolution of dynamic cell assemblies (3, 4). Moreover, intermittent epochs of global synchronization may enable otherwise segregated network elements access to a cognitive global workspace, which may be necessary for effortful processing (49). Transient exploration of this workspace may allow the brain to efficiently balance segregated and integrated neural dynamics.

Epochs at which more functional connections transitioned (i.e., crossed their median value) than expected by chance alone were sporadically distributed in time and found in all 10 individuals as well as in a replication data set and simulated rsfMRI data. We hypothesize that these “transition epochs” mark changeover points between distinct metastates (33, 50). Our observation of concomitant changes in a complex network property—network efficiency—suggests that dynamic fluctuations at the pairwise level might be coordinated in such a way as to achieve a topological “objective.”

Our time-resolved analysis of network efficiency revealed that multiple spatially distributed regions simultaneously increased, for brief intervals, their topological efficiency and, by inference, their capacity to transfer information. However, these intervals of high efficiency are supported by long anatomical connections and thus likely carry an extra metabolic cost (44). We argue that intermittent periods of high efficiency may hence be a dynamic strategy that has evolved to minimize metabolic requirements, analogous to intermittent search strategies that constitute an optimal solution in settings as diverse as food foraging in animals and eye movements in humans (51). Some of the dynamic regions that we identified have been described as “transmodal” (e.g., default mode parietal areas), owing to their functional associations with multiple intrinsic connectivity networks (ICNs) (52). On the basis that some transmodal regions are also highly dynamic, we suggest that their multiple functional associations may be realized through a dynamic process of time-division multiplexing, where the region is connected with specific ICNs for a fraction of time. A time-averaged analysis may thus reveal the “echoes” (52) of multiple ICNs at transmodal regions.

In this study, we also investigated the role of dynamic behavior with reference to a conventional time-averaged modular decomposition of functional brain networks. We found that the most dynamic connections were intermodular and localized to known hubs of default mode and fronto-parietal systems (36), suggesting that these hubs were connectors linking multiple modules with one another (53). Our results suggest that time-averaged modular decompositions may be explained by differences in the topographic layout of dynamic and static connections. In particular, we found that dynamic connections were typically connections with time-averaged correlations near zero and were thus more likely to straddle two modules. Intramodular connections composing the orbitofrontal-limbic module were relatively dynamic compared with the other three modules identified. This “transitional” module (29) may therefore be more flexible than the others, possibly supporting transient psychological states (19) and complementing frontoparietal regions in supporting adaptive, context-dependent control (28, 54).

### Methodological Considerations.

First, we used traditional volume-based parcellations of cortical and subcortical regions, as defined by established volumetric atlases. Recent investigations have demonstrated the benefit of surface-based parcellation (30, 55), which may reduce heterogeneity in region size. Larger regions are more likely to encapsulate multiple temporally independent modes (24) that might cancel each other out when averaged. The dynamic behavior in regionally averaged activity may therefore be influenced by region size, although we found no evidence of a linear relation between volume and our index of consistency. Second, we used what is effectively the Pearson (full) correlation coefficient to estimate functional connectivity within each time window. A criticism of full correlation is its sensitivity to indirect functional relations between pairs of regions that are mediated by a third region (38). Third, methodological options for preprocessing rsfMRI data are many and varied. We used preprocessing options recommended by the HCP and replicated our main findings using a different noise removal method applied to an independent data set as well as simulated rsfMRI data that were necessarily free from any noise confounds. We replicated our findings using the “scrubbing” procedure (47) to correct for head motion (*SI Appendix*, Fig. S10), noting that scrubbing should be used cautiously because it introduces variation in the degrees of freedom per time window.

### Conclusion.

Time-averaged characterizations of functional brain networks are inherently static and as such reduce the rich temporal dynamics of the resting brain to temporal averages. This represents an oversimplification. With an abundance of high temporal resolution rsfMRI data to be released over the next few years, time-resolved analysis of functional brain networks and their topological organization will become a feasible and likely widespread analysis. Here we have shown that dynamic fluctuations in functional connectivity at the pairwise level appear to be coordinated across the brain so as to realize globally coordinated variations in network efficiency over time, which might represent a balance between optimizing information processing and minimizing metabolic expenditure.

## Materials and Methods

### Functional MRI Data.

Minimally preprocessed rsfMRI data for 10 healthy, unrelated adults (age: 22–35, 4 males) were obtained from the Human Connectome Project (13). Data were obtained under the Q2 Data Release and comprised 1,200 frames of multiband, gradient-echo planar imaging acquired during a period of 14 min and 33 s with the following parameters: relaxation time, 720 ms; echo time, 33.1 ms; flip angle, 52°; field of view, 280 × 180 mm; matrix, 140 × 90; and voxel dimensions, 2 mm isotropic. Individuals were fixated on a projected bright crosshair on a dark background during data acquisition. Only one of the four runs acquired for each individual was analyzed in this study (left-right encoded, second session).

### Replication Data Set.

An independent data set comprising 10 healthy, unrelated adults was obtained under the Q3 Data Release (age: 26–35, 5 males). Right-left encoding runs acquired in the first session were analyzed. Data were already preprocessed using a metaclassification approach applied to independent components (independent component analysis-based X-noisifier) (30). *SI Appendix*, Table S1, provides details of the individuals composing the main and replication data sets.

### Time-Resolved Functional Connectivity.

Time-resolved functional connectivity was estimated between pairs of regional time series using a tapered sliding-window approach. Tapering provides better suppression of spurious correlations and may reduce sensitivity to outliers. An exponentially tapered window (56) spanning *N* time points was defined by the weight vector *w*_{τ} = *w*_{0}*e*^{(τ} ^{−} ^{N)/θ}, τ = 1, …, *N*, θ > 0, and *w*_{0} = (1 − *e*^{−1/θ})/(1 − *e*^{−N/θ}). The weighted Pearson product-moment correlation between region *i* and region *j* at time *t* ≥ *N* was then computed as *SI Appendix*. The window length was set to 60 s and the exponent, θ, was set to a third of the window length (56). Statistical analysis was performed on the time-resolved correlation coefficients

### Nonstationarity Test Statistic.

A univariate test statistic was developed to measure the extent of time-varying (nonstationary) fluctuations in the time-resolved correlation coefficients for each pair of regions. For a given pair of regions *i* and *j*, median crossing points were defined as the set of solutions *J* denotes the total number of median crossing points. A pair of consecutive crossing points (*t*_{n}, *t*_{n+1}) defined an excursion from the median value. The longer and larger the excursions were from the median value, the greater the evidence for nonstationary behavior. A statistic that increased as a function of excursion length and/or height was therefore devised. The length of the *n*th excursion was the time interval *l*_{n} = *t*_{n}_{+}_{1} − *t*_{n}. The height of the *n*th excursion was *SI Appendix*, Fig. S11.

### Vector Autoregressive Null Model.

Stable (stationary) VAR models were fitted to the regional time series using maximum-likelihood estimation. VAR model responses were then simulated to generate surrogate regional time series, satisfying the null hypothesis of a linearly correlated, stationary, multivariate stochastic process. This was repeated to generate 250 independent null data sets. The surrogate regional time series comprising each null data set were then processed in the same way as the actual data. This enabled empirical estimation of the null distribution of the test statistic and graph measures evaluated in this study. It was computationally infeasible to fit a single multidimensional VAR model with a covariance structure with dimensions equal to the number of regions. Two-dimensional VAR models were therefore independently fitted to each pair of regional time series, implying that the model response for a given region was conditional for the pair of regions under consideration. See *SI Appendix* for details of simulation of VAR model responses. This study was approved by the QIMR Human Research Ethics Committee (P1331).

## Acknowledgments

Data were provided by the Human Connectome Project, WU–Minn Consortium (1U54MH091657; Principal Investigators: David Van Essen and Kamil Ugurbil) funded by the 16 National Institutes of Health (NIH) institutes and centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. We acknowledge support provided by the Australian National Health and Medical Research Council [Career Development Fellowship GNT1047648 (to A.Z.)], Australian Research Council, Monash Larkins Award (to A.F.), Queensland Health Fellowship (to M.B.), and the James S. McDonnell Foundation (M.B. and L.L.G.).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: azalesky{at}unimelb.edu.au.

Author contributions: A.Z., A.F., L.C., L.L.G., and M.B. designed research; A.Z., A.F., L.C., L.L.G., and M.B. performed research; A.Z. contributed new reagents/analytic tools; A.Z., A.F., L.C., L.L.G., and M.B. analyzed data; and A.Z., A.F., L.C., L.L.G., and M.B. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1400181111/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- van den Heuvel MP,
- Sporns O

- ↵
- van den Heuvel MP,
- Kahn RS,
- Goñi J,
- Sporns O

- ↵
- Bassett DS,
- Bullmore E

- ↵
- ↵
- ↵
- Biswal BB,
- et al.

- ↵
- Fox MD,
- et al.

- ↵
- ↵
- Greicius MD,
- Supekar K,
- Menon V,
- Dougherty RF

- ↵
- Sporns O

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Smith SM,
- et al.

- ↵
- Allen EA,
- et al.

- ↵
- ↵
- Bassett DS,
- et al.

- ↵
- Kitzbichler MG,
- Henson RN,
- Smith ML,
- Nathan PJ,
- Bullmore ET

- ↵
- Fornito A,
- Harrison BJ,
- Zalesky A,
- Simons JS

- ↵
- ↵
- Honey CJ,
- Kötter R,
- Breakspear M,
- Sporns O

- ↵
- Deco G,
- Jirsa VK

- ↵
- ↵
- ↵
- ↵
- ↵
- Raichle ME,
- et al.

- ↵
- ↵
- Beggs JM,
- Plenz D

- ↵
- Freyer F,
- et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- Liang X,
- Zou Q,
- He Y,
- Yang Y

- ↵
- ↵
- ↵
- ↵
- Dehaene S,
- Kerszberg M,
- Changeux JP

- ↵
- Van de Ville D,
- Britz J,
- Michel CM

- ↵
- ↵
- Braga RM,
- Sharp DJ,
- Leeson C,
- Wise RJ,
- Leech R

- ↵
- ↵
- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Biological Sciences
- Neuroscience