Graph-theoretical analysis for energy landscape reveals the organization of state transitions in the resting-state human cerebral cortex
Jiyoung Kang aff001; Chongwon Pae aff002; Hae-Jeong Park aff001
Authors place of work:
Center for Systems and Translational Brain Sciences, Institute of Human Complexity and Systems Science, Yonsei University, Seoul, Republic of Korea
aff001; Department of Nuclear Medicine, Yonsei University College of Medicine, Seoul, Republic of Korea
aff002; BK21 PLUS Project for Medical Science, Yonsei University College of Medicine, Seoul, Republic of Korea
aff003; Department of Cognitive Science, Yonsei University, Seoul, Republic of Korea
Published in the journal:
PLoS ONE 14(9)
The resting-state brain is often considered a nonlinear dynamic system transitioning among multiple coexisting stable states. Despite the increasing number of studies on the multistability of the brain system, the processes of state transitions have rarely been systematically explored. Thus, we investigated the state transition processes of the human cerebral cortex system at rest by introducing a graph-theoretical analysis of the state transition network. The energy landscape analysis of brain state occurrences, estimated using the pairwise maximum entropy model for resting-state fMRI data, identified multiple local minima, some of which mediate multi-step transitions toward the global minimum. The state transition among local minima is clustered into two groups according to state transition rates and most inter-group state transitions were mediated by a hub transition state. The distance to the hub transition state determined the path length of the inter-group transition. The cortical system appeared to have redundancy in inter-group transitions when the hub transition state was removed. Such a hub-like organization of transition processes disappeared when the connectivity of the cortical system was altered from the resting-state configuration. In the state transition, the default mode network acts as a transition hub, while coactivation of the prefrontal cortex and default mode network is captured as the global minimum. In summary, the resting-state cerebral cortex has a well-organized architecture of state transitions among stable states, when evaluated by a graph-theoretical analysis of the nonlinear state transition network of the brain.
Physical sciences – Chemistry – Physical chemistry – Reaction dynamics – Transition state – Mathematics – Systems science – Nonlinear dynamics – Nonlinear systems – Probability theory – Probability distribution – Biology and life sciences – Anatomy – Brain – Cerebral cortex – Neuroscience – Neural networks – Brain mapping – Functional magnetic resonance imaging – Neuroimaging – Medicine and health sciences – Diagnostic medicine – Diagnostic radiology – Magnetic resonance imaging – Radiology and imaging – Computer and information sciences – Network analysis – Research and analysis methods – Imaging techniques
A dynamic complex system can possess several stable states (attractors) for a given set of system parameters [1–6]. If a system has multiple coexisting stable states and can switch among them in response to noise or intrinsic perturbations to the system, it is generally referred to as a multistable system [7, 8]. In this respect, the brain at rest can be considered as a system showing multistability [1–7, 9, 10]. From the perspective of the multistable brain, the conventional term “resting state” is not a homogeneous state but a period of switching among multiple micro-states (or sub-states). Here, we will refer to a brain state as a sub-state during the resting-state period. Of note, this multistablity perspective differs from studies on functional connectivity dynamics [11–20], which have described the dynamic nature of the brain in terms of temporal changes in its interactions (connectivity parameters). In contrast, multistability in the complex system is an emergent property of nonlinear interactions among nodes in the system without any changes in their connectivity.
In a multistable system, stable states and the transition processes among them characterize the dynamics of the system. To explore the multistability and state transitions in the dynamic brain, energy landscape analysis has recently been applied to fluctuations of blood oxygenation level-dependent (BOLD) functional magnetic resonance imaging (fMRI) [21–27]. Prior to its introduction to the brain research field, energy landscape analysis had already shown its utility in understanding the dynamics of multi-dimensional complex systems, such as protein dynamics and the thermodynamics of liquids [28–33]. In the studies of brain dynamics using energy landscape analysis, distributed activity patterns across brain regions have often been used to define brain states, one of which the brain belongs to at each measurement time point [21–27]. In the energy landscape analysis, the energy of a state is the negative log probability of the occurrence of the state (thus, frequent states have low energy) according to the Boltzmann distribution of the state. The (inverse) frequency distribution of all possible brain states (patterns of brain activities across the brain regions) is called an energy landscape (see Fig 1).
The energy landscape of the system consists of several valleys with local minima (called “stable states” or “attractors”, abbreviated as LM) that have energies lower (more frequent) than their neighbors do in the valleys. Thus, the dynamics of the system can be divided into intravalley (within the basin of a local minimum) and intervalley (between local minima) motions. In the former case, a state of the system wanders around a local minimum of the energy surface that the state belongs to, whereas in the latter case, a state transits from one local minimum to another, surpassing an energy barrier.
For a state transition between two local minima states, an optimal pathway refers to the path with the lowest maximal energy barrier among all possible paths. The optimal path may contain “intermediate states” (a type of local minima) and “transition states” (saddle point states) along the path. Among many transition states along the path, the transition state with the highest energy on the pathway determines the transition rate. Therefore, for brevity, we refer to this rate-determining transition state (having the highest energy on the pathway) as the transition state (TS) between two states (See Fig 1). Table 1 summarizes the terminology used in the current paper.
Despite a growing number of studies on the multistability of the resting brain systems [1, 3–6, 21, 34–37], the state transition processes between local minima in the brain systems have not yet been sufficiently investigated. In the present study, we explored the multistability and state transitional properties of the human cerebral cortex system. We estimated an energy landscape of brain states using a pairwise maximum entropy model (MEM) of the resting-state fMRI (rs-fMRI) data from the Human Connectome Project (HCP) database . We extracted local minima and optimal pathways among them in the energy landscape, and then explored the characteristics of the brain state transition processes from the network-theoretical perspective using the state transition network, where states are represented as nodes, transitions between two states (nodes) as edges, and transition rates as edge weights.
We analyzed the state transition network using the degree of state nodes (occurrence frequency during state transition) and path lengths (how many transient states are needed to arrive at a final state). From this transition network analysis, we also tested whether hub-like TSs exist, similar to spatial hubs found in the conventional network analysis of the brain [39–42]. We further examined state transitions toward the global minimum to identify intermediate stable states (a type of local minima) that mediate those multi-step transitions. The hierarchy in the brain state transition was investigated by clustering state transitions into intra- and inter-group transitions according to the transition rate. We finally investigated the organizational properties of the resting-state brain by comparing the transitional properties of the baseline cortical system with those of a virtual system by altering the MEM parameters.
The results of our analysis suggest that the cerebral cortex system at rest contains multiple stable states that are clustered into two major state groups. The transition between brain states across the two state groups was mediated by a frequent TS, which operated as a hub of the transition network. When we removed this hub state, which bridges most transition processes across the two groups, between-group transitions occurred via an alternative TS, indicating redundancy in state transition. State transition in the brain appears to involve multi-step state transitions, with some stable states serving as intermediate states for the complete transition. We also found that the baseline cerebral cortex at rest shows a more complex and organized state transition network than those of artificially altered systems. This network approach to the state transition in the brain may provide a new framework for brain exploration and become an effective tool for understanding healthy and abnormal brain systems, concerning brain state dynamics.
2. Materials and methods
2.1 Resting-state fMRI data set
The pairwise MEM (explained in the following section) was generated using rs-fMRI data of 470 participants (192 males, 278 females, age: 29.19 ± 3.51 years) from the HCP database , which was used in our previous study . Briefly, all data were sampled at TR = 0.72 s, during 4 runs, with 1200 time points per run. We projected the voxel-wise time courses from within a region on the eigenvector linked to the largest eigenvalue (obtained using principal component analysis) to generate the time series of interest for each region.
The effects of rigid motion and their derivatives were regressed out, followed by linear detrending and despiking of the extracted signals [43–46]. Although there is an ongoing debate concerning filtering and global regression, we regressed out global signal changes in the whole-brain mask to emphasize local and short-term state fluctuation (high or low) at each region in representing specific brain states, followed by a high-pass filtering at 0.01 Hz. Indeed, the current analysis is based on the assumption of the resting-state being in an equilibrium, i.e., without long-term statistical (temporal) changes in the dynamic properties.
Since computational cost dramatically increases with the number of degrees of freedom of the system (2N, N: number of nodes), we extracted the rs-fMRI time series of only 19 regions of interest (ROIs) out of 33 cortical regions defined in the Desikan-Killiani Atlas . In choosing 19 ROIs to define a cortical state, we included brain regions associated with the default mode network (DMN), the salience network, the parietal lobe and prefrontal cortices to focus on resting-state dynamics at the higher cognitive brain areas  in the cognitive hierarchy. Due to the complexity, we excluded primary/secondary sensory, motor cortical regions and subcortical brain regions in the evaluation. We also confined ROIs to a hemisphere since previous studies showed strong symmetric activities (e.g., symmetric independent components found in many previous studies, including Smith, Beckmann ) and strong interhemispheric connectivity; e.g., Kang, Pae . Coactive regions across hemispheres exhibit similar time courses and, thus, are considered to be less informative in defining diverse brain states.
The ROIs used in this study are the precuneus (PC), parahippocampal gyrus (PH), caudal middle frontal gyrus (cMF), fusiform gyrus (FG), inferior parietal lobe (IP), isthmus cingulate gyrus (IC), lateral orbitofrontal gyrus (lOF), medial orbitofrontal (mOF), pars-opercularis (Op), pars-orbitalis (Or), pars-triangularis (Tr), rostral anterior cingulate gyrus (AC), rostral middle frontal gyrus (rMF), superior frontal gyrus (SF), superior parietal gyrus (SP), supramarginal gyrus (SM), frontal pole (FP), temporal pole (TP), and insula (IN) in the left hemisphere (Fig 1A). The ROIs in the left-hemisphere were mainly evaluated and presented in the current study because the default mode network and many frontal cortices of interest have stronger connections (hubs) in the left hemisphere than the right hemisphere . However, we confirmed that similar results were obtained from the ROIs in the right-hemisphere (see Section B in S1 File).
For each ROI, signals were thresholded to represent deactive (0) and active (1) states. Since the number of local minima was maximized when the threshold was zero in the empirical evaluation (unpresented) and in our previous analysis , we selected zero as the threshold to binarize regional states after global regression (Fig 1B). A brain state was defined by merging all (binarized) 19 regional states into a state vector (the number of elements of a state vector is 19). Due to a high sample size demand to estimate brain states (for all 219 possible states), we concatenated all brain state samples from four sessions of 470 participants into a group-level sample data set (the total number of state samples, 1200 samples × 4 sessions × 470 participants) and estimated parameters of the group-level pairwise MEM using the method described in the following section.
2.2 Construction of pairwise maximum entropy model (MEM)
To analyze resting-state activity in the cerebral cortex, we utilized the pairwise MEM estimation approach described in previous studies (Fig 1B) [21–23].
The estimation process consists of a step for defining brain states and a pairwise MEM model for state dynamics, and an optimization step for MEM model parameters to fit probability distributions of empirical brain states and states generated by the model. Details of the model construction are provided in our previous study  (for the mathematical details, see review ref. ).
Suppose the brain state space is represented by
where the value of σi(i = 1, 2, …, N) is either 0 (deactive) or 1 (active), indicating a local activity at a node (brain region) i. N denotes the total number of nodes (or ROIs). For a probability distribution p(Vk) for all possible brain states Vk(k = 1, 2, …, 2N), the entropy S can be defined as
In the pairwise MEM, the average of each node activity,
and the averages of all pairwise products,
derived from the experimental data (i.e., σi(t), σj(t)), play as constraints that should be fit in the model estimation. With these constraints, maximizing the entropy S derives the probability p(Vk) of a brain state Vk as a Boltzmann distribution
where E(Vk) indicates the energy of the state Vk and can be described as the following equation:
where the parameters Hi and Jij represent the activation tendency (baseline sensitivity) of node i and the pairwise interaction between nodes i and j, respectively. σi(Vk) indicates an activity at a node i (i.e., σi) in a brain state Vk. A gradient ascent algorithm was employed to estimate MEM parameters, Hi and Jij. These parameters were iteratively updated to minimize the difference of average values in Eqs (3) and (4) between model-generated signals (with parameters) and observed signals using the following equations,
Here, 〈σi〉m and 〈σiσj〉m were calculated as follows:
The scale parameter ag was initially set to 0.1. The parameters were optimized until the gradients reached a value lower than 10−5.
To evaluate the effectiveness of the pairwise MEM, we calculated the accuracy value, rD,
Here, Dk is the Kullback-Leibler divergence between the probability distributions of the k-th order model network and the empirical network,
where pN represents the empirical distribution of the network state. In the calculation for the empirical probability distribution of brain states, we calculated the frequency of each state in the group-level sample data set described above.
We also evaluated the reliability parameter ER,
Here, rs and Sk are given by
The measures rD and rS evaluate the adequacy of the pairwise MEM over the independent MEM in explaining time series, in two different aspects; rD and rS use Kullback-Leibler divergences and difference of the entropy between independent (1st order) and pairwise (2nd order) MEMs. The reliability ER was defined to compare those two different measures. If Hi and Jij are estimated without error, ER is equal to 1.
2.3 Energy landscape analysis
To describe the dynamics of the cerebral cortex system at rest, we performed the energy landscape analysis. More specifically, first, we elucidated the local minima (attractors), and then evaluated energy barriers between pairs of attractors, following the procedure described in previous works [21, 23].
To construct an energy landscape, the distance between two states should be first defined. Based on this distance, neighbor states can be defined to extract local minima. Following previous energy landscape studies, we defined the distance between two states as the number of elements (bits) that differ between two state vectors. We also assumed a gradual state transition, and the energy landscape was examined by changing one element of the state vector for each step.
The local minima (also called stable states) were defined as states that have lower energy (more frequent) relative to their neighbors. To evaluate the energy barrier for each local minima pair, the lowest energy pathways were extracted by using disconnectivity graph analysis . Specifically, for each possible pair of local minima, we recorded the shortest path connecting the two local minima. The highest energy on this path was selected as a threshold to remove states that exhibited higher energy than the threshold. We repeated this step until the two local minima had been disconnected. The highest energy value of the last connected path was assigned to the threshold of the local minimum pair. The energy barrier, EB, between two local minima i and j, was defined as the lower value between Eth(Vi,Vj)–E(Vi) and Eth (Vi,Vj) − E(Vj), where Eth (Vi,Vj) represents the threshold as defined above. These disconnectivity graph calculations were performed using the i-graph library .
2.4 Construction of the state transition networks
In the present study, we constructed three types of state transition networks; a full state transition network composed of all possible states along transition pathways among local minima (STN-FS), a state transition network of states from local minima toward the global minimum (STN-GM), and a state transition network among TSs and local minima states (STN-LM).
We first constructed a STN-FS as a directional weighted network (Fig 2A). For all possible pairs of local minima, state transition pathways were identified as described in the above section. All states on the state transition pathways among local minima were considered nodes of the STN-FS. Since the forward and backward state transition pathways were identical for a pair of local minima, we only considered the state transitions from the higher to lower local minimum. For all edges (transitions between pairs of states), we assigned weights using a transition rate, which is given by
Here, EB represents the energy barrier (or activation energy) which is the energy difference between the highest energy on the transition path and the energy of the initial state (see Fig 1C).
We then constructed the STN-GM to focus on the details of the state transition processes toward the global minimum (Fig 2B). To construct this network, all nodes and edges that were not connected to the global minimum in the STN-FS were removed. In the STN-GM analysis, we particularly identified intermediate states that mediate multi-step transitions toward the global minimum.
In order to focus on the transition rate among local minima, we constructed a STN-LM (Fig 3A), considering local minima and TSs as nodes and transition rates as edges. According to the transition state theory, a transition rate between two states is determined solely by the energy difference of the initial state and the rate-determining transition state (saddle point), i.e., the TS (Fig 1C). In other words, the transition rate between two states only depended on its energy barrier (maximal energy difference) along the pathway.
Based on the transition rates in STN-LM, we evaluated the architecture of the state transitions of the system by applying a cluster analysis of transition rates to differentiate intra- and inter-group transitions. The clustering of transitions in the STN-LM was conducted using the UPGMA (Unweighted Pair Group Method with Arithmetic Mean) algorithm  with Euclidian metric; distances between nodes were defined by transition rates.
The constructed network was analyzed in terms of network theory using node degree (frequency of occurrence during state transition) and path length (number of transitions to reach the target state). We also analyzed the effective path length between local minima, which is defined as the difference between the total path length (number of transitions) and the Hamming distance (number of different elements between two state vectors, i.e., number of regions that showed active/deactive differences) of the initial and final states.
We finally investigated the organizational properties of the resting-state brain by comparing the transitional properties of the baseline cortical system with those of a virtual system by altering the MEM parameters. The comparison between the resting-state brain and the altered (virtual) system was conducted over the STN-FS and STN-LM.
3.1 Maximum entropy model for the cerebral cortex system at rest
In order to generate the energy landscape of the brain state, we estimated the first and second order interaction parameters (i.e., baseline sensitivity Hi and pairwise interactions Jij) of the MEM using binarized rs-fMRI activation patterns. The activation patterns of rs-fMRI data were reproduced with a high accuracy of fit (rD = 86.3%) and reliability (ER = 99.9%) (Figure A in S1 File). Baseline sensitivity parameters Hi and pairwise interaction, Jij, are displayed in Figure A in S1 File and Fig 1B. Details of the obtained MEM parameters are described in Figure A in S1 File.
3.2 Multiple stable states in the resting-state of the cerebral cortex system
Analysis of the energy landscape identified 14 local minima of the cerebral cortex system at rest (Figs 2 and 3). Complementary states (active versus deactive for each brain region) of five local minima were also found to be local minima (Fig 3). Two pairs, local minimum (LM) 1 and LM12, LM7 and LM8, were nearly complementary states of each other. In these pairs, all regions were complementary except for one brain region, in each: the deactive precuneus (PC) in the LM7 and LM8, and deactive fusiform gyrus (FG) in the LM1 and LM12.
The most stable local minimum (i.e., global minimum) was LM11, where most cortical regions were active except for the insula, supramarginal gyrus, superior parietal lobe, and fusiform gyrus (see LM11 map in Fig 2C).
3.3 The state transition network among full states (STN-FS)
Utilizing disconnectivity graph analysis , 91 transition pathways were extracted for all possible pairs of the 14 local minima. All states on the 91 transition pathways were regarded as nodes and transition rates between pairs of nodes as edges in the STN-FS (Fig 1C and 1D). As a result, a total of 219 nodes and 1201 edges composed an STN-FS (Fig 2A). When we evaluated node degrees for all nodes in the STN-FS, three (state) nodes showed a significantly higher node degree than the rest (arrows in Fig 2A). Most (86.8%) effective path lengths, i.e., the difference between the total path length and the Hamming distance of two initial and final local minima, were less than 8 (Figure B in S1 File). For half of the total number of pathways (49 pathways), effective path lengths had the shortest value of 0.
3.4 Analysis of state transition network from local minima toward the global minimum (STN-GM)
For 14 local minima, 13 transition processes toward the global minimum were considered in the STN-GM with 82 nodes and 141 edges. A total of eight TSs, which determine transition rate, were found in the STN-GM. The STN-GM also showed that intermediate local minima (e.g., LM7, LM9, LM12, and LM14) were involved in the transition processes of other local minima transitioning toward the global minimum (Fig 2B and 2C). For instance, the transition pathways that started from the LM1, LM3, LM4, LM6, and LM12 passed through LM7 before reaching the global minimum (LM11). The rates for these transitions were determined by energy differences between TS2 and the initial local minima. The state transition from LM6 to the global minimum (LM11) contained an intermediate state (LM7) and the energy of the rate-determining transition state between LM6 and LM7 was smaller than that of LM7 and LM11 (i.e., energy of TS2), and, thus, the rate-determining transition state was TS2 (upper Fig 2D).
However, LM5, LM8, LM9, LM10, LM12, LM13, and LM14 had their own rate-determining transition states along transition paths toward the global minimum (Fig 2B and 2C). Indeed, in the state transition from LM14 to the global minimum (LM11), which contained an intermediate state (LM12), the energy of the rate-determining transition state between LM11 and LM12 (energy of TS8) was larger than that of LM12 and LM11 (energy of TS6), and, thus, the rate-determining transition state was TS8 (lower Fig 2D).
In this way, by analyzing this reduced state transition network (STN-GM), we could identify the characteristics of all transition processes on their way to the global minimum.
3.5 Analysis of a state transition network among rate-determining transition states and local minima states (STN-LM)
The STN-LM was composed of 27 nodes (13 TSs plus 14 local minima) and 90 edges (Fig 3A). We found a clustered structure in the STN-LM: one cluster containing six local minima (LM8, LM9, LM11, LM12, LM13, and LM14), and the other containing their complementary local minima (LM1, LM2, LM3, LM4, LM6, and LM7). A similar clustering result was found by using energy barriers (Fig 3B). Interestingly, only one rate-determining transition state, TS2, bridged two clusters. TS2 is composed of active regions in the FP, SF, AC, mOF, PC, IP, IC, TP, PH, FG, which overlap mostly with coactivation of the default mode network  and the anterior and medial temporal lobe (Fig 3D). Since TS2 has a high node degree (Fig 2A), we can refer to TS2 as a hub in the transition network.
To investigate the effects of TS2 on the transition process, we removed TS2 and explored the state transition pathways. After the removal of TS2, we found 36 alternative pathways. In these 36 pathways, the complementary state of TS2, namely ~TS2, bridged the state transition processes between clusters, instead of TS2. The energy difference between TS2 and ~TS2 was very small, 0.00737. Since the transition rate is proportional to exp(-Ebarrier) in state transition theory, the estimated ratio of the transition rate between the original and alternative pathway was 99.27%. Thus, we added the ~TS2 node to the state transition network, STN-LM. Moreover, the property of clustered transitions was also observed for the transition processes in the ~TS2 system (Figure C in S1 File). Since TS2 and ~TS2 have high node degrees (i.e., measure of frequency of appearance), both TS2 and ~TS2 play as hub TSs.
We extracted the factors that determined transition rates (i.e., energy barriers) for both the TS2 system and TS2 removed (~TS2) system (Fig 4); the Hamming distances between initial and final states were positively correlated with energy barriers (r = 0.608, p = 1.654 × 10−10 for the TS2 system, and r = 0.611 p = 1.276 ×10−10 for the ~TS2 system, Fig 4B). However, there was no such relation between the energy barriers and effective path lengths (Figure C in S1 File).
We further investigated the transition processes by separating the inter-group and intra-group processes (Fig 4D, 4E and 4F). For the intra-group transitions, positive correlations between Hamming distances and path lengths were observed for both TS2 and ~TS2 systems.
However, for the inter-group transitions, we could not find such associations. Thus, we further separated the inter-group transitions and found negative and positive correlations between Hamming distances and path lengths for the transitions from group 1 to 2 (r = -0.865, p = 3.123 × 10−5), and from group 2 to 1 (r = 0.921, p = 3.240 × 10−9) in the TS2 system (Fig 4D).
Interestingly, in the ~TS2 system, the correlations were reversed; positive and negative correlations were found for the transitions from group 1 to 2 (r = 0.865, p = 3.123 × 10−5), and from group 2 to 1 (r = -0.921, p = 3.240 × 10−9) (Fig 4D). Since the distances between TS2 and local minima of group 1 were longer than those of group 2 in the TS2 system, but in the ~TS2 system, an inverse association was found (i.e., the distances between ~TS2 and local minima of group 1 were shorter than those of group 2) (Fig 4C), the cause of the inverse correlations could be the distance between the transition state and initial states. In the inter-group transitions, path lengths from an initial local minimum to the other group local minima depended on Hamming distances from the initial local minimum to the hub TS (Fig 4E and 4F).
3.6 Effects on the state transition processes following the alteration of the system
The effects of the global strengths of pairwise interactions on the transition process were investigated by scaling all Jij parameters (αJij, α = 0.0, 0.1, …, and 5.0). Both increases and decreases of the scale factor α tended to reduce the total number of local minima (Fig 5A). When a markedly small or large-scale factor, α < 0.7 or α > 1.7, was used, only one local minimum was found.
Here, we further analyzed two representative examples: α = 0.8 and α = 1.2. In both cases, the total number of local minima was reduced by perturbations, from 14 (baseline resting-state) to 7 (for α = 0.8) and 9 (for α = 1.2), and thus 21 and 36 transition processes were considered in their state transition networks, respectively.
In the state transition network of the weak coupling system (α = 0.8), 137 nodes and 258 edges, which were decreased compared to the baseline resting-state, were found. Under the strong pairwise interaction (α = 1.2), compared to the baseline resting-state, the total number of nodes was decreased (225 nodes) and that of edges was increased (456 edges). Since most pairwise parameters were positive, the energy of the states increased and decreased for weak and strong alterations in the scale parameter, respectively. In both cases, a positive correlation between node degree and energy was found for the nodes of the transient and transition states (Figure D in S1 File).
In the baseline resting-state, nodes were densely connected to other nodes (Figure D in S1 File), and, the maximum node degree was 44, which was larger than that of the altered systems (8 and 10, for α = 0.8, and 1.2, respectively). For all cases, more than half of the effective path lengths had a value of 1 (Fig 4). In the weak coupling system (α = 0.8), the longest effective path length was 15, which was smaller than that of the others (21 and 25 for α = 1.0, and 1.2, respectively).
In contrast to the baseline resting-state, in these altered systems, simple and deep energy valleys were found (Fig 5C). Indeed, the state transition processes were simpler than those of the baseline resting-state (Fig 5D). Except for LM6 in the weak coupling system, all local minima directly transitioned to their global minimum (Figure D in S1 File).
The brain at rest has been considered a highly dynamic complex system operating at a critical value of coupling that maximizes multistability [21, 34–36, 55]. Beyond the multistability of the resting-state cortical system, we systematically investigated the architecture of the state transition processes by applying a graph-theoretical analysis to state transition. State transition network analysis suggests a well-organized state transition process embedded in the resting-state human cerebral cortex system. The characteristics of the state transition in the resting-state cortex system are discussed in the subsequent paragraphs.
4.1 Intermediate states in brain state dynamics
The resting-state brain has intermediate states in state dynamics. Some state transitions in the cerebral cortex system toward the global local minimum occurred in multi-steps via several intermediate stable states (or intermediate local minima) (Figs 2C and 6B). This phenomenon is similarly found in biochemical reactions by enzymes in biological systems [30, 56–59]. For instance, during the rebinding of ligand (CO or O2 molecules) to myoglobin after photolysis, several intermediate states were observed in spectroscopic experiments, and these intermediates have often been explained in terms of the regulation of ligand binding mechanisms [57–59]. The state transitions during membrane fusion processes occur via intermediate states, which were explored in computational and experimental studies [60–62]. Similar to these phenomena in biochemical systems, the intermediate stable states during brain state transition may also play a role in lowering energy barriers. We speculate that this lowered energy barrier may regulate and expedite transitions along certain pathways of brain state transitions in the resting-state whole-brain system. It should also be noted that transitions between some pairs of local minima are straightforward without any intermediate transition states.
4.2 Characteristics of brain state transitions in the resting-state dynamics
Current network analysis of stable states of the cerebral cortex suggests several characteristics of brain state dynamics.
First, in the brain state transition network, local minima are highly clustered mainly into two groups, and the manner in which state transitions among distributed local minimum occurred was different between inter-group and intra-group transitions.
Second, most inter-group state transitions (from a local minimum at a cluster to a local minimum at the other) occurred via some hub transition states (saddle points) (e.g., TS2 and ~TS2) in the transition pathway (Fig 6D). This phenomenon makes the inter-group transition different from the intra-group state transition, where the transition state along the path between two states differed according to the initial state. Those hub transition states are analogous to hubs found in the conventional network analysis of the brain connectome. Brain connectome studies have shown a hub-like structural architecture in the brain, which is considered to mediate efficient information exchange [39–42]. Similar to network analysis focusing on the spatial geometry of the connectome, the current result suggests that inter-group brain state transitions occur mostly via hub states (more frequently occurring states) in the temporal geometry (Fig 6C).
Third, path lengths (number of transitions to reach the target state) were positively correlated with the Hamming distances (i.e., number of different state bits (or regions)) between the initial and final states within the intra-group transitions. This result implies that the transition states of intra-group transitions take advantage of shorter transition paths, i.e., efficient transition from a brain state to the other state with minimal transition numbers. In the inter-group state transitions, path lengths were determined by the Hamming distance between the initial state and the hub transition state on the path, not the final state (Figs 4D, 4E, 4F and 6D). We speculated that for significant state changes in the cortical system, the brain may minimize transitioning costs by traveling via hub transition states, not simply following short transition paths.
4.3 Redundant inter-group state transition in the cortical system
The current simulation study suggests that the cerebral cortex has redundant transition pathways. A transition state (e.g., TS2) mediates most of the inter-group state transition processes, serving as a transition hub in the resting-state transition network (STN-LM). When we excluded this hub transition state, its complementary state ~TS2 appeared to serve as a detour for inter-group transitions with similar transition rates. This alternative hub ~TS2 was the complementary state of the hub state (TS2) of the baseline resting-state. The energy level of ~TS2 was similar to that of TS2 and the rates of the transition processes were similar to each other. We considered pathways via ~TS2 as “redundant” pathway in inter-group transition processes, as those were chosen after removing TS2 as a transition state.
The existence of multiple pathways has been reported in some complex systems where a reaction occurs among multiple units cooperatively. For example, in some biomolecule systems, e.g. F1-ATPase and myosin V motor, two reaction pathways have been reported [63, 64]. Multiple transition pathways may be associated with “degeneracy” or “redundancy” in the complex brain system [65, 66]. The redundant pathways could be particularly advantageous in maintaining effective state transitions when a certain transient state cannot play its role in the state transition.
4.4 Baseline brain state dynamics compared to those of altered systems
To understand the organization principle of the baseline configuration for brain state dynamics as done in Kang et al. , we compared the baseline (observed) network configuration with virtual (artificially altered) network configurations by scaling the pairwise interactions, and analyzed their state transition dynamics. Compared to the virtual networks, the baseline resting-state transition network (STN-FS) contained a bigger number of states with high node degrees and relatively longer path lengths. Furthermore, neither the clustered structure nor the intermediate states of the baseline system were observed in the altered virtual systems. For example, in the STN-LM of the baseline system, which contains nodes on the pathway toward the global minimum (LM11) from other local minima (Fig 2), four types of transition processes were identified with four local minima (LM7, LM9, LM12, and LM14) as intermediate states. This property was not found in the virtual systems, which showed a much simpler transition network (Fig 2 and S3). We also conducted energy landscape analysis of network configurations with random values and with randomly permutated Hi and Jij. The results presented in the Supporting Information suggest simpler energy landscape and simpler transition network in the random network configurations than in the real resting-state network (see Section C in S1 File).
Furthermore, hub intermediate states and hub transition states were found only in the baseline resting-state system, not in the virtual systems. In the virtual systems, altered by scaling pairwise interactions from the baseline resting-state system, local minima having a high node degree disappeared; the hub-like local minima of the baseline system (having relatively high energy) were eliminated first by scaling pairwise interactions. This phenomenon, of which higher energy was eliminated first by network alteration, was consistent with that reported in our previous study on the subcortical system . Meanwhile, low-energy local minima tended to persist even after network alteration. If we consider some task performances as states deviating from the baseline system, those sustaining local minima may act as common bases from which diverse functions arise or as fundamental elements of maintenance of dynamic brain systems. It may be that the baseline cerebral cortex system is configured to allow network systems to effectively transition among diverse brain states, which may be a necessary element in the workings of the complex network systems exhibiting multiple stable states.
All these transition characteristics may possibly be embedded in the nonlinear coupling over the structural network. We speculate that network topology may provide a biased playground of multistability, and endogenous fluctuations during resting-state may drive state transitions over the structural playground. This interpretation about the interplay between network topology and noise is in line with the dynamic nature of the brain [9, 21]. In the current study, we showed that the multistable nature of brain states and the well-organized properties of the transition processes can emerge from nonlinear interactions over the cortical brain network.
4.5 Nonlinear dynamic brain systems
Resting-state brain dynamics and non-stationary functional connectivity have recently been explored by evaluating brain connectivity in sliding window fashion from the viewpoint of a linear system [11–20]. For example, Park and colleagues assumed linear interactions (connectivity) between brain regions change by time, and these time-dependent interactions were estimated for consecutive windows . In contrast, multistability in the complex system is an emergent property of nonlinear interactions among nodes in the system, without any changes in the internal connectivity. From the perspective of the nonlinear system, Hansen, Battaglia  suggested that non-stationary functional connectivity (FC) (particularly, rapid transitions switching between a few discrete FC states) can be explained by the non-linearity of the nodal activity that derives the structural brain system. Spiegler, Hansen  attributed the nonstationary FC to the criticality of the nonlinear brain system embedded in the structural network topology. Similarly, Cabral et al. also showed that dynamic functional connectivity can emerge from a static structural connectivity with various non-linear dynamic models of the brain . They showed that diverse FC states are emergent when the brain is operating at the edge of criticality. Pillai and Jirsa  also showed that multiple sub-states undergo structured flows on the manifold of the low-dimensional state spaces (functional subspaces) and this emergent behavior is attributable to the synaptic coupling level over the nonlinear interactions. Rabinovich and colleagues  argues that both flexible and reproducible transitions among multiple meta-stable states can emerge in the nonlinear system, which may explain state transitions in the decision-making process. All those studies [37, 55, 67–69] are based on a model with nonlinear temporal dynamics described using a differential equation. In terms of the nonlinear interaction and its consequent emergence of multiple states in the complex brain, our approach using pairwise MEM is in line with those studies [37, 55, 67–69]. However, the current evaluation differs from those studies in that the analysis of brain dynamics with a pairwise MEM is based on the statistical mechanics, which deals with the emergence of stable states and their transitions in terms of probability. Nevertheless, the two approaches (analysis with a differential equation and analysis in the statistical mechanics) are known to be equivalent since ensemble probability distributions of each state (in statistical mechanics) can be derived from a large number of trajectories as solutions of a differential equation (for a rigorous explanation, see ergodic theorem in the statistical physics).
A small number of major transition paths in the current study may serve as manifold-like transitions found in Pillai and Jirsa , where a lower dimensional manifold of state transitions was induced by an asymmetric interaction (due to task). The organized state transitions explored in this study may be correspondent to reproducible transitions among multiple meta-states in Rabinovich et al. .
4.6. Neurobiological interpretation
In this study, the transition state TS2 appears to mediate most transitions, particularly from inter-group transitions between complementary states. Neurobiologically, TS2 (Fig 7) corresponds to the typical default mode network as shown in [70, 71]. In particular, the DMN in  has a distributed pattern of activations at the posterior cingulate cortex, medial prefrontal cortex, anterior cingulate cortex, superior frontal cortex, frontal pole and inferior parietal lobe. Many of those DMN regions are known to play hubs in the structural brain network . The current transition network analysis newly suggests that DMN also serves as a transition hub in the state transition. In the nonlinear systems perspective, the brain has been considered to be in a saddle point of phase transition, ready to perform diverse tasks rather than in a stable state [21, 34–36, 55]. In this respect, it is not surprising that DMN, a core of the brain network, mediates most state transitions.
The global minimum LM11 represents a combined activation of TS2 (DMN) and the middle and inferior prefrontal cortices. Why the coactivation at the prefrontal and default mode networks appears frequently is an open question but some previous works support this phenomenon. The distributed metabolic activity pattern in a group average fluorodeoxyglucose (FDG) positron emission tomography (PET) in resting state  is highly consistent with the pattern of LM11. Considering that FDG-PET indicates accumulated energy metabolism of the brain at rest, increased metabolism at brain regions corresponding to LM11 explains frequent occurrence of LM11, i.e., global minimum. The coactivation among distributed regions of the DMN and prefrontal cortices may be attributable to anatomical neural projections. For example, the regions of DMN are highly interconnected via anatomical projections observed in fiber tracking of diffusion tensor imaging .
Note that TS2 (i.e., DMN) mediates sequential state transitions particularly between complementary states (in polarity). When brain states in Group 2 (Fig 3) transit to their complementary states in Group 1 or eventually to the global minimum state (LM11), those states progressively (but statistically) transit to destinations via the DMN transition state (TS2). This polarity switching has also been observed in previous MEM studies [25, 26, 75]. All this transition occurs progressively, via transition states or intermediate local minima. The progressive spatio-temporal transition may be associated with appearance of microstate modes detected in EEG  and some phenomena of travelling waves [77, 78] or metastable waves . Considering (semi-)cyclic waves of spatial activation patterns during resting state [76, 78, 79], the progressive transition between complementary states or global minimum is not surprising.
The local minima found in the current study may be associated with diverse introspective thoughts or unconscious reactions to intrinsic sensations and the state transitions among local minima may reflect mind-wandering. Although the functional role of each local minimum and the mechanism for state transition are unknown, the anatomical configuration of the brain network may contribute to form the typical energy landscape of brain states, which induces organized transitions in this study. All these speculations need to be further researched.
4.7 Limitations and challenges
The current study has several limitations and challenges. Due to the high computational cost and the requirement of a large sample size, we evaluated the dynamics of a reduced brain system at each hemisphere (focusing on the left hemisphere in the main text) but did not evaluate those of the whole-brain system. In spite of strong symmetry (e.g., 21), the interaction between two hemispheres and its effect on the dynamic system in the whole-brain system remain to be explored. For the same reason, we did not include sensory and subcortical brain regions in the current analysis. We speculate that the dynamic properties of the whole-brain nonlinear system would be much more complex. In spite of technical challenges, exploration of the state transition properties of the whole-brain system with more precisely parcellated brain regions would greatly expand our understanding of the brain system.
In the preprocessing step, we decided to conduct a global regression since we focus on the dynamics of regional brain states induced by pairwise interactions within the system (not induced by any other modulatory effects) in an equilibrium. To define a regional brain state with an activity of a region that fluctuates locally, we applied MEM to fMRI signals after global regression. The global regression emphasized properties of the state dynamics when it was applied to the subcortical brain system in our previous study .
Like any other modeling approach, the modelling with MEM simplifies the dynamic phenomenon of the brain by binarizing regional activities into two binarized states, “1” and “0”. This simplification may not be precise enough to describe microlevel states at a region but this simplification enables us to explore nonlinear spatio-temporal dynamics of the complex brain system, which may otherwise not be easily achieved. Since binarization of regional signals was done after global regression, the regional states “1” and “0” are defined relative to the global mean at each time point. Thus, we denoted “0” as deactive (suppressed) compared to the global mean although it may denote inactivation at some regions. Nevertheless, the neurophysiological underpinnings are not clear whether the low fMRI amplitudes (denoted as “0”) indicate inactive (silent) or deactive (suppressed) neural states. Due to the complex neuro-vascular coupling in the BOLD-fMRI , even high fMRI amplitudes (denoted as “1”) can possibly be induced by increased metabolic demands after firing of inhibitory populations (a suppressed brain state). Therefore, the definition of brain states should be further clarified with more pieces of knowledge on the neural states reflected in the fMRI signals.
Despite this simplification (binarization), the pairwise MEM provides a new framework of brain research by considering the distributed brain activity as a result of nonlinear pairwise interactions among brain regions in a system. Some computational models have nonlinear formulations of interactions comparable to the pairwise MEM. Deco and colleagues  introduced a nonlinear dynamic brain model based on a dynamical mean field reduction. The method in , however, does not estimate interaction parameters directly. Instead, they simplified the model by adopting structural connectivity as interaction parameters . In the task-based fMRI analysis, a nonlinear dynamic causal model has been proposed . The nonlinear dynamic causal model , which is practically limited to a small size circuit with an external perturbation (input), does not provide a direct way of evaluating the probability of state occurrence during resting state.
In contrast to these nonlinear approaches to the brain, Vidaurre and colleagues [83, 84] have introduced a hidden Markov modelling (HMM) to identify brain states and to explore state transitions in the brain. However, a state in the HMM is defined in terms of a whole-brain connectivity pattern (embedded in the cross-spectral density, corresponding to the system parameter Jij in the MEM) that changes along the time course. Similarly, the dynamic nature of the brain has been explored in terms of temporal changes in its functional connectivity [11–20]. Meanwhile, MEM assumes a nonlinear system that generates multistablity without changes in the system configuration and the state transition is a result of nonlinear interactions. In this respect, MEM is a practical way of exploring the brain’s complex dynamics, which is not easily replaced. In line with this notion, many studies have shown the usefulness of MEM in the analysis of brain dynamics with respect to a nonlinear system [21, 25, 26, 75, 85].
Previous studies have revealed alterations in the dynamics of networks associated with brain disorders such as schizophrenia [86, 87] and autism . A growing number of studies are showing altered dynamics in other brain disorders as well. However, the dynamic properties in brain diseases have not been thoroughly researched. The current frameworks for dynamic brain states can be used to identify altered dynamic architectures in neuropsychiatric disorders. Research using clinical data will yield results that can validate the usefulness of the proposed approach.
1. Deco G, Tononi G, Boly M, Kringelbach ML. Rethinking segregation and integration: contributions of whole-brain modelling. Nature reviews Neuroscience. 2015;16(7):430–9. doi: 10.1038/nrn3963 26081790
2. Deco G, Jirsa VK. Ongoing cortical activity at rest: criticality, multistability, and ghost attractors. The Journal of neuroscience: the official journal of the Society for Neuroscience. 2012;32(10):3366–75.
3. Cabral J, Kringelbach ML, Deco G. Exploring the network dynamics underlying brain activity during rest. Progress in neurobiology. 2014;114:102–31. doi: 10.1016/j.pneurobio.2013.12.005 24389385
4. Tognoli E, Kelso JA. The metastable brain. Neuron. 2014;81(1):35–48. doi: 10.1016/j.neuron.2013.12.022 24411730
5. Freyer F, Roberts JA, Ritter P, Breakspear M. A canonical model of multistability and scale-invariance in biological systems. PLoS Comput Biol. 2012;8(8):e1002634. doi: 10.1371/journal.pcbi.1002634 22912567
6. Freyer F, Roberts JA, Becker R, Robinson PA, Ritter P, Breakspear M. Biophysical mechanisms of multistability in resting-state cortical rhythms. The Journal of neuroscience: the official journal of the Society for Neuroscience. 2011;31(17):6353–61.
7. Kelso JA. Multistability and metastability: understanding dynamic coordination in the brain. Philos Trans R Soc Lond B Biol Sci. 2012;367(1591):906–18. doi: 10.1098/rstb.2011.0351 22371613
8. Schwartz JL, Grimault N, Hupe JM, Moore BC, Pressnitzer D. Multistability in perception: binding sensory modalities, an overview. Philos Trans R Soc Lond B Biol Sci. 2012;367(1591):896–905. doi: 10.1098/rstb.2011.0254 22371612
9. Breakspear M. Dynamic models of large-scale brain activity. Nature neuroscience. 2017;20(3):340–52. doi: 10.1038/nn.4497 28230845
10. Rabinovich MI, Varona P. Robust transient dynamics and brain functions. Front Comput Neurosci. 2011;5:24. doi: 10.3389/fncom.2011.00024 21716642
11. Monti RP, Hellyer P, Sharp D, Leech R, Anagnostopoulos C, Montana G. Estimating time-varying brain connectivity networks from functional MRI time series. NeuroImage. 2014;103:427–43. doi: 10.1016/j.neuroimage.2014.07.033 25107854
12. Park H-J, Friston KJ, Pae C, Park B, Razi A. Dynamic effective connectivity in resting state fMRI. NeuroImage. 2017.
13. Jeong SO, Pae C, Park HJ. Connectivity-based change point detection for large-size functional networks. NeuroImage. 2016;143:353–63. doi: 10.1016/j.neuroimage.2016.09.019 27622394
14. Hutchison RM, Womelsdorf T, Gati JS, Everling S, Menon RS. Resting-state networks show dynamic functional connectivity in awake humans and anesthetized macaques. Hum Brain Mapp. 2013;34(9):2154–77. doi: 10.1002/hbm.22058 22438275
15. Handwerker DA, Roopchansingh V, Gonzalez-Castillo J, Bandettini PA. Periodic changes in fMRI connectivity. NeuroImage. 2012;63(3):1712–9. doi: 10.1016/j.neuroimage.2012.06.078 22796990
16. Chang C, Glover GH. Time-frequency dynamics of resting-state brain connectivity measured with fMRI. NeuroImage. 2010;50(1):81–98. doi: 10.1016/j.neuroimage.2009.12.011 20006716
17. Allen EA, Damaraju E, Plis SM, Erhardt EB, Eichele T, Calhoun VD. Tracking whole-brain connectivity dynamics in the resting state. Cereb Cortex. 2014;24(3):663–76. doi: 10.1093/cercor/bhs352 23146964
19. Calhoun VD, Miller R, Pearlson G, Adali T. The chronnectome: time-varying connectivity networks as the next frontier in fMRI data discovery. Neuron. 2014;84(2):262–74. doi: 10.1016/j.neuron.2014.10.015 25374354
20. Preti MG, Bolton TA, Van De Ville D. The dynamic functional connectome: State-of-the-art and perspectives. NeuroImage. 2017;160:41–54. doi: 10.1016/j.neuroimage.2016.12.061 28034766
21. Kang J, Pae C, Park HJ. Energy landscape analysis of the subcortical brain network unravels system properties beneath resting state dynamics. Neuroimage. 2017;149:153–64. doi: 10.1016/j.neuroimage.2017.01.075 28159684
22. Watanabe T, Hirose S, Wada H, Imai Y, Machida T, Shirouzu I, et al. A pairwise maximum entropy model accurately describes resting-state human brain networks. Nat Commun. 2013;4:1370. doi: 10.1038/ncomms2388 23340410
23. Watanabe T, Hirose S, Wada H, Imai Y, Machida T, Shirouzu I, et al. Energy landscapes of resting-state brain networks. Frontiers in neuroinformatics. 2014;8:12. doi: 10.3389/fninf.2014.00012 24611044
24. Watanabe T, Kan S, Koike T, Misaki M, Konishi S, Miyauchi S, et al. Network-dependent modulation of brain activity during sleep. NeuroImage. 2014;98:1–10. doi: 10.1016/j.neuroimage.2014.04.079 24814208
25. Watanabe T, Masuda N, Megumi F, Kanai R, Rees G. Energy landscape and dynamics of brain activity during human bistable perception. Nat Commun. 2014;5:4765. doi: 10.1038/ncomms5765 25163855
26. Ezaki T, Sakaki M, Watanabe T, Masuda N. Age-related changes in the ease of dynamical transitions in human brain activity. Hum Brain Mapp. 2018;39(6):2673–88. doi: 10.1002/hbm.24033 29524289
27. Gu S, Cieslak M, Baird B, Muldoon SF, Grafton ST, Pasqualetti F, et al. The Energy Landscape of Neurophysiological Activity Implicit in Brain Network Structure. Sci Rep. 2018;8(1):2507. doi: 10.1038/s41598-018-20123-8 29410486
28. Rao F, Karplus M. Protein dynamics investigated by inherent structure analysis. Proceedings of the National Academy of Sciences of the United States of America. 2010;107(20):9152–7. doi: 10.1073/pnas.0915087107 20435910
29. Li CB, Yang H, Komatsuzaki T. Multiscale complex network of protein conformational fluctuations in single-molecule time series. Proceedings of the National Academy of Sciences of the United States of America. 2008;105(2):536–41. doi: 10.1073/pnas.0707378105 18178627
30. Frauenfelder H, Sligar SG, Wolynes PG. The energy landscapes and motions of proteins. Science. 1991;254(5038):1598–603. doi: 10.1126/science.1749933 1749933
31. Gfeller D, De Los Rios P, Caflisch A, Rao F. Complex network analysis of free-energy landscapes. Proceedings of the National Academy of Sciences of the United States of America. 2007;104(6):1817–22. doi: 10.1073/pnas.0608099104 17267610
32. Delvenne JC, Yaliraki SN, Barahona M. Stability of graph communities across time scales. Proceedings of the National Academy of Sciences of the United States of America. 2010;107(29):12755–60. doi: 10.1073/pnas.0903215107 20615936
33. Goldstein M. Viscous Liquids and the Glass Transition: A Potential Energy Barrier Picture. The Journal of Chemical Physics. 1969;51(9):3728–39.
34. Golos M, Jirsa V, Dauce E. Multistability in Large Scale Models of Brain Activity. PLoS Comput Biol. 2015;11(12):e1004644. doi: 10.1371/journal.pcbi.1004644 26709852
35. Deco G, Senden M, Jirsa V. How anatomy shapes dynamics: a semi-analytical study of the brain at rest by a simple spin model. Front Comput Neurosci. 2012;6:68. doi: 10.3389/fncom.2012.00068 23024632
36. Deco G, Jirsa VK, McIntosh AR. Resting brains never rest: computational insights into potential cognitive architectures. Trends in neurosciences. 2013;36(5):268–74. doi: 10.1016/j.tins.2013.03.001 23561718
37. Hansen EC, Battaglia D, Spiegler A, Deco G, Jirsa VK. Functional connectivity dynamics: modeling the switching behavior of the resting state. NeuroImage. 2015;105:525–35. doi: 10.1016/j.neuroimage.2014.11.001 25462790
38. Van Essen DC, Ugurbil K, Auerbach E, Barch D, Behrens TE, Bucholz R, et al. The Human Connectome Project: a data acquisition perspective. NeuroImage. 2012;62(4):2222–31. doi: 10.1016/j.neuroimage.2012.02.018 22366334
40. Senden M, Reuter N, van den Heuvel MP, Goebel R, Deco G. Cortical rich club regions can organize state-dependent functional network formation by engaging in oscillatory behavior. NeuroImage. 2016.
41. van den Heuvel MP, Sporns O. Rich-club organization of the human connectome. The Journal of neuroscience: the official journal of the Society for Neuroscience. 2011;31(44):15775–86.
42. Honey CJ, Sporns O. Dynamical consequences of lesions in cortical networks. Hum Brain Mapp. 2008;29(7):802–9. doi: 10.1002/hbm.20579 18438885
43. Weissenbacher A, Kasess C, Gerstl F, Lanzenberger R, Moser E, Windischberger C. Correlations and anticorrelations in resting-state functional connectivity MRI: a quantitative comparison of preprocessing strategies. NeuroImage. 2009;47(4):1408–16. doi: 10.1016/j.neuroimage.2009.05.005 19442749
44. Power JD, Barnes KA, Snyder AZ, Schlaggar BL, Petersen SE. Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion. NeuroImage. 2012;59(3):2142–54. doi: 10.1016/j.neuroimage.2011.10.018 22019881
45. Thomas JB, Brier MR, Bateman RJ, Snyder AZ, Benzinger TL, Xiong C, et al. Functional connectivity in autosomal dominant and late-onset Alzheimer disease. JAMA neurology. 2014;71(9):1111–22. doi: 10.1001/jamaneurol.2014.1654 25069482
46. Taylor JS, Rastle K, Davis MH. Interpreting response time effects in functional imaging studies. NeuroImage. 2014;99:419–33. doi: 10.1016/j.neuroimage.2014.05.073 24904992
47. Desikan RS, Segonne F, Fischl B, Quinn BT, Dickerson BC, Blacker D, et al. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. NeuroImage. 2006;31(3):968–80. doi: 10.1016/j.neuroimage.2006.01.021 16530430
48. Power JD, Cohen AL, Nelson SM, Wig GS, Barnes KA, Church JA, et al. Functional network organization of the human brain. Neuron. 2011;72(4):665–78. doi: 10.1016/j.neuron.2011.09.006 22099467
49. Smith SM, Beckmann CF, Andersson J, Auerbach EJ, Bijsterbosch J, Douaud G, et al. Resting-state fMRI in the Human Connectome Project. NeuroImage. 2013;80:144–68. doi: 10.1016/j.neuroimage.2013.05.039 23702415
50. Nielsen JA, Zielinski BA, Ferguson MA, Lainhart JE, Anderson JS. An evaluation of the left-brain vs. right-brain hypothesis with resting state functional connectivity magnetic resonance imaging. PLoS One. 2013;8(8):e71275. doi: 10.1371/journal.pone.0071275 23967180
51. Yeh FC, Tang AN, Hobbs JP, Hottowy P, Dabrowski W, Sher A, et al. Maximum Entropy Approaches to Living Neural Networks. Entropy. 2010;12(1):89–106.
52. Becker OM, Karplus M. The topology of multidimensional potential energy surfaces: Theory and application to peptide structure and kinetics. The Journal of Chemical Physics. 1997;106(4):1495–517.
53. Csárdi G, Nepusz T. The igraph software package for complex network research. Inter Journal Complex Systems. 2006:1695.
54. Sokal R, Michener C. A statistical method for evaluating systematic relationships. University of Kansas Science Bulletin. 1958;38:1409–38.
55. Cabral J, Kringelbach ML, Deco G. Functional connectivity dynamically evolves on multiple time-scales over a static structural connectome: Models and mechanisms. NeuroImage. 2017;160:84–96. doi: 10.1016/j.neuroimage.2017.03.045 28343985
56. Stagno JR, Liu Y, Bhandari YR, Conrad CE, Panja S, Swain M, et al. Structures of riboswitch RNA reaction states by mix-and-inject XFEL serial crystallography. Nature. 2016.
57. Yuan Y, Tam MF, Simplaceanu V, Ho C. New look at hemoglobin allostery. Chem Rev. 2015;115(4):1702–24. doi: 10.1021/cr500495x 25607981
58. Vesper MD, de Groot BL. Collective dynamics underlying allosteric transitions in hemoglobin. PLoS Comput Biol. 2013;9(9):e1003232. doi: 10.1371/journal.pcbi.1003232 24068910
59. Mihailescu MR, Russu IM. A signature of the T—> R transition in human hemoglobin. Proceedings of the National Academy of Sciences of the United States of America. 2001;98(7):3773–7. doi: 10.1073/pnas.071493598 11259676
60. François-Martin C, Rothman JE, Pincet F. Low energy cost for optimal speed and control of membrane fusion. Proceedings of the National Academy of Sciences of the United States of America. 2017;114(6):1238–41. doi: 10.1073/pnas.1621309114 28115718
61. Ryham RJ, Klotz TS, Yao L, Cohen FS. Calculating Transition Energy Barriers and Characterizing Activation States for Steps of Fusion. Biophysical journal. 2016;110(5):1110–24. doi: 10.1016/j.bpj.2016.01.013 26958888
62. Smirnova YG, Marrink S-J, Lipowsky R, Knecht V. Solvent-Exposed Tails as Prestalk Transition States for Membrane Fusion at Low Hydration. J Am Chem Soc. 2010;132(19):6710–8. doi: 10.1021/ja910050x 20411937
63. Shimabukuro K, Muneyuki E, Yoshida M. An alternative reaction pathway of F1-ATPase suggested by rotation without 80 degrees/40 degrees substeps of a sluggish mutant at low ATP. Biophys J. 2006;90(3):1028–32. doi: 10.1529/biophysj.105.067298 16258036
64. Uemura S, Higuchi H, Olivares AO, De La Cruz EM, Ishiwata S. Mechanochemical coupling of two substeps in a single myosin V motor. Nat Struct Mol Biol. 2004;11(9):877–83. doi: 10.1038/nsmb806 15286720
66. Edelman GM, Gally JA. Degeneracy and complexity in biological systems. Proceedings of the National Academy of Sciences of the United States of America. 2001;98(24):13763–8. doi: 10.1073/pnas.231499798 11698650
67. Spiegler A, Hansen EC, Bernard C, McIntosh AR, Jirsa VK. Selective Activation of Resting-State Networks following Focal Stimulation in a Connectome-Based Network Model of the Human Brain. eNeuro. 2016;3(5).
68. Pillai AS, Jirsa VK. Symmetry Breaking in Space-Time Hierarchies Shapes Brain Dynamics and Behavior. Neuron. 2017;94(5):1010–26. doi: 10.1016/j.neuron.2017.05.013 28595045
70. Lopez-Persem A, Verhagen L, Amiez C, Petrides M, Sallet J. The Human Ventromedial Prefrontal Cortex: Sulcal Morphology and Its Influence on Functional Organization. J Neurosci. 2019;39(19):3627–39. doi: 10.1523/JNEUROSCI.2060-18.2019 30833514
71. Margulies DS, Ghosh SS, Goulas A, Falkiewicz M, Huntenburg JM, Langs G, et al. Situating the default-mode network along a principal gradient of macroscale cortical organization. Proc Natl Acad Sci U S A. 2016;113(44):12574–9. doi: 10.1073/pnas.1608282113 27791099
72. van den Heuvel MP, Sporns O. Rich-club organization of the human connectome. J Neurosci. 2011;31(44):15775–86. doi: 10.1523/JNEUROSCI.3539-11.2011 22049421
73. Park HJ, Lee JD, Chun JW, Seok JH, Yun M, Oh MK, et al. Cortical surface-based analysis of 18F-FDG PET: measured metabolic abnormalities in schizophrenia are affected by cortical structural abnormalities. Neuroimage. 2006;31(4):1434–44. doi: 10.1016/j.neuroimage.2006.02.001 16540349
75. Watanabe T, Rees G. Brain network dynamics in high-functioning individuals with autism. Nature Communications. 2017;8:16048. doi: 10.1038/ncomms16048 28677689
76. Michel CM, Koenig T. EEG microstates as a tool for studying the temporal dynamics of whole-brain neuronal networks: A review. Neuroimage. 2018;180(Pt B):577–93. doi: 10.1016/j.neuroimage.2017.11.062 29196270
77. Alexander DM, Trengove C, van Leeuwen C. Donders is dead: cortical traveling waves and the limits of mental chronometry in cognitive neuroscience. Cogn Process. 2015;16(4):365–75. doi: 10.1007/s10339-015-0662-4 26139038
78. Zhang H, Watrous AJ, Patel A, Jacobs J. Theta and Alpha Oscillations Are Traveling Waves in the Human Neocortex. Neuron. 2018;98(6):1269–81 e4. doi: 10.1016/j.neuron.2018.05.019 29887341
79. Roberts JA, Gollo LL, Abeysuriya RG, Roberts G, Mitchell PB, Woolrich MW, et al. Metastable brain waves. Nat Commun. 2019;10(1):1056. doi: 10.1038/s41467-019-08999-0 30837462
80. Lecrux C, Hamel E. Neuronal networks and mediators of cortical neurovascular coupling responses in normal and altered brain states. Philos Trans R Soc Lond B Biol Sci. 2016;371(1705).
81. Deco G, Ponce-Alvarez A, Mantini D, Romani GL, Hagmann P, Corbetta M. Resting-state functional connectivity emerges from structurally and dynamically shaped slow linear fluctuations. J Neurosci. 2013;33(27):11239–52. doi: 10.1523/JNEUROSCI.1091-13.2013 23825427
82. Stephan KE, Kasper L, Harrison LM, Daunizeau J, den Ouden HE, Breakspear M, et al. Nonlinear dynamic causal models for fMRI. Neuroimage. 2008;42(2):649–62. doi: 10.1016/j.neuroimage.2008.04.262 18565765
83. Vidaurre D, Quinn AJ, Baker AP, Dupret D, Tejero-Cantero A, Woolrich MW. Spectrally resolved fast transient brain states in electrophysiological data. Neuroimage. 2016;126:81–95. doi: 10.1016/j.neuroimage.2015.11.047 26631815
84. Vidaurre D, Smith SM, Woolrich MW. Brain network dynamics are hierarchically organized in time. Proc Natl Acad Sci U S A. 2017;114(48):12827–32. doi: 10.1073/pnas.1705120114 29087305
85. Ezaki T, Watanabe T, Ohzeki M, Masuda N. Energy landscape analysis of neuroimaging data. Philos Trans A Math Phys Eng Sci. 2017;375(2096).
86. Loh M, Rolls ET, Deco G. A dynamical systems hypothesis of schizophrenia. PLoS Comput Biol. 2007;3(11):e228. doi: 10.1371/journal.pcbi.0030228 17997599
87. Cabral J, Fernandes HM, Van Hartevelt TJ, James AC, Kringelbach ML, Deco G. Structural connectivity in schizophrenia and its impact on the dynamics of spontaneous functional networks. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2013;23(4):046111.