Ghost Attractors in Spontaneous Brain Activity: Recurrent Excursions Into Functionally-Relevant BOLD Phase-Locking States

Functionally relevant network patterns form transiently in brain activity during rest, where a given subset of brain areas exhibits temporally synchronized BOLD signals. To adequately assess the biophysical mechanisms governing intrinsic brain activity, a detailed characterization of the dynamical features of functional networks is needed from the experimental side to constrain theoretical models. In this work, we use an open-source fMRI dataset from 100 healthy participants from the Human Connectome Project and analyze whole-brain activity using Leading Eigenvector Dynamics Analysis (LEiDA), which serves to characterize brain activity at each time point by its whole-brain BOLD phase-locking pattern. Clustering these BOLD phase-locking patterns into a set of k states, we demonstrate that the cluster centroids closely overlap with reference functional subsystems. Borrowing tools from dynamical systems theory, we characterize spontaneous brain activity in the form of trajectories within the state space, calculating the Fractional Occupancy and the Dwell Times of each state, as well as the Transition Probabilities between states. Finally, we demonstrate that within-subject reliability is maximized when including the high frequency components of the BOLD signal (>0.1 Hz), indicating the existence of individual fingerprints in dynamical patterns evolving at least as fast as the temporal resolution of acquisition (here TR = 0.72 s). Our results reinforce the mechanistic scenario that resting-state networks are the expression of erratic excursions from a baseline synchronous steady state into weakly-stable partially-synchronized states – which we term ghost attractors. To better understand the rules governing the transitions between ghost attractors, we use methods from dynamical systems theory, giving insights into high-order mechanisms underlying brain function.

Introduction

For healthy human cognition, the brain needs to engage in functionally meaningful activity through an integration of information incoming from various segregated brain areas (Tononi and Edelman, 1998; Sporns et al., 2000). At rest, brain activity has been shown to reveal the spontaneous activation of meaningful functional subsystems, sharing spatial features with networks of brain areas typically activated during task (Beckmann et al., 2005; Fox et al., 2005; Damoiseaux et al., 2006). These spatially activated coalitions of brain regions, dubbed resting-state networks (RSNs), have been remarkably consistent across neuroimaging studies and utilized in describing functional changes in disruptions to the healthy brain functioning (Greicius, 2008; Fox and Greicius, 2010; van den Heuvel and Hulshoff Pol, 2010; Vargas et al., 2013; Kaiser et al., 2015). However, while RSNs represent spatially meaningful information, in order to further investigate the generative mechanisms of RSNs and their functional role, it is important to further characterize their behavior in the temporal domain (Preti et al., 2016; Cabral et al., 2017b).

Indeed, recent advances have focused on how these spatially coherent functional patterns can explain the complex brain dynamics evolving in time (Chang and Glover, 2010; Hutchison et al., 2013; Allen et al., 2014). However, the most appropriate way to characterize network dynamics at the whole brain level is still unclear. The most common approach to dynamic functional connectivity (dFC) has been the sliding-window method, which describes statistical relationship between brain regions in successive intervals of time and generates recurrent states of functional connectivity using unsupervised learning (Hutchison et al., 2013; Allen et al., 2014; Calhoun et al., 2014). However, the choice of the “window” size introduces limitations which hinders the temporal resolution as well as statistical validation (Hindriks et al., 2016; Preti et al., 2016). To overcome these caveats, recent development has focused on describing single frame functional connectivity [FC(t)] either by considering BOLD co-activations (Karahanoğlu and Van De Ville, 2015; Tagliazucchi et al., 2016) or BOLD phase coherence (Glerean et al., 2012; Cabral et al., 2017b). Framewise co-activation analysis considers the brain regions with BOLD signal above a certain threshold before clustering into distinct FC patterns (Karahanoğlu and Van De Ville, 2015; Tagliazucchi et al., 2016). While it allows for higher temporal resolution, it is still dependent on the choice of the threshold as well as limited to describing simultaneous (in-phase) activations. On the other hand, phase coherence techniques represent the time instances as relative phase relationships between brain regions and thus do not require thresholding and are sensitive to phase-shifted patterns (Glerean et al., 2012; Cabral et al., 2017b).

To overcome issues with high data dimensionality, Cabral and colleagues have proposed to represent the instantaneous relationships between brain regions using the largest magnitude eigenvector of BOLD phases (a 1xN vector for each time point) instead of the NxN phase synchronization matrix (Cabral et al., 2017b). Notably, Leading Eigenvector Dynamic Analysis (LEiDA) has been shown not only to improve clustering performance, but to consistently capture meaningful BOLD phase-locking states (PL-states) that closely overlap with previously-described functional subsystems (Cabral et al., 2017b; Figueroa et al., 2019; Lord et al., 2019). By representing whole-brain activity over time as a succession of discrete PL states, it is possible to quantify the fractional occupancy, the probability of transition as well as the Dwell Time of individual states. Importantly, these measures have shown to be significantly related with cognitive performance (Cabral et al., 2017b), to be altered in clinical populations of patients suffering with major depressive disorder (Figueroa et al., 2019), as well as to describe the network-specific modulation of resting-state activity by the psychoactive compound psilocybin (Lord et al., 2019). As such, LEiDA opens up as a useful tool to quantitatively characterize individual fingerprints in dynamic functional connectivity, reinforcing a mechanistic scenario proposed by theoretical works where RSNs are the expression of a repertoire of BOLD FC configurations emerging from complex non-linear interactions in the whole-brain network (Ghosh et al., 2008; Cabral et al., 2011; Deco and Jirsa, 2012; Deco et al., 2013; Haimovici et al., 2013; Hansen et al., 2015).

Here, we explore this mechanistic hypothesis using the mathematical formalism from dynamical systems theory and Markov chains in order to characterize the spatio-temporal dynamics of spontaneous brain activity in terms of probabilistic trajectories between recurrent BOLD phase-locking patterns. We validate the functional role of the patterns obtained by comparing them with known RSNs. Furthermore, we evaluate the stability of BOLD phase locking states based on their Fractional Occupancy, Dwell Times and Transition probabilities. While previous works have applied LEiDA to condition-specific datasets, with reduced sample sizes, here we make secondary use of a large open source dataset of healthy participants, demonstrating the reliability of the yielded metrics across subjects and consecutive fMRI recording sessions.

Materials and Methods

Data

All data used in this work comes from a publicly available database – the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essene and Kamil Ugurbil; 1U54MH091657) with funding from the sixteen NIH Institutes and Centers supporting the NIH Blueprint for Neuroscience Research; and by the McDonell Center for Systems Neuroscience at Washington University.

Participants

100 unrelated subjects [mean age 29.5 years old, 55% females (Glasser et al., 2013)].

Neuroimaging HCP Acquisition

Each participant underwent four resting-state fMRI sessions lasting 14 min 30 s with a repetition time (TR) of 0.72 s, on a 3-T connectome Skyra scanner (Siemens) – two during the first day and two during the second day. The 2 fMRI sessions acquired on the same day differ only in the oblique axial acquisition phase encoding, one being from Left to Right (LR) and the other from Right to Left (RL). The acquisition and pre-processing of the data is fully described in detail at the HCP website https://www.humanconnectome.org/. Here, we used the fMRI data acquired on the first day of scanning. One subject was excluded because one session was missing. In total, two same-day resting-state fMRI sessions from 99 of the 100 unrelated subjects’ sessions were used for the analysis.

Parcellation

To reduce the dimensionality of the voxel-based data (Voxels × Time), the Anatomic Automatic Labeling (AAL) atlas was used to define N = 90 anatomically distinct cortical and sub-cortical regions covering the whole brain, excluding the cerebellum. Data was reduced to size N × Time, with Time = 1200 TR per session, by averaging the BOLD signals in all the voxels associated to each brain region (Figure 1A).

FIGURE 1

Figure 1. Time-evolving patterns of BOLD Phase Dynamics. (A) BOLD signals from a representative fMRI scan of the HCP dataset averaged over all voxels within each region of interest (ROI). ROIs were defined using an anatomically-based parcellation scheme (AAL) covering the entire brain (here excluding the cerebellum). (B) To illustrate BOLD phase dynamics, we select a representative interval of TRs. At each TR, blue circles represent the brain areas whose BOLD phase projects into the main BOLD phase direction (captured by the leading eigenvector, see (C), and red dots represent the brain areas whose BOLD phase project into the opposite direction of the main BOLD phase orientation. It serves to illustrate that the phase-shifted signals (in red) do not directly correspond to supra-threshold BOLD increases. (C) Phase portraits of the analytic BOLD signal at each TR, where the real and imaginary axis represent the cosines and sines of the Hilbert phase at each TR. (D) Representation of the brain patterns captured by the signs of the leading eigenvector at each TR, illustrating how phase-locking patterns evolve smoothly over several TRs, whereas the corresponding BOLD signals (shown in panel B) exhibit significantly different activation patterns over the same range of TRs. (E) Representation of the instantaneous phase coherence matrices obtained at each TR as the cosine of the phase difference at each instant of time.

Analysis

BOLD Phase Dynamics

To compute the phase relationship between brain regions, for each region n with n = 1…N, a BOLD phase θ(n,t) varying in time t, was calculated via Hilbert transform (Glerean et al., 2012). The analytical signal expresses the regional signal x(t) as x(t) = A(t)*θ(t)) with A and θ representing the time-varying amplitude and phase respectively (Figure 1C). The first and last time points were removed from each time series, to exclude the boundary artifacts induced by the Hilbert transform. Subsequently, for every pair of brain regions n and m at time t the phase coherence matrix dPC is calculated as follows: dPC(n,m,t) = cos(θ(n,t)−θ(m,t)), where cos(0) = 1 represents the case when the two brain areas n and m are aligned at time t (Figure 1E). Conversely cos(π) = −1 indicates the two brain areas nandm to be anti-aligned at time t. Lastly, cos(π/2) = 0 shows the two brain areas n and m at time t to be orthogonal to each other and therefore their phase relationship being 0.

Phase Dynamics Leading Eigenvector

We used LEiDA, where only the 1xN leading eigenvector V₁(t) of the dPC is considered in the analysis, to describe the phase coherence pattern of the (NxN) dPC(t) at every time-point t with reduced dimensionality. In other words, we calculated the eigendecomposition of dPC(t) at every time t [dPC(t) = V(t)D(t)V⁻¹(t), where columns of V(t) are the corresponding eigenvectors of dPC(t) and D(t) is the diagonal matrix of the eigenvalues of dPC(t), and we took the first (most dominant) eigenvectorV₁(t) to represent the BOLD PL pattern at each time point with size 1xN. Since dPC(t) is symmetric, its eigenvectors are orthogonal (V^–1(t) = V(t)^T) and the eigenvalues are real. Each element in the eigenvector can be associated to a specific brain area [i.e., in Figures 1B–D each brain area is colored according to its sign in V₁(t)]. The NxN dominant connectivity pattern at every time t captured by V₁(t) can simply be retrieved by calculating the matrix product of the eigenvector with its transpose as $V_1(t)*V_1^T(t)$ (Cabral et al., 2017b).

With the aforementioned reduction, whole-brain activity at each time point t is represented by a 1xN vector, where N is the number of brain regions defined by the applied parcellation. Each vector V₁(t) can be seen as an observation of the dynamical system and can be represented as a point in a N-dimensional space R^N (in Figure 2A, represented in R³ for illustration). Each fMRI experiment is thus characterized by a trajectory of the leading eigenvector V₁ in this N-dimensional space. To get a graphical representation (Figures 2A,B), we project each vector V₁(t) on the space determined by the first three principal components of all V₁s, i.e., the x, y, and z coordinates are given by the cosine distance between each 1xN V₁(t) and the first 3 Nx1 eigenvectors of the NxN covariance matrix of all V₁s (with size NxT).

FIGURE 2

Figure 2. Partition of the N-dimensional phase space into a K-dimensional state-based space. (A) Representation of all BOLD PL patterns captured at each TR in a reduced 3-dimensional perspective of the phase space. Since each observation is a 1xN vector – corresponding to the leading eigenvector of BOLD phases at each TR – the full phase space is N-dimensional, where N = 90 is the number of ROIs used to parcellate the brain). Each dot corresponds to one fMRI volume recorded over time (TR = Repetition Time = 0.72 s). Dots are placed according to their cosine distance with respect to the first three principal components (i.e., the first 3 eigenvectors of the covariance matrix) of all observations. (B) Partition of the Phase Space using K-means clustering decomposes the space of observations into k clusters, where each observation from an fMRI experiment is assigned to a cluster given its closest proximity to the corresponding centroid. The centroids obtained for K = 5 are represented by coloring each brain area using the same color scheme as in Figure 1, representing distinct whole-brain BOLD phase-locking patterns. (C) Illustration of the partition of a 2D plane into k Voronoi cells, where each point in a given cell is closer to its centroid than to any other centroid.

Partition of Phase Space

In order to achieve a state-based representation, the leading eigenvectors obtained from all 99 participants in the Left-Right (LR) fMRI scanning session – corresponding to a total of T = 118602 observations (99 × 1198 TRs) with N = 90 dimensions each – are partitioned into a set of discrete states. Importantly, we do not include in this partition the Right-Left (RL) fMRI scanning sessions from the same 99 participants recorded on the same day, which will serve to test the validity and consistency of the results, as described in the following section.

Given the large number of observations in this dataset, clustering algorithms relying on the TxT similarity matrix had to be discarded because of limited computational resources (i.e., computing our TxT matrix requires >100 GB of RAM). Instead, we use the k-means algorithm, which relies on an iterative process to find the solution that minimizes the distance between each (1xN) observation and the closest 1xN cluster centroid (we note that, given the large number of dimensions, we use the Cosine distance, which significantly reduces the computation time with respect to City Block or Euclidean distances). As such, k-means algorithm is used to iteratively cluster the leading eigenvectors into k = 2 to k = 20 clusters (resulting in 19 partitions), repeating each calculation 100 times to ensure stability in the results. Since each observation represents a time point, the output vector of cluster assignments – where each observation is assigned to its closest 1xN cluster centroidα = 1…k – can be approached as a trajectory $\underset{¯}{x}\left(t\right)$ in state space.

In Figure 2 we show how the k-means clustering algorithm divides the phase-space (here represented in only three dimensions for illustration) into k = 5 (Voronoi) cells, where each location in the phase space is assigned to the closest centroid. Using colors to represent the regions of the phase space assigned to each cluster R^α, we represent in Figure 2B the same observations from panel A, but highlighting the cluster assignment at each time point. Although the 3-dimensional representation serves to illustrate the partition of the phase-space, there is a clear overlap of colors given that the phase-space, defined in N = 90 dimensions, cannot be adequately represented in three dimensions only. To illustrate the decomposition into k-means clustering algorithm we show, in Figure 2C, the partition of a 2D plane into k Voronoi cells, where each region in space is assigned to the closest centroid.

Each cluster R^α(with α = 1…k) is now represented by its cluster centroid V_cα, each corresponding to a distinct BOLD phase-locking state (which will be described in detail in the section “Results”). To assess the quality of the cluster separation, the silhouette value is computed for each k, which estimates how similar each observation V₁(t) is to its own cluster compared to other clusters.

Comparison to Reference Intrinsic Functional Networks

The existence of functionally interconnected subsystems, where subsets of brains areas consistently activate together even during rest, has been widely explored in studies of resting-state functional connectivity. Intrinsic Functional Networks, typically assessed using correlation analysis, have been consistently detected in large cohorts of resting-state fMRI experiments (Yeo et al., 2011), but the analysis of their temporal dynamics has been hindered by the methods used to assess them, namely sliding window methods with their choice of the window over which connectivity is computed (Hindriks et al., 2016).

Here, we verify if the centroids obtained from clustering BOLD phase leading eigenvectors obtained at TR resolution share spatial similarities with the seven cerebral intrinsic functional networks estimated by Yeo et al. (2011) clustering correlation-based functional connectivity between 1175 regions of interest from 1000 participants.

To do so – and since our BOLD PL centroids V_α are defined in AAL parcellation – we take the mask of the Yeo parcellation into seven non-overlapping functional networks defined in MNI152 space¹ and the mask of the AAL parcellation in the same MNI152 space, and calculate, for each of the 90 AAL brain areas, the proportion of voxels assigned to each of the seven functional networks, obtaining in this way 7 1 × 90 vectors representing the intrinsic functional networks in AAL space.

Subsequently, we compute the Pearson’s correlation (with associated p-values) between these seven networks and the centroids V_α obtained from our clustering analysis across the whole range of k explored (setting all negative values of the centroids’ vectors to zero, to consider only the areas whose BOLD phase is shifted from the main orientation).

Projection of the Validation Dataset Into the Same State Space

We used the second fMRI scanning session from each of the 99 participants recorded on the same day as the primary dataset – differing only in the oblique axial acquisition phase encoding, being Right to Left (RL) instead of Left to Right (LR) – to verify the validity and consistency of the partition performed in the previous session. To do so, we obtained all the 1xN leading eigenvectors from the Right-Left (RL) scanning session – totaling 118602 observations (1198 TRs × 99 subjects) – using the same methodology as before, but instead of running the k-means algorithm, we compute the cosine distance between each 1xN eigenvector V₁(t) and the k 1xN cluster centroids V_cα obtained from the previous analysis, and define the trajectory vector $\underset{¯}{x}\left(t\right)$ by assigning each V₁(t) to its closest cluster centroid V_cα.

Fractional Occupancy

Following the cluster partition into k PL states evolving in time t, the probabilities – or fractional occupancies – ${\mathrm{\Pi }}_{\mathrm{\alpha }}^{(S)}$ associated to each PL state α and each scan S, can be calculated as follows:

$${\mathrm{\Pi }}_{\mathrm{\alpha }}^{(S)}=\frac{1}{T}\sum _{t=1}^T\mathrm{\chi }[\underset{¯}{x}\left(t\right)\in R^{\mathrm{\alpha }}](1)$$

where χ is the indicator function – χ(A) = 1 if the event A is true, and χ(A) = 0 otherwise, and T = 1198 is the number of time points (TRs) corresponding to each fMRI scan (S). In other words, the equation counts the number of times when the trajectory $\underset{¯}{x}\left(t\right)$ is assigned to each of the defined clusters R^(α), divided by the total number of time points T. Furthermore, given that participants are constantly in resting state – i.e. without performing any task -, we assume stationarity in the data within each scan (justifying the time average in Eq. 1). Cluster probabilities are estimated separately for each individual fMRI scan.

Dwell Time

To describe the average time periods when a given PL state α is being visited in each fMRI scan S, the Dwell Time $D{T}_{\mathrm{\alpha }}^{(S)}$ is defined as the mean of all the consecutive periods of each state, i.e.,

$$D{T}_{\mathrm{\alpha }}^{(S)}=\frac{1}{p_{\mathrm{\alpha }}}\sum _1^{p_{\mathrm{\alpha }}}{C}_{p_{\mathrm{\alpha }}}(2)$$

where DT_α is the Dwell Time of PL state α, p_α is the number of consecutive periods assigned to PL state α and C_pα is the duration of each consecutive period.

Markov Chain Transition Probabilities

Following the same rational as in Eq.1, the definition of the probability Π_αβ to be in the PL state α at time bin t and in the PL state β at time bin t + 1 can be written as follows:

$${\mathrm{\Pi }}_{\mathrm{\alpha }\mathrm{\beta }}^{(S)}=\frac{1}{T-1}\sum _{t=1}^{T-1}\mathrm{\chi }[\underset{¯}{x}\left(t\right)\in R^{\mathrm{\alpha }},\underset{¯}{x}(t+1)\in R^{\mathrm{\beta }}](3)$$

and thus, the transition probability matrix $W_{\mathrm{\alpha }\mathrm{\beta }}^{(S)}$ of each fMRI scan S is defined as:

$$W_{\mathrm{\alpha }\mathrm{\beta }}^{(S)}=P[\underset{¯}{x}(t+1)\in R^{\mathrm{\beta }}|\underset{¯}{x}\left(t\right)\in R^{\mathrm{\alpha }}]=\frac{{\mathrm{\Pi }}_{\mathrm{\alpha }\mathrm{\beta }}^{(S)}}{{\mathrm{\Pi }}_{\mathrm{\alpha }}^{(S)}}.(4)$$

W_αβ defines the transition matrix from state alpha to state beta. This defines an homogeneous Markov chain, characterizing the transition between BOLD phase locking states. The transition probability matrix $W_{\mathrm{\alpha }\mathrm{\beta }}^{(S)}$ is estimated separately for each scan S. To each matrix $W_{\mathrm{\alpha }\mathrm{\beta }}^{(S)}$ is associated a transition graph with an oriented arrow from α to β if $W_{\mathrm{\alpha }\mathrm{\beta }}^{(S)}$ > 0 (see Supplementary Figures S2–S4). To illustrate the transitions at the group level, we represent the transition graph of the average transition matrix W_αβ across all scans in the LR sessions.

Intra-Class Correlation

In order to calculate the reliability of the computed measures between the LR and RL fMRI sessions recorded on the same day, we calculated the Inter-Class Correlation (ICC) (Landis and Koch, 1977; Xing and Zuo, 2018). ICC describes the proportion of within-subject variability versus between-subject variability across recording conditions as follows:

$$ICC=\frac{MS{E}_b-MS{E}_w}{MS{E}_b+MS{E}_w},(5)$$

where MSE_w and MSE_b are the within-subject and between-subject mean squared errors, respectively (Xing and Zuo, 2018). Positive ICC values (i.e., when within-subject MSE is smaller than the between-subject MSE) indicate individual reliability, which, depending on its value, is categorized as low (0 < ICC < 0.2), fair (0.2 < ICC < 0.4), moderate (0.4 < ICC < 0.6), substantial (0.6 < ICC < 0.8) and almost perfect (0.8 < ICC < 1) (Landis and Koch, 1977).

Effects of Low-Pass Temporal Filtering

A typical step in the pre-procession of fMRI resting-state data is the application of a low-pass filter to exclude high frequency noise in the BOLD signal (typically < 0.1 Hz). However, given that BOLD signals are already averaged over all voxels within each brain area – which should improve the signal-to-noise ratio-, and given the instantaneous nature of our dynamic analysis, we performed our first analysis directly on the unfiltered BOLD signals recorded at a TR of 0.72 s, corresponding to a Nyquist frequency of fNq = 1/2TR = 0.694 Hz.

To verify whether the higher frequency components in the BOLD signal are meaningful for the dynamic analysis of functional networks, we apply a 2nd order Butterworth band-pass filter to the ROI-averaged BOLD signals – before computing LEiDA and clustering into k = 5 states – varying the low-pass cut-off frequency to 0.07, 0.1, 0.2 or 0.6944 Hz, while keeping the lower high-pass frequency limit at 0.01 Hz (to exclude only the ultra-slow signal drifts from the scanner). ICC measures were subsequently calculated for the corresponding Dwell Times, Fractional Occupancy and Transition Probabilities obtained.

Results

Phase-Locking States Reveal Relevant Functional Networks

We obtain a set of BOLD phase-locking patterns from the first session of resting-state fMRI of 99 unrelated subjects using the LEiDA approach (see section “Materials and Methods” for details). Each BOLD phase-locking pattern is represented as a vector with N elements, each element representing the projection of the BOLD phase of a brain area into the leading eigenvector of all BOLD phases (here N = 90 since we use the 90 non-cerebellar brain regions from the AAL atlas).

Firstly, we verify the overlap between the BOLD phase-locking states obtained across clustering solutions (with 2 < k < 20) to seven intrinsic functional networks defined in the literature (Yeo et al., 2011). In Figure 3, we report for all partitions into k states (rows), the k cluster centroids obtained (columns). The cluster centroids V_cα(representing BOLD phase-locking states) are represented in the brain by coloring only the brain areas whose BOLD phase projects in the opposite direction from the main orientation of BOLD phases (negative elements in V_cα). BOLD phase locking states are color-coded according to the most significantly correlated RSN used as reference (shown in panel B), given a corrected threshold of p < 0.05/k, and in black otherwise. The same Figure is shown from a top view perspective in Supplementary Figure S1.

FIGURE 3

Figure 3. Centroids overlap with Functional Brain Networks. (A) Representation of the centroids obtained for each k-means clustering solution with K ranging between 2 and 20. Centroids are represented in cortical space, rendering only the ROIs whose BOLD phase is shifted > ±π/2 with respect to the leading direction. ROIs are colored according to the reference Functional Brain Network (shown in panel B) to which they most significantly relate with. Pearson’s r and corresponding p-value are reported as a title only when surviving a conservative threshold of p < 0.05/K, to correct for the number of independent hypotheses tested in each partition model. Centroids not significantly overlapping with any of the reference functional networks are colored in black. Side views of the same centroids are reported in Supplementary Figure S1. (B) Reference functional brain networks estimated from 1000 subjects from correlation-based intrinsic functional connectivity (Yeo et al., 2011). (C) Silhouette value, used to evaluate clustering performance, shows a peak for partitions into 2 to 6 clusters.

Sorting the states according to their probability of occurrence, we find consistently across clustering solutions a most prevalent state, occupying the first column of Figure 3, in which the BOLD phase of all brain areas project into the same direction. Since it does not reveal the separation of any particular subsystem, and does not significantly overlap with any reference functional network, this so-called global state (state 1) is represented as a transparent brain.

The remaining states are all characterized by a phase shift in the BOLD signal of a given subset of brain areas, which are highlighted as colored patches. Notably, most of the obtained cluster centroids demonstrate a close statistical similarity to reference functional networks, revealing a strong and highly significant overlap (up to r = 0.89, with p-values down to10⁻³⁰) with the different RSNs used as reference. We also find that some partitions show different PL states overlapping with the same reference RSN. When no significant overlap is found, the patches are colored in black. One example is the second most prominent state appearing in all clustering solutions with k > 7, which involves regions of basal ganglia, which have been omitted in the analysis of functional networks from Yeo et al. (2011).

We chose to focus on the clustering solution with k = 5 within the range of best clustering performance according to the silhouette value (Figure 3C), as it reveals a meaningful partition of the BOLD PL patterns into four representative functional networks. For k = 5 we found State 2 to correlate with the Default Mode Network (r = 0.71, p = 10^–14), State 3 to correlate with the Fronto-parietal Network (r = 0.84, p = 10^–21), State 4 to correlate with the Visual Network (r = 0.88, p = 10^–29) and finally State 5 to mostly correlate with the Ventral Attention Network (r = 0.71, p = 10^–14). In Supplementary Figures S5, S6, we also report the overlap of cluster centroids with reference functional brain networks obtained for the filtered series (0.04–0.07 Hz) both from top and side view.

Exploration of a Repertoire of BOLD Phase Locking States

In Figure 4, we show the different representations of the BOLD PL states and their properties. Each PL state is represented in two ways: on the left we plot the N = 90 vector elements as arrows representing the magnitude of projection of each brain area into the leading eigenvector of BOLD phases V₁ and on the right by rendering and coloring the brain regions shifted from the main orientation (corresponding to the red arrows on the left) according to the relevant functional system to which they maximally overlap with (Figures 4A,E).

FIGURE 4

Figure 4. Repertoire of BOLD Phase Locking States obtained using K = 5. (A) BOLD Phase Locking states represent recurrent patterns of BOLD phase alignment across the whole brain. Each centroid is a vector V_C of size 1xN whose elements indicate how each brain area projects into it. Each centroid is represented in the brain in two different ways: (left) by placing an arrow at the centre of gravity of each brain area and setting its size, direction and color according to the magnitude and sign of the corresponding element in V_C (coloring in red for positive projections into V_C, and blue otherwise). (right) Rendering all brain areas with positive values in V_C colored according to the functional network to which they show maximal overlap (see Figure 3 and/or panel E below). (B) Phase-locking matrices computed as the outer product of the centroid vectors V_C. (C) Scatter plots of state Fractional Occupancy, plotting the values obtained for the 99 fMRI scans in the LR scanning session versus the values from the RL session. (D) Scatter plots of mean state Dwell Times, plotting the values obtained for the 99 fMRI scans in the LR scanning session versus the values from the RL session. (E) Correlation between each BOLD PL state and the seven networks of intrinsic functional connectivity from Yeo et al. (2011). * indicates significance corrected for multiple comparisons with Pearson’s p-value < 0.05/7.

The PL states can also be represented in the form of a matrix by computing the matrix product of each centroid’s vector V_cα and its transpose, describing the pairwise relationship between individual brain regions in each PL state (Figure 4B).

Assuming stationarity of the brain’s dynamical regime during rest, we compute the probability of occurrence of PL states as well as their mean Dwell Time within each fMRI scan (see section “Materials and Methods”). In Figures 4C,D, we show the probabilities of occurrence and mean Dwell Times of each PL state obtained for each participant, plotting the values obtained from the first fMRI session (LR) versus the values obtained from the second same-day fMRI session (RL). We find that, in both LR and RL sessions, State 1 shows high variability both in terms of probability of occurrence (mean = 0.51, standard deviation (std) = 0.16) and Dwell Times (mean = 3.94 s, std = 1.73 s), with some subjects spending as little as 20% of the time in this globally coherent state, whereas others spend up to 80% of the time, with some occurrences lasting up to 10 s (the reliability of metrics across recordings will be addressed in a following section). Interestingly, the other 4 states show consistently lower probabilities of occurrence, with state 2 (overlapping with the DMN) occurring on average 16.6 ± 7.6% of the time (mean ± std), being slightly more prevalent than the other states (state 3: 12.7 ± 6.2%; state 4: 9.9 ± 4.7%; state 5: 9.5 ± 5.5%). Not only do these functionally relevant PL patterns occur less often, but they also show, consistently across subjects, much shorter Dwell Times, lasting on average around 2 TRs (state 2, 1.71 ± 0.34 s; state 3 1.57 ± 0.37 s; state 4, 1.4 ± 0.34 s; state 5, 1.3 ± 0.22 s).

In Figure 4E, we report the correlation between each PL state and the seven intrinsic functional networks used as reference (see section “Materials and Methods”). We observed State 2 to correlate only with the Default Mode Network (r = 0.71, p = 10^–14), State 3 to correlate with the Fronto-parietal Network (r = 0.84, p = 10^–24), State 4 to correlate with the Visual Network (r = 0.88, p = 10^–29) and finally State 5 to mostly correlate with the Ventral Attention Network (r = 0.71, p = 10^–14) but also with the Somatomotor Network (r = 0.53, p = 10^–8). Supplementary Figure S7 of the states’ measures for the filtered data (0.04–0.07 Hz) is added in the Supplementary Material.

Recurrent Excursions Into BOLD PL States

Similar to the probability of occurrence of a given state, we can quantitatively characterize the temporal trajectories by the probabilities of transition between the different BOLD PL states. In Figure 5A we show the average transition matrix, W_αβ, as the probability of switching from state α to state β. We noted that the highest probabilities of transition (W_αβ > 0.5) were along the diagonal (representing the probability to remain in the same state) as well as along the first column (representing the transitions back to the state 1). The characteristic self-transitions (α → α) along the diagonal are a distinctive feature of the system, indicating the relative stability of each state. State 1 reveals the highest stability (with 77% probability of remaining in it in the following TR), whereas the probability of remaining in the other states is close to chance levels. The scatter plots in Figure 5B show the transition probabilities obtained for each of the 99 participants (LR session vs. RL session), revealing consistency of the results across participants and scanning sessions.

FIGURE 5

Figure 5. Trajectories of brain activity in state space. (A) Transition matrix quantifying the probability of the trajectory transiting from one state to another as defined in Equation 5, averaged across the scans from all 99 participants in the LR session. (B) Transition probabilities estimated for each fMRI scan, each dot corresponding to one participant, plotting the probabilities of switching in the LR fMRI session versus the RL fMRI session. (C) Transition Graph is constructed from transition matrix W where edges α→β are directed and weighted with weight W_αβ.

Another relevant feature is the asymmetry of the transition matrix, which is indicative of an imbalance in the reciprocity of transitions both to and from a given state, as can be observed with the apparent proclivity for switching into the (global) state 1, whereas the probability to leave from it is much smaller.

In Figure 5C, a Transition Graph is constructed from the transition matrix W shown in panel A, where edges α → β are directed and weighted with weight W_αβ. This gives a good insight into the spontaneous transition dynamics and motivates the use of the Markov chain transition matrix beyond the probability of occurrence alone. Supplementary Figure S8 of the transition graph and matrix for the filtered data (0.04–0.07) is added in the Supplementary Material. Furthermore, the comparison of LR and RL sessions for the probabilities of transition is added in Supplementary Figure S10.

Reliability of Individual Metrics

To assess the metric’s reliability across same-subject same-day recordings when compared to other subjects, we computed the Intra-Class Correlation (see section “Materials and Methods”) for each measure above, namely, the Fractional Occupancies, the Dwell Times and the Transition Probabilities. In Figure 6A, we show the Fractional Occupancy for all 5 states to have moderate reliability values (State1: ICC = 0.59, State 2: ICC = 0.47, State 3: ICC = 0.42, State 4: ICC = 0.39 and State 5: ICC = 0.51). The Dwell Times for the first three states had moderate values of ICC and States 4 and 5 showed poor values (State1: ICC = 0.55, State 2: ICC = 0.37, State 3: ICC = 0.4, State 4: ICC = 0.32, and State 5: ICC = 0.28, Figure 6B).

FIGURE 6

Figure 6. Individual reliability of the Phase-locking states’ measures: (A) Intra-Class Correlation (ICC) calculated for the Fractional Occupancy, showing positive ICC values for all five PL states, meaning that the within-subject error is smaller than the between-subject error. All values are within a moderate range of within-subject reliability (i.e., 0.4 < ICC < 0.6) according to the categorization by Landis and Koch (1977). (B) Intra-Class Correlation for Dwell Times for all five states showing the states 1,2 and 3 to be in the moderate reliability range. (C) ICC for the normalized Probability Transition Matrix showing positive ICC values in all transitions, with the probabilities of transition from state 1 being the most reliable, whereas other transitions, particularly between states 2 to 5, showing lower reliability. The ICC is categorized, based on (Landis and Koch, 1977) as low (0 < ICC < 0.2), fair (0.2 < ICC < 0.4), moderate (0.4 < ICC < 0.6), substantial (0.6 < ICC < 0.8) and almost perfect (0.8 < ICC < 1).

Regarding the transition matrix, Figure 6C shows that the probability of remaining in the (global) State 1 has one of the highest ICC values (ICC = 0.61), with most of the transition to and from the State 1 showing a range of fair to moderate ICC values (0.29 < ICC < 0.62). States 2 and 3 have border-line moderate values of ICC in the self-transitions (State 2: ICC = 0.39, State 3: ICC = 0.40) and some of the transitions to other states were also in the moderate range. States 4 and 5 seem to have relatively poor, but still positive, ICC values for the probability of transitions metric (Figure 6C).

Taken overall, the ICC results show that all the measures evaluated have smaller within-subject error than the between-subject error (given ICC values are positive for all measures), indicating that the measures proposed herein capture individual fingerprints of dynamic functional connectivity. To improve the assertion of individual landscapes and reliability of the methodology we added scatter plots for all the three measures (Fractional Occupancy, Dwell Times and Transition Probabilities) of the two sessions (LR and RL) in Figures 4, 5.

Effect of the Temporal Filtering

All the results shown so far were obtained directly from the ROI-averaged BOLD signals from the HCP dataset, without applying any temporal frequency filter. Temporal filtering is a typical pre-processing step in resting-state fMRI analysis to remove frequency components regarded as noise. In this section, we evaluate whether the inclusion of the higher frequency components in the BOLD signal improve the analysis of dynamic functional connectivity by evaluating its effects on the reliability (ICC) of the measures across sessions, which should be maximized if assuming stationarity in individual resting-state dynamics.

As shown in Figures 7A,C, filtering has a crucial effect on the Dwell Times, with lower cut-off frequencies leading to longer Dwell Times for all states, and especially for state 1. Notably, when reaching up to the Nyquist frequency, the mean Dwell Times of states 2 to 5 approach the duration of 2TRs consistently for all subjects and in both LR and RL fMRI sessions, while state 1 lasts for slightly longer periods. When evaluating the Dwell Times ICC (Figure 7E), we find that the ICC is maximal for states 1, 2, and 3 when the high frequency components of the BOLD signal are included. However, it is important to take into account that the accurate estimation of Dwell Times is limited by the temporal resolution of the current fMRI dataset, minimizing the difference between subjects and hence affecting the ICC estimation.

FIGURE 7

Figure 7. Temporal Filtering effect on Dwell Time and Fractional Occupancy. (A,C) Dwell Times obtained across four different band-pass filters applied to the BOLD signals, keeping the lower bound fixed at 0.01 Hz and varying the higher (lowpass) cutoff frequency between 0.07, 0.1, 0.2, and 0.6944 Hz, which corresponds to the Nyquist frequency fNq = 1/(2TR). Results are shown for LR (A) and RL (C) recording sessions. (B,D) Fractional Occupancies across the four different filtrations for LR (B) and RL (D) recording sessions. (E,F) Inter-Class Correlation (ICC) across four different filtrations for Dwell Times (E) and Fractional Occupancy (F). The red dashed line represents the threshold for moderate reliability based on the Landis and Koch (1977) ICC scale.

Regarding the Fractional Occupancy of the states, it appears from Figures 7B,D that filtering does not affect the overall values estimated across subjects and in the different sessions. However, when looking at the ICC values for Fractional Occupancy (Figure 7F), we find that this measure is much more reliable within individuals (and across all states) when the high frequency components of the BOLD signal are included. Similarly, the ICC of the Transition Probabilities (shown in Supplementary Figure S9) shows a slight improvement in the reliability of most transitions. Actually, when considering only frequencies <0.2 Hz, a few transition probabilities have even a negative ICC, indicating no individual reliability at all. Overall, these findings point to the direction that it is important to consider the high frequency components of resting-state BOLD signals when assessing individual fingerprints in dynamic functional connectivity.

Discussion

The challenge of describing dynamic functional connectivity for a mechanistic understanding of the brain processing as well as for its potential use in clinical research, has been of great interest to the neuroimaging community (Hutchison et al., 2013; Preti et al., 2016). With the advent of open multimodal neuroimaging data, it is possible to address and validate these approaches in representative datasets (Van Essen et al., 2013; Poldrack and Gorgolewski, 2014). In this work, we apply, formalize and validate the Leading Eigenvector Dynamics Analysis (LEiDA) to a large cohort of 99 healthy unrelated HCP subjects (Glasser et al., 2013). We describe brain activity during rest as a time evolving trajectory in a low-dimensional state space, where states are defined according to characteristic whole-brain configurations of BOLD phase-locking. Furthermore, we validate these BOLD phase-locking states to reference networks of intrinsic functional connectivity (Yeo et al., 2011) and compute their properties of fractional occupancy, Dwell Times and transition probabilities. We subsequently assess the reliability of these measures across the two same-day fMRI recordings (using Intra-Class Correlation) and show that all measures have a smaller within-subject error than the between-subject error (ICC values > 0), with the highest reliability values being detected when including the high frequency components (>0.1 Hz) of ROI-averaged BOLD signals in the analysis. We argue that such interpretation of brain activity, validated with reliability analysis, has the potential to identify individual-specific fingerprints in the brain’s dynamical landscape and thus serve personalized clinical applications in diagnostics and therapeutics of patients with cognitive disorders.

Concepts and methods from dynamical systems theory are proving useful in the analysis of brain activity at the macroscopic scale, as they serve to formally characterize the complex dynamics emerging from the collective behavior of billions of interacting neurons, exhibiting features such as multi-stability, meta-stability and self-organized criticality, that may serve helpful to identify the underlying principles coordinating cognition at the whole-brain level (Deco and Jirsa, 2012; Tognoli and Kelso, 2014; Cocchi et al., 2017; Roberts et al., 2019). Here, we aimed at a different characterization of the dynamical properties of the intrinsic functional networks emerging spontaneously and consistently during rest. Our analysis revealed a repertoire of BOLD phase-locking states through which the trajectory of brain activity consistently returns in time and across subjects.

By analyzing the PL state’s fractional occupancy, Dwell Time and probability of transitions, our results revealed that the BOLD phase-locking states can be divided in two groups according to their dynamical properties: On one hand, our algorithm consistently detects a state where all the ROI-averaged BOLD signals project into the same direction captured by the leading eigenvector (state 1 for all clustering solutions). This state exhibits longer Dwell Times and shows high between-subject variability but also high within-subject reliability. On the other hand, we detect for all k-means clustering solutions, a set of k-1 states where the BOLD signals of some brain areas project into the opposite direction from the main BOLD phase orientation. These states occur consistently less often and last for shorter times than the global state, but reoccur consistently across subjects and sessions. Given the reduced stability of these BOLD phase-locking states with respect to the meta-stable globally synchronized state, we refer to this second group as “ghost” attractor states. In other words, ghost attractors in this framework refer to short-lived (or weakly stable) network configurations that consistently reoccur across fMRI recordings.

Regarding the functional relevance of these “ghost” phase-locking states, our results show a clear and highly significant overlap of most cluster centroids (obtained for the whole range of partitions explored) with a set of seven previously identified networks of intrinsic functional connectivity used as reference. This finding indicates that these patterns of BOLD phase locking, despite being obtained from a different analytic perspective than more conventional correlation-based analyzes, are closely related to the so-called resting-state networks. Yet, unlike correlation-based analyses that reveal only the spatial map of these functional networks, the LEiDA approach allows characterizing their properties over the temporal dimension. As the reference RSNs are computed from the correlation-based static functional connectivity, a perfect match to the BOLD phase-locking states detected herein is not expected. Rather, they can be considered for validation of the functional relevance of the PL states and served to guide the choice of the number of states or further analysis. The number of states chosen is a trade-off between more fine-grained but less robust state solutions as demonstrated by the increasing specificity of functional subsystems for higher k. Here, the clustering solution with k = 5 was chosen for being within the range of maximal Silhouette value and for revealing a separation into distinct functionally meaningful systems such as the Default Mode Network, the Frontoparietal Network, the Ventral Attention Network and the Visual Network. However, a partition into a higher number of states may prove necessary when addressing particular conditions that affects a particular subsystem optimally defined for higher k. For instance, in a previous work using LEiDA, the partition into k = 10 was chosen for detecting the network that most significantly distinguished patients in remission from major depressive disorder and controls (Figueroa et al., 2019), whereas another study found the solution with k = 7 to optimally highlight the effects of psilocybin (Lord et al., 2019).

For all partitions into k > 7, our algorithm consistently detected a functional subsystem involving the basal ganglia (colored in black in Figure 3 and Supplementary Figure S1) for not overlapping with any of the reference RSNs) as the second most prevalent BOLD PL state. This indicates that resting-state activity also involves connectivity to subcortical areas, which appears particularly important for the study of psychiatric disorders, such as anxiety-related disorders involving the basal ganglia. Following previous LEiDA studies, we chose here a coarse parcellation into N = 90 brain areas and did not include the BOLD signal detected in the cerebellum. The Anatomic Automatic Labeling Atlas has been validated in many studies and has shown consistency in the LEiDA results across datasets (Cabral et al., 2017b; Figueroa et al., 2019; Lord et al., 2019). However, it is based on an anatomic definition of the brain regions and as such might not generalize adequately to the dynamic functional connectivity analysis. We expect to extend to finer-grained and fMRI-derived parcellations in future studies, potentially including other substructures such as the cerebellum, in order to gain a wider insight into the network configurations observed in brain activity at the macroscopic scale (Cammoun et al., 2012; Glasser et al., 2016; Schaefer et al., 2018).

The mechanistic interpretation of the empirical data proposed herein serves as a great candidate for further theoretical exploration by whole-brain computational models (Ghosh et al., 2008; Deco et al., 2009; Honey et al., 2009). To this date, many models have demonstrated well-matched dynamics to the brain activity as represented by static functional connectivity in a critical range of parameters where the brain is poised between noisy and oscillatory activity. Furthermore, different properties were shown to have an impact on the emerging dynamics such as propagation delays, local vs. global connections, signal-to-noise ratio, and local inhibitory rules (Deco et al., 2009, 2014; Honey et al., 2009; Cabral et al., 2011). Extending such modeling endeavors away from static functional connectivity to a dynamic representation of the experimental data is currently becoming a possible new avenue into understanding the underlying principles governing dynamic functional connectivity (Hansen et al., 2015; Cabral et al., 2017a; Deco et al., 2018, 2019). Recently, Deco et al. have shown how the dynamic representation of resting-state data in wakefulness and sleep (characterized using LEiDA) can serve to explore how a whole-brain model can be perturbed to identify the brain regions responsible for the transition between awake and sleep state (Deco et al., 2019).

Representing dynamic functional connectivity through the prism of dynamical system theory hypothesizes the existence of attractors in N-dimensional space through which the functional activity evolves in time. Assuming this hypothetical scenario, it describes a state-based propagation of the data, rendering the underlying dynamics in a discrete sense (Baker et al., 2014; Karahanoğlu and Van De Ville, 2015; Preti et al., 2016; Cabral et al., 2017a). However, other methods have considered dynamic functional connectivity from a continuous point-of-view, such as the spatio-temporal connectome where brain activity is described as a temporal graph (Griffa et al., 2017; Vohryzek et al., 2019) and auto-regressive models (Liégeois et al., 2017). We acknowledge that looking at the brain activity in a discrete sense is only one of the interpretations currently proposed in trying to describe the emergent complex phenomena observed in whole-brain dynamics.

It is to be noted that the applied clustering algorithm is just one amongst many decomposition methods that can partition the LEiDA results into meaningful states. Indeed (Cabral et al., 2017b) compared the results from k-means algorithm to the Hidden Markov Model (HMM), in their paper on cognitive performance of patients, showing similar results with both approaches (Cabral et al., 2017b). However, k-means was chosen here for its relatively simple implementation and its relatively low computational cost, revealing functionally meaningful cluster centroids.

New imaging methods benefit greatly from the reliability analysis that investigates individual variabilities across recordings sessions. Especially in clinical applications, reliability is crucial to obtain stable measures across time for individual subjects (i.e., low within-subject variability) and at the same time distinguishable differences between subjects (i.e., high between-subject variability) (Xing and Zuo, 2018; Zuo et al., 2019). In this work, Intra-Class Correlation is used to calculate the desirable ratio between between-subject variability and within-subject variability across recording sessions. One of the Intra-Class Correlation scales proposed by Landis and Koch (1977) to assess reliability for clinical applications suggests that values 1.0 > ICC > 0.8 have excellent reliability, 0.8 > ICC > 0.6 substantial reliability, 0.6 > ICC > 0.4 have moderate reliability and, 0.4 > ICC > 0.2 poor substantial reliability. In other words, it is desirable to obtain high reliability values for the method’s possible clinical application.

In the last part of our study, we show that including the high frequency components up to the Nyquist frequency maximizes the ICC values (reaching a mainly moderate range of ICC values). As such, it is likely that the temporal resolution of the fMRI acquisition might have hindered further increase in reliability. Although the Dwell Times become significantly shorter if no smoothing is applied – which may decrease the detection of RSNs in correlation-based analysis – we find that the occurrence of these states is intrinsically short, given that the measures become more reliable. Although the hemodynamic response function (HRF) is intrinsically slow, the capacity of the BOLD signal to detect faster frequency components is still highly debated in the literature (Glerean et al., 2012; Deco et al., 2019). Nevertheless, it is likely that resting-state dynamics occurs at a faster time scale than captured with the BOLD signal, as suggested by MEG studies that point to a duration of around 200 ms (Baker et al., 2014; Vidaurre et al., 2016). Here, we show that LEiDA allows detecting meaningful dynamic network configurations occurring at relatively short time-scales for fMRI analysis, which may serve useful not only for resting-state analysis but also for the detection of task related patterns [as in Stark et al. (2019)], that may not be captured with conventional general linear models using the HRF. Overall, we expect that novel insights into BOLD signal temporal characteristics and improvements in fMRI temporal resolution might increase the ICC reliability of these measures.

Conclusion

In summary, we combine novel analytic tools to quantitatively characterize brain activity in each fMRI scan as a trajectory through a discrete set of BOLD phase-locking states. Given the dynamical properties of these states (fractional occupancy, Dwell Time and transition probability) we propose that RSNs behave as ghost attractors, emerging spontaneously and for brief periods, but recurring consistently across subjects and sessions. Our study corroborates previous theoretical works that put forward an interpretation of brain activity as a trajectory evolving in time in an energy landscape (Deco et al., 2011; Ashourvan et al., 2017). By demonstrating the functional relevance of the BOLD phase-locking states detected and the reliability of the measures across same-subject sessions we go further by revealing the existence of individual-specific energy landscapes in brain activity with potential application in patient-specific diagnostics and therapeutics.

Data Availability Statement

All data and codes used in this study are publicly available on https://github.com/jvohryzek/GhostAttractors.

Ethics Statement

The studies involving human participants were reviewed and approved by The Washington University review board including all study protocols. Informed consent was obtained for all participants. The participants provided their written informed consent to participate in this study.

Author Contributions

JV and JC carried out the analysis, wrote the main manuscript and came up with the paper’s ideas. BC and GD verified and advised the theoretical methods. MK and JC supervised the whole project. All authors participated in the discussion of the ideas and contributed in the final writing of the manuscript.

Funding

This work has been funded by FEDER through the Competitiveness Factors Operational Program (COMPETE), by National funds through the Foundation for Science and Technology (FCT) under the scope of the project UID/Multi/50026; and by the projects NORTE-01-0145-FEDER-000013 and NORTE-01-0145-FEDER-000023, supported by the Northern Portugal Regional Operational Programme (NORTE 2020), under the Portugal 2020 Partnership Agreement, through the European Regional Development Fund (FEDER). JC was supported by Portuguese Foundation for Science and Technology CEECIND/03325/2017, Portugal. GD acknowledges funding from the European Union’s Horizon 2020 FET Flagship Human Brain Project under Grant Agreement 785907 HBP SGA2, the Spanish Ministry Project PSI2016-75688-P (AEI/FEDER) and the Catalan Research Group Support 2017 SGR 1545. MK was supported by the European Research Council Consolidator Grant: CAREGIVING (615539), Pettit Foundation, Carlsberg Foundation and Center for Music in the Brain, funded by the Danish National Research Foundation (DNRF117). BC was supported by the French Government through the UCA-Jedi project managed by the National Research Agency (ANR-15-IDEX-01) and, in particular, by the interdisciplinary Institute for Modeling in Neuroscience and Cognition (NeuroMod) of the Université Côte d’Azur.

Conflict of Interest

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

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnsys.2020.00020/full#supplementary-material

Footnotes

1. ^ surfer.nmr.mgh.harvard.edu/fswiki/CorticalParcellation_Yeo2011

References

Allen, E. A., Damaraju, E., Plis, S. M., Erhardt, E. B., Eichele, T., and Calhoun, V. D. (2014). Tracking whole-brain connectivity dynamics in the resting state. Cerebral. Cortex 24, 663–676. doi: 10.1093/cercor/bhs352

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashourvan, A., Gu, S., Mattar, M. G., Vettel, J. M., and Bassett, D. S. (2017). The energy landscape underpinning module dynamics in the human brain connectome. NeuroImage 157, 364–380. doi: 10.1016/j.neuroimage.2017.05.067

PubMed Abstract | CrossRef Full Text | Google Scholar

Baker, A. P., Brookes, M. J., Rezek, I. A., Smith, S. M., Behrens, T., Smith, P. J. P., et al. (2014). Fast transient networks in spontaneous human brain activity. ELife 2014, 1–18. doi: 10.7554/eLife.01867

PubMed Abstract | CrossRef Full Text | Google Scholar

Beckmann, C. F., DeLuca, M., Devlin, J. T., and Smith, S. M. (2005). Investigations into resting-state connectivity using independent component analysis. Philos. Trans. R. Soc. BBiol. Sci. 360, 1001–1013. doi: 10.1098/rstb.2005.1634

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, J., Hugues, E., Sporns, O., and Deco, G. (2011). Role of local network oscillations in resting-state functional connectivity. NeuroImage 57, 130–139. doi: 10.1016/j.neuroimage.2011.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, J., Kringelbach, M. L., and Deco, G. (2017a). Functional connectivity dynamically evolves on multiple time-scales over a static structural connectome: models and mechanisms. NeuroImage 160, 84–96. doi: 10.1016/j.neuroimage.2017.03.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, J., Vidaurre, D., Marques, P., Magalhães, R., Silva Moreira, P., Miguel Soares, J., et al. (2017b). Cognitive performance in healthy older adults relates to spontaneous switching between states of functional connectivity during rest. Sci. Rep. 7:5135. doi: 10.1038/s41598-017-05425-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Calhoun, V. D., Miller, R., Pearlson, G., and Adali, T. (2014). The chronnectome: time-varying connectivity networks as the next frontier in fmri data discovery. Neuron 84, 262–274. doi: 10.1016/j.neuron.2014.10.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Cammoun, L., Gigandet, X., Meskaldji, D., Thiran, J. P., Sporns, O., Do, K. Q., et al. (2012). Mapping the human connectome at multiple scales with diffusion spectrum MRI. J. Neurosci. Methods 203, 386–397. doi: 10.1016/j.jneumeth.2011.09.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, C., and Glover, G. H. (2010). Time-frequency dynamics of resting-state brain connectivity measured with fMRI. NeuroImage 50, 81–98. doi: 10.1016/j.neuroimage.2009.12.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Cocchi, L., Gollo, L. L., Zalesky, A., and Breakspear, M. (2017). Criticality in the brain: a synthesis of neurobiology, models and cognition. Prog. Neurobiol. 158, 132–152. doi: 10.1016/j.pneurobio.2017.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Damoiseaux, J. S., Rombouts, S. A. R. B., Barkhof, F., Scheltens, P., Stam, C. J., Smith, S. M., et al. (2006). Consistent resting-state networks across healthy subjects. Proc. Natl. Acad. Sci. U.S.A. 103, 13848–13853. doi: 10.1073/pnas.0601417103

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Cruzat, J., Cabral, J., Tagliazucchi, E., Laufs, H., Logothetis, N. K., et al. (2019). Awakening: predicting external stimulation to force transitions between different brain states. Proc. Natl. Acad. Sci. U.S.A. 116, 18088–18097. doi: 10.1073/pnas.1905534116

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Cruzat, J., and Kringelbach, M. L. (2019). Brain Songs: Discovering the Relevant Timescale of the Human Brain. Nature Communications. Berlin: Nature Publishing Group.

Google Scholar

Deco, G., Jirsa, V., Mcintosh, A. R., Sporns, O., and Kotter, R. (2009). Key role of coupling, delay, and noise in resting brain fluctuations. Proc. Natl. Acad. Sci. U.S.A. 106, 10302–10307. doi: 10.1073/pnas.0901831106

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., and Jirsa, V. K. (2012). Ongoing cortical activity at rest: criticality, multistability, and ghost attractors. J. Neurosci. 32, 3366–3375. doi: 10.1523/JNEUROSCI.2523-11.2012

CrossRef Full Text | Google Scholar

Deco, G., Jirsa, V. K., and McIntosh, A. R. (2011). Emerging concepts for the dynamical organization of resting-state activity in the brain. Nat. Rev. Neurosci. 12, 43–56. doi: 10.1038/nrn2961

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Ponce-Alvarez, A., Hagmann, P., Romani, G., Mantini, D., and Corbetta, M. (2014). How local excitation-inhibition ratio impacts the whole brain dynamics. J. Neurosci. 34, 7886–7898. doi: 10.1523/JNEUROSCI.5068-13.2014

CrossRef Full Text | Google Scholar

Deco, G., Ponce-Alvarez, A., Mantini, D., Romani, G. L., Hagmann, P., and Corbetta, M. (2013). Resting-state functional connectivity emerges from structurally and dynamically shaped slow linear fluctuations. J. Neurosci. 33, 11239–11252. doi: 10.1523/JNEUROSCI.1091-13.2013

CrossRef Full Text | Google Scholar

Deco, G., Vidaurre, D., and Kringelbach, M. L. (2019). Revisiting the global workspace: orchestration of the functional hierarchical organisation of the human brain. BioRxiv [preprint]. doi: 10.1101/859579

CrossRef Full Text | Google Scholar

Deco, G., Cruzat, J., Cabral, J., Knudsen, G. M., Carhart-Harris, R. L., Whybrow, P. C., et al. (2018). Whole-brain multimodal neuroimaging model using serotonin receptor maps explains non-linear functional effects of LSD. Curr. Biol. 28, 3065.e6–3074.e6. doi: 10.1016/j.cub.2018.07.083

PubMed Abstract | CrossRef Full Text | Google Scholar

Figueroa, C. A., Cabral, J., Roel, M. J. T., Rapuano, K. M., Van Hartevelt, T. J., and Ruhé, H. G. (2019). Altered ability to access a clinically relevant control network in patients remitted from major depressive disorder. Hum. Brain Mapp. 40, 2771–2786. doi: 10.1002/hbm.24559

PubMed Abstract | CrossRef Full Text | Google Scholar

Fox, M. D., and Greicius, M. (2010). Clinical applications of resting state functional connectivity. Front. in Syst. Neurosci. 4:19. doi: 10.3389/fnsys.2010.00019

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghosh, A., Rho, Y., Mcintosh, A. R., Kö Tter, R., and Jirsa, V. K. (2008). Noise during rest enables the exploration of the brain’s dynamic repertoire. PLoS Comput. Biol. 4:1000196. doi: 10.1371/journal.pcbi.1000196

PubMed Abstract | CrossRef Full Text | Google Scholar

Glasser, M. F., Coalson, T. S., Robinson, E. C., Hacker, C. D., Harwell, J., Yacoub, E., et al. (2016). A multi-modal parcellation of human cerebral cortex. Nature 536, 171–178. doi: 10.1038/nature18933

PubMed Abstract | CrossRef Full Text | Google Scholar

Glasser, M. F., Sotiropoulos, S. N., Wilson, J. A., Coalson, T. S., Fischl, B., Andersson, J. L., et al. (2013). The minimal preprocessing pipelines for the human connectome project. NeuroImage 80, 105–124. doi: 10.1016/j.neuroimage.2013.04.127

PubMed Abstract | CrossRef Full Text | Google Scholar

Glerean, E., Salmi, J., Lahnakoski, J. M., Jääskeläinen, I. P., and Sams, M. (2012). Functional magnetic resonance imaging phase synchronization as a measure of dynamic functional connectivity. Brain Connectivity 2, 91–101. doi: 10.1089/brain.2011.0068

PubMed Abstract | CrossRef Full Text | Google Scholar

Greicius, M. (2008). Resting-state functional connectivity in neuropsychiatric disorders. Curr. Opin. Neurol. 21, 424–430. doi: 10.1097/WCO.0b013e328306f2c5

PubMed Abstract | CrossRef Full Text | Google Scholar

Griffa, A., Ricaud, B., Benzi, K., Bresson, X., Daducci, A., Vandergheynst, P., et al. (2017). Transient networks of spatio-temporal connectivity map communication pathways in brain functional systems. NeuroImage 155, 490–502. doi: 10.1016/j.neuroimage.2017.04.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Haimovici, A., Tagliazucchi, E., Balenzuela, P., and Chialvo, D. R. (2013). Brain organization into resting state networks emerges at criticality on a model of the human connectome. Physical. Rev. Lett. 110:178101. doi: 10.1103/PhysRevLett.110.178101

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, E. C. A., Battaglia, D., Spiegler, A., Deco, G., and Jirsa, V. K. (2015). Functional connectivity dynamics: modeling the switching behavior of the resting state. NeuroImage 105, 525–535. doi: 10.1016/j.neuroimage.2014.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Hindriks, R., Adhikari, M. H., Murayama, Y., Ganzetti, M., Mantini, D., Logothetis, N. K., et al. (2016). Can sliding-window correlations reveal dynamic functional connectivity in resting-state fMRI? NeuroImage 127, 242–256. doi: 10.1016/j.neuroimage.2015.11.055

PubMed Abstract | CrossRef Full Text | Google Scholar

Honey, C. J., Sporns, O., Cammoun, L., Gigandet, X., Thiran, J.-P., Meuli, R., et al. (2009). Predicting human resting-state functional connectivity from structural connectivity. Proc. Natil. Acad. Sci. U.S.A. 106, 2035–2040. doi: 10.1073/pnas.0811168106

PubMed Abstract | CrossRef Full Text | Google Scholar

Hutchison, R. M., Womelsdorf, T., Allen, E. A., Bandettini, P. A., Calhoun, V. D., Corbetta, M., et al. (2013). Dynamic functional connectivity: promise, issues, and interpretations. NeuroImage 80, 360–378. doi: 10.1016/j.neuroimage.2013.05.079

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaiser, R. H., Andrews-Hanna, J. R., Wager, T. D., and Pizzagalli, D. A. (2015). Large-scale network dysfunction in major depressive disorder. JAMA Psychiatry 72, 603611. doi: 10.1001/jamapsychiatry.2015.0071

PubMed Abstract | CrossRef Full Text | Google Scholar

Karahanoğlu, F. I., and Van De Ville, D. (2015). Transient brain activity disentangles fMRI resting-state dynamics in terms of spatially and temporally overlapping networks. Nat. Commun. 6:7751. doi: 10.1038/ncomms8751

PubMed Abstract | CrossRef Full Text | Google Scholar

Landis, J. R., and Koch, G. G. (1977). The measurement of observer agreement for categorical data. Biometrics 33, 159–174. doi: 10.2307/2529310

CrossRef Full Text | Google Scholar

Liégeois, R., Laumann, T. O., Snyder, A. Z., Zhou, J., and Yeo, B. T. T. (2017). Interpreting temporal fluctuations in resting-state functional connectivity MRI. NeuroImage 163, 437–455. doi: 10.1016/J.NEUROIMAGE.2017.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lord, L. D., Expert, P., Atasoy, S., Roseman, L., Rapuano, K., Lambiotte, R., et al. (2019). Dynamical exploration of the repertoire of brain networks at rest is modulated by psilocybin. NeuroImage 199, 127–142. doi: 10.1016/j.neuroimage.2019.05.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Poldrack, R. A., and Gorgolewski, K. J. (2014). Making big data open: data sharing in neuroimaging. Nat. Neurosci. 17, 1510–1517. doi: 10.1038/nn.3818

PubMed Abstract | CrossRef Full Text | Google Scholar

Preti, M. G., Bolton, T. A., and Ville, V. D. D. (2016). The dynamic functional connectome: state-of-the-art and perspectives. NeuroImage 160, 41–45. doi: 10.1016/j.neuroimage.2016.12.061

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, J. A., Gollo, L. L., Abeysuriya, R. G., Roberts, G., Mitchell, P. B., Woolrich, M. W., et al. (2019). Metastable brain waves. Nat. Commun. 10:1056. doi: 10.1038/s41467-019-08999-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Schaefer, A., Kong, R., Gordon, E. M., Laumann, T. O., Zuo, X.-N., Holmes, A. J., et al. (2018). Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity MRI. Cerebral. Cortex 28, 3095–3114. doi: 10.1093/cercor/bhx179

CrossRef Full Text | Google Scholar

Sporns, O., Tononi, G., and Edelman, G. M. (2000). Connectivity and Complexity: The Relationship Between Neuroanatomy and Brain Dynamics. Neural Networks. Amsterdam: Elsevier Science Ltd.

Google Scholar

Stark, E. A., Cabral, J., Riem, M. M. E., Van Ijzendoorn, M. H., Stein, A., and Kringelbach, M. L. (2019). The power of smiling: the adult brain networks underlying learned infant emotionality. Cerebral. Cortex 00, 1–11. doi: 10.1093/cercor/bhz219

PubMed Abstract | CrossRef Full Text | Google Scholar

Tagliazucchi, E., Siniatchkin, M., Laufs, H., and Chialvo, D. R. (2016). The voxel-wise functional connectome can be efficiently derived from co-activations in a sparse spatio-temporal point-process. Front. Neurosci. 10:38. doi: 10.3389/fnins.2016.00381

CrossRef Full Text | Google Scholar

Tognoli, E., and Kelso, J. A. S. (2014). The metastable brain. Neuron 81, 35–48. doi: 10.1016/j.neuron.2013.12.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Tononi, G., and Edelman, G. M. (1998). Consciousness and complexity. Science 282, 1846–1851. doi: 10.1126/science.282.5395.1846

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Essen, D. C., Smith, S. M., Barch, D. M., Behrens, T. E. J., Yacoub, E., and Ugurbil, K. (2013). The WU-minn human connectome project: an overview. NeuroImage 80, 62–79. doi: 10.1016/j.neuroimage.2013.05.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Vargas, C., López-Jaramillo, C., and Vieta, E. (2013). A systematic literature review of resting state network-functional MRI in bipolar disorder. J. Affect. Disord.. 150, 727–735. doi: 10.1016/j.jad.2013.05.083

PubMed Abstract | CrossRef Full Text | Google Scholar

Vidaurre, D., Quinn, A. J., Baker, A. P., Dupret, D., Tejero-Cantero, A., and Woolrich, M. W. (2016). Spectrally resolved fast transient brain states in electrophysiological data. NeuroImage 126, 81–95. doi: 10.1016/j.neuroimage.2015.11.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Vohryzek, J., Griffa, A., Mullier, E., Friedrichs-Maeder, C., Sandini, C., Schaer, M., et al. (2019). Dynamic spatio-temporal patterns of brain connectivity reorganize across development. Network Neurosci. 4, 115–133. doi: 10.1162/netn_a_00111

PubMed Abstract | CrossRef Full Text | Google Scholar

Xing, X. X., and Zuo, X. N. (2018). The anatomy of reliability: a must read for future human brain mapping. Sci. Bull. 63, 1606–1607. doi: 10.1016/j.scib.2018.12.010

CrossRef Full Text | Google Scholar

Yeo, B. T. T., Krienen, F. M., Sepulcre, J., Sabuncu, M. R., Lashkari, D., Hollinshead, M., et al. (2011). The organization of the human cerebral cortex estimated by intrinsic functional connectivity. J. Neurophysiol. 106, 1125–1165. doi: 10.1152/jn.00338.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Zuo, X. N., Xu, T., and Milham, M. P. (2019). Harnessing reliability for neuroscience research. Nat. Hum. Behav. 3, 768–771. doi: 10.1038/s41562-019-0655-x

PubMed Abstract | CrossRef Full Text | Google Scholar