The dynamics of resting fluctuations in the brain: metastability and its dynamical cortical core

Abstract

In the human brain, spontaneous activity during resting state consists of rapid transitions between functional network states over time but the underlying mechanisms are not understood. We use connectome based computational brain network modeling to reveal fundamental principles of how the human brain generates large-scale activity observable by noninvasive neuroimaging. We used structural and functional neuroimaging data to construct whole- brain models. With this novel approach, we reveal that the human brain during resting state operates at maximum metastability, i.e. in a state of maximum network switching. In addition, we investigate cortical heterogeneity across areas. Optimization of the spectral characteristics of each local brain region revealed the dynamical cortical core of the human brain, which is driving the activity of the rest of the whole brain. Brain network modelling goes beyond correlational neuroimaging analysis and reveals non-trivial network mechanisms underlying non-invasive observations. Our novel findings significantly pertain to the important role of computational connectomics in understanding principles of brain function.

Introduction

“When we take a general view of the wonderful stream of our consciousness, what strikes us first is the different pace of its parts. Like a bird’s life, it seems to be made of an alternation of flights and perchings.” William James ¹

Survival remains the perhaps most important problem faced by brains and a key challenge is how to segregate and integrate relevant information over different timescales when faced with hostile, often constantly changing environments². Reconciling different speeds of information processing, from fast to slow, is especially important, and could be key to the relative evolutionary success of mammals whose sophisticated brains are able to combine prior information from past memories with current stimuli to predict the future and to adapt behaviour accordingly^3,4,5.

This was recognized well over a century ago by William James, generally acknowledged as one of the fathers of modern cognitive psychology¹. Speaking of this problem using the apt metaphor of the stream of consciousness, James noted that there is a different pace to its parts, comparing it to the life of a bird whose journey consists of an “alternation of flights and perchings”. In the language of today’s dynamical systems, the flights are akin to fast, segregative tendencies and the perchings to slower, integrative tendencies of the dynamic brain in action^{2, 6, 7}. In addition, motivated by recent experimental and modelling work of other labs^{8, 9}, we investigate cortical heterogeneity across areas. By optimizing the spectral characteristics of each local brain node (in the coupled network), this allowed us to discover a dynamical core of the brain, i.e. the set of brain regions, which through their oscillations are driving the rest of the brain. Furthermore, with regards to balancing the different speeds of processing, a large body of psychological research has focused on what is known as dual process theories^{10, 11}, identifying competing fast and slow systems which have to co-exist and function on multiple time-scales in order for the brain to efficiently allocate the resources necessary for survival^{12, 13}.

Yet, the temporal dynamics and underlying neural mechanisms of this temporal processing on multiple timescales are poorly understood. Here we aim to provide a better understanding of the dynamics using computational brain network modelling which has emerged as a powerful tool for investigating the causal dynamics of the human brain, when carefully constrained by functional (FC) and structural connectivity (SC) obtained from empirical neuroimaging data^{14,15,16,17,18}. This theoretical framework has been largely successful in explaining the highly structured dynamics arising from spontaneous brain activity in the so-called resting-state-networks (RSN)^19,20,21, even if the resting brain never truly rests²⁰. Efficient task-related brain activity has been shown to rely on metastability of spontaneous brain activity allowing for optimal exploration of the dynamical repertoire²² but it is not known if this metastability is maximally metastable⁶. In dynamic systems, metastability refers to a state that falls outside the natural equilibrium state of the system but persists for an extended period of time. One example of a metastable dynamical system is a winnerless competition²³, however metastability can arise from a number of underlying mechanisms and it is in this broader sense that we use the term metastability. Here, we use metastability to denote the variability of the global synchronization as measured by the Kuramoto order parameter, following the work of Wildie and Shanahan²⁴. This is motivated by previous research which has shown that this variability can be linked to the underlying metastable cluster synchronization²⁵.

We investigated the dynamics of the brain network system through a local node neural mass description based on the most general form of expressing both noisy asynchronous dynamics and oscillations, namely a normal form of a Hopf bifurcation^26,27,28. Previous research has shown the usefulness, richness and generality of this type of model for describing EEG dynamics at the local node level^{27, 28}. This normal form allowed us to fit the model to neuroimaging data over time, i.e. not only by fitting the grand average FC but also by fitting the temporal structure of the fluctuations, functional connectivity dynamics²⁹ (FCD, Fig. 1A,B).

Figure 1

Methods for measuring fit between simulated and empirical data. (A) The fitting of the FC is measured by the Pearson correlation coefficient between corresponding elements of the upper triangular part of the matrices. (B) For comparing the FCD statistics, we collected the upper triangular elements of the matrices (over all participants or sessions) and compared the simulated and empirical distribution by means of the Kolmogorov-Smirnov distance between them. The Kolmogorov–Smirnov distance quantifies the maximal difference between the cumulative distribution functions of the two samples. (C) We measure the metastability as the standard deviation of the Kuramoto order parameter across time. The Kuramoto order parameter measures the global level of synchronization of the n oscillating signals. Under complete independence, the n phases are uniformly distributed and thus R is nearly zero, whereas R = 1 if all phases are equal (full synchronization). For calculating the metastability of the empirical and simulated BOLD signals, we first band-pass filtered within the narrowband 0.04–0.07 Hz and computed the instantaneous phase φ_k(t) of each narrowband signal k using the Hilbert transform. The Hilbert transform yields the associated analytical signals. The analytic signal represents a narrowband signal, s(t), in the time domain as a rotating vector with an instantaneous phase, φ(t), and an instantaneous amplitude, A(t). Bottom panel visualizes a single example scenario (of many possible others) where the model system’s metastability increases as a function of G. We also indicate the metastability measured in empirical data. Part of figure B is based on the work of Allen and colleagues⁴⁴.

We further explored if the optimal working point where FC and FCD are fitted corresponds to a dynamical region where the global metastability of the whole brain is maximized⁶. In addition, motivated by recent experimental and modelling work of other labs^{8, 9}, we investigate cortical heterogeneity across areas. By optimizing the spectral characteristics of each local brain node (in the coupled network), this allowed us discovering a dynamical core of the brain, i.e. the set of brain regions, which through their oscillations is driving the rest of the brain. As such this investigation was designed to provide an empirical, scientific footing for James’ metaphorical speculations of the flights and perchings of human brain dynamics, and to demonstrate the potential of sophisticated brain network computational modelling to provide new insights into the causal mechanisms of neuroimaging results.

Results

The results arose from using brain network computational models for the analysis of empirical neuroimaging data characterising the functional and structural connectivity of 24 healthy human participants acquired using standard MRI techniques¹⁷ (see Methods). In particular, we were able to gain new insights on the emergence of transiently spatiotemporal structured networks among segregated brain regions by examining a whole-brain network model using a very general neural mass model known as the normal form of a Hopf bifurcation (also known as Landau-Stuart Oscillators), which is the canonical model for studying the transition from noisy to oscillatory dynamics²⁶ (Fig. 2). Here, we extended previous research on local node dynamics^{27, 28} by studying the whole-brain network dynamics, i.e. by investigating how those local noisy oscillators interact, and how the emerging whole-brain network activity relates to fMRI resting state dynamics. Within this model, each node of the network is modeled by a normal Hopf bifurcation, with an intrinsic frequency ω _i in the 0.04–0.07 Hz band (i = 1, …,n). The intrinsic frequencies were estimated directly from the data, as given by the averaged peak frequency of the narrowband BOLD signals of each brain region (see Methods). The state of each node i is determined by its phase, φ _i (t), and the interaction between nodes depends both on the structural couplings and the phase difference between the nodes. The model has only two types of control parameters, namely: one single global parameter, G, that represents the global scaling of the anatomical connectivity matrix, and the bifurcation parameters a _j for each node (see Fig. 3 and methods for the general structure and strategy of the brain network model).

Figure 2

Construction of individual brain network models. (A) The brain network model was based on individual structural connectivity (SC) matrices from 24 participants derived from tractography of DTI (left) between the 68 regions of the Desikan-Kahilly parcellation (middle). The control parameters of the models were tuned using the grand average FC and FCD derived from fMRI BOLD data (right). (B) For modelling local neural masses we used the normal form of a Hopf bifurcation, where depending on the bifurcation parameter, the local model generates a noisy signal (left), a mixed noisy and oscillatory signal (middle) or an oscillatory signal (right). It is at the border between noisy and oscillatory behaviour (middle), where the simulated signal looks like the empirical data, i.e. like noise with an oscillatory component around 0.05 Hz.

Figure 3

Fitting of the empirical data by the brain network Hopf model for different working points. (A) Level of fitting of the FC, FCD and metastability as a function of the global scaling parameter G for three different bifurcation parameters a = [−0.2 0 0.2], namely at the noisy oscillatory region, at the edge of the bifurcation and at the oscillatory regime. (B) The three measures for assessing fitting between simulated and empirical data are shown color-coded as a function of bifurcation parameter a and global scaling parameter, G. The best fitting of the three measures is obtained for a region at the brink of the Hopf bifurcation, i.e. for bifurcation parameter a, at the edge of zero on the negative side. In this region not only the correlation between the empirical and simulated FC is maximized (upper panel), but also the statistics of the rapid switching between FC(t) across time (FCD) is minimized in Kolmogorov-Smirnov sense (middle panel), and the level of metastability of the data is perfectly reproduced (bottom panel).

Maximal metastability at the optimal working point of model

Using the Hopf model, we were able discern the dynamical properties of the optimal working point of the system that is able to fit the characteristics of the empirical fMRI data. We were able to distinguish the origin of resting activity between the two hypothesized scenarios, namely: 1) noisy excursions at the edge of a critical bifurcation^{19, 20, 30, 31} or 2) metastable oscillations¹⁶. The first scenario refers to the entrainment of noisy dynamics through the underlying anatomical connectivity matrix, i.e. inducing correlations of the local noise because of the underlying SC connections. The second scenario refers to the structuring of metastable cluster synchronizations of the underlying local oscillatory dynamics through the underlying anatomical SC connections. We define metastability as the standard deviation of synchrony at the network level described by order parameter R(t), where R(t) measures the phase uniformity and varies between 0 for a fully desynchronized network and 1 for a fully synchronized network (see methods and Fig. 1C)²⁴. The present model is able to describe both types of dynamics, and the smooth transitions from one to the other, i.e. the transition from noisy to oscillatory dynamics (Fig. 2). In order to distinguish the dynamical scenario, we investigated the capabilities of the model for fitting the grand average FC and also the time dependent characteristics of the RSN as reflected in the FCD in the different dynamical working regions (i.e. as a function of the control parameters). The grand average FC describes the mean spatial structure of the resting activity, whereas the FCD captures the statistical characteristic of the temporal structure of those spatial correlations (see Methods and ref. 29).

Figure 3 shows that the best fit to the empirical data of Hopf model is found at the brink of the Hopf bifurcation. We equalized all local bifurcation parameters to a common value i.e. a _j  = a, in order to reduce the investigations to just two parameters, namely global bifurcation parameter (a) and global coupling strength (G). Figure 3 shows how the empirical data are fitted in the Hopf model for different working points. The right column of Fig. 3 shows the level of fitting of the FC, FCD and metastability. As can be seen, the best fitting of the three measures is obtained at the region on the brink of the Hopf bifurcation, i.e. for bifurcation parameter a, at the edge of zero on the negative side, such that the oscillators remain damped still. In this region not only the correlation between the empirical and simulated FC is maximized, but also the statistics of the rapid switching between FC(t) across time (FCD) is minimized in Kolmogorov-Smirnov sense, and the level of metastability of the data is reproduced. The fitting of the FC was measured by the Pearson correlation coefficient between corresponding elements of the upper triangular part of the matrices (see Fig. 1 and Methods). For comparing the FCD statistics, we collected the upper triangular elements of the matrices (over all participants or sessions) and compared the simulated and empirical distribution by means of the Kolmogorov-Smirnov distance between them (see Methods).

Furthermore, the results showed that only in the region at the border between noisy and oscillatory behaviour, is where the signals resembles the data, i.e. like noise with an oscillatory component around 0.05 Hz (Fig. 2). The first three columns of Fig. 3 show the dependence of those measurements as a function of the global scaling parameter G for three specific values of the bifurcation parameter a, namely at the noisy region, at the edge of the bifurcation and at the oscillatory regime. Clearly, the best results are obtained for the second column (at the edge of the bifurcation). The same panel shows that the FCD is the best constraining measure. There is a broad range of G where the FC and the metastability is well fitted, but only a relative narrow range where the FCD statistics is minimal, i.e. maximally fitted. In other words, the spatiotemporal structure of the FC is more informative than the grand average of the FC (i.e. the “classical” RSN). This is important, because until now, brain network models have always been fitted with the grand average FC - but see also the work of Hansen and colleagues²⁹.

We would like to remark that Fig. 3 characterizes some of the bifurcation behaviour of the whole system. Indeed, the metastability for example serves as a network metric and characterizes the variability of this global synchronization as a function of those two control parameters. All three parameter spaces in Fig. 3B, in conjunction, present a full picture of the spatiotemporal organization of the system. The three metrics characterize computationally the bifurcation properties of the full network dynamics.

Perhaps most importantly, as shown in Fig. 3, the brain network model shows maximal metastability at the optimal working point of the model (a = 0 and G = 2.85), where the metastability is reflecting the variability of the synchronization between different nodes, i.e. the fluctuations of the states of phase configurations as a function of time²⁴. Further characterisation of these results is shown in Fig. 4 which shows the optimal working point at the edge of the Hopf bifurcation (i.e. bifurcation parameter a = 0), the FC, FCD and FCD statistics for three levels of global coupling G namely low, optimal and large. For comparison, the same matrices and distributions are plotted on the rightmost column for the empirical data (Fig. 4B). Only the FCD and its statistics (bottom row) are constraining enough for optimizing the working point. Please note that for low G the FCD statistics does not show any switching between states in the RSN and that for very large G there are too much switching between states.

Figure 4

Fitting to the grand average FC is a necessary but not sufficient condition for best empirical fitting. (A) The figure shows the result of fitting the model to the empirical as a function of the global coupling parameter, G, at the optimal working point at the edge of the Hopf bifurcation (i.e. bifurcation parameter a = 0). Three different coupling points were selected (low, optimal and large in the three columns) and we show the resulting FC correlation, FCD correlations and FCD histogram. Note that for low G the FCD statistics does not show any switching between RSN and that for very large G there are too much switching between states. (B) For comparison, the same matrices and distributions are plotted for the empirical data. Note how only the FCD (row 2) and its statistics (row 3) are constraining enough for optimizing the working point the model to fit the empirical data (compare the distributions in row 3 and compare plots for FCD and FC fitting in row 1).

Dynamical core: contribution of individual brain regions to dynamics

In order to obtain information about the dynamical characteristics of each single brain area and to generate a heterogeneous brain network model (i.e. with different dynamics at each node), we optimized each single bifurcation parameter a _j independently by fitting for each value of global coupling G the spectral characteristics of the simulated and empirical BOLD signals at each brain area (see Methods). The main results are plotted in Fig. 5, where Fig. 5A shows the evolution of the fitting of the FC and FCD statistics as a function of G. For large enough value of the global coupling a good fitting of both is obtained, i.e. large correlation between the empirical and simulated grand average FC and low difference in the statistics of the empirical and simulated FCD (Kolmogorov-Smirnov distance). Please note that Fig. 5A is generated in a different way than Fig. 4A (which uses only the optimum fit a = 0 for all regions and G = 2.85). Instead, in Fig. 5A, for each G we optimize the bifurcation value, a, for each region (shown in Fig. 5B). As can be seen at a critical value of G, the bifurcation values remain the same, only scaled. Thus the FCD fit in Fig. 5A will asymptote as G increases.

Figure 5

Spectral characteristics of the dynamical core of the human brain. To generate a heterogeneous brain network model (i.e. with different dynamics at each node), we optimized each single bifurcation parameter independently by fitting for each value of global coupling G the spectral characteristics of the simulated and empirical BOLD signals at each brain area. (A) The evolution of the fitting of the FC and FCD statistics as a function of G. For large enough value of the global coupling a good fitting of both is obtained, i.e. large correlation between the empirical and simulated grand average FC and low difference in the statistics of the empirical and simulated FCD (Kolmogorov-Smirnov distance). (B) The evolution of the single values of the local bifurcations parameters as a function of the global coupling G. For low values of G homogeneous local bifurcation parameters around zero are obtained. When the level of fitting improves for larger values of G a more heterogeneous distribution of is obtained. (C) The local bifurcation parameters for each region for the uncoupled network (i.e. G = 0) and for the optimal coupling (G = 5.4). If the network is uncoupled, each single brain area fitted the spectral characteristics of the empirical BOLD signals in a very homogeneous way by local bifurcations parameters at the edge of the local Hopf bifurcation, i.e. at zero. (D) When the whole-brain network is coupled, we can discover the “true” intrinsic local dynamics that fits the local empirical BOLD characteristics and the global quantities FC, FCD and metastability.

For optimizing a _j values, we use a greedy optimisation strategy, where we iteratively increase or decrease the a _j value according to the local power of the signal in a given region j. Greedy algorithms exploit local optima, but often approximate optimal solutions well in reasonable time and produce good results as shown in Fig. 5. The local bifurcation parameters for each region for the uncoupled network (i.e. G = 0) and for the optimal coupling (G = 5.4) can be seen in Fig. 5C. If the network is uncoupled, each single brain area fitted the spectral characteristics of the empirical BOLD signals in a very homogeneous way by local bifurcations parameters at the edge of the local Hopf bifurcation, i.e. at zero. When the brain network is coupled, the “true” intrinsic local dynamics for the profile of optimal local bifurcation parameters a _j observed at that point that fit the local empirical BOLD characteristics and the global quantities FC, FCD and metastability (Fig. 5D).

Brain regions, for which best predictions were achieved in an oscillatory mode, i.e. with bifurcation parameters a > 0.1 are visualised in Fig. 6. We found that the dynamical core within this parcellation consisted of eight lateralised brain regions: medial orbitofrontal cortex, posterior cingulate cortex and transverse temporal gyrus in the right hemisphere, and caudal middle frontal gyrus, precentral gyrus, precuneus cortex, rostral anterior cingulate cortex and transverse temporal gyrus in the left hemisphere. Those nodes working at the edge of the bifurcation are highlighted as a “dynamical core” whose perturbations can propagate in an optimal way to the rest of the network.

Figure 6

Dynamical core in the human brain. The figure shows the dynamical core regions on the edge of bifurcation (location of neural masses shown in light blue and transparent blue for the full region). These are the nodes with the ability to react immediately to changes in the predicted input and thus likely to drive the rest of the brain networks. The eight regions are clearly lateralised; and in the right hemisphere encompass medial orbitofrontal cortex, posterior cingulate cortex and transverse temporal gyrus, while in the left hemisphere include caudal middle frontal gyrus, precentral gyrus, precuneus cortex, rostral anterior cingulate cortex and transverse temporal gyrus. Interestingly, some of these regions are part of the default mode network (medial orbitofrontal cortex, posterior cingulate cortex and precuneus cortex) while others have been implicated in memory processing (parahippocampal and transverse temporal gyrus), auditory processing (transverse temporal gyrus), selection for action (rostral anterior cingulate cortex and caudal middle frontal gyrus) and motor execution (precentral gyrus).

Discussion

We provide mechanistic explanations of the complex spatiotemporal dynamics of brain function arising from James’ early speculations¹ to much more detailed scientific enquiry^{2, 32,33,34}. This confirms that brain results from complex interactions in a system of non-linearly coupled, non-linear oscillatory processes which display dynamical system phenomena such as multiple stable states, instability, state transitions and metastability, of which the latter has been proposed to form a core dynamical description of coordinated brain and behavioral activity⁶.

In the 1980s the physicist Hermann Haken suggested to mechanistically interpret brain processes of segregation and integration as a sequence of semistable states, so-called saddle states³⁵. He proposed to view the complex integrative and segregative tendencies as expressions of emergent lower-dimensional behavior of collective variables, which he termed ‘order parameters’. Scott Kelso popularized this concept using the term ‘metastability’ based on his brain-behaviour experiments and drawing inspiration from other researchers including Rodolfo Llinás and Francisco Varela^{33, 34}. He generalized metastability to include the oscillatory states of brain processes found in between complete synchronization and independence^{32, 36}. Later research has formalized these concepts more rigorously, e.g. via the heteroclinic channel^{37, 38} and Structured Flows on Manifolds (SFM)^{39, 40}.

We shed new causal light on the mechanisms underlying RSNs by extending previous research which has demonstrated the existence of RSNs, i.e. brain networks correlated within the grand average FC during resting state^{21, 41, 42}. FC has become routinely used as a biomarker in various clinical applications, even though its predictive value holds only for group analyses, and not currently for the individual⁴³. This problem arises most likely from the lack of taking time into account, i.e. the non-stationary nature of the resting state dynamics^{44, 45}. Hansen and colleagues demonstrated that the grand average FC is more closely linked to the SC and linear models of FC²⁹. When non-linearities are considered in the network models, the spatiotemporally dynamic repertoire of the network is significantly enhanced and the resting state dynamics shows the non-stationary FCD, which expresses itself as the switching dynamics of the FC. While Hansen and colleagues proposed FCD as a novel biomarker and demonstrated that all known RSNs can be derived from the non-linear network dynamics of FCD, they did not fit the model to the empirical functional time series data. The patterns in the FCD matrix arise from what is essentially a random process and thus different for different measurements. This renders the fitting process for brain network models more complex than fitting with the grand average FC, for which a Pearson correlation across empirical and simulated FC matrices is sufficient.

We have addressed this issue through a systematic fitting approach of the random process in FCD to the empirical data. The conjunction of using sophisticated fitting and systematic parameter analysis allowed us to test the mechanistic hypotheses underlying the resting state, i.e. whether the brain at rest operates close to the edge of a bifurcation and/or occupies a metastable state. Both scenarios can be mechanistically realized by non-linearly coupling Hopf bifurcators²⁶. Hopf oscillators have been used previously in connectome-based modelling of resting state dynamics in EEG/MEG and fMRI¹⁴, as well as for the modelling of the detailed temporal dynamics in EEG/MEG^{27, 28}. The usage here though is different from the previous research, since the Hopf oscillators act as the sources of BOLD signal in the connectome based network model. Ghosh and colleagues used the Hopf oscillators as the sources of the electrophysiological signal and employed the Balloon Windkessel to derive the BOLD signal⁴⁶. Given this interpretation, they needed to include all the signal transmission delays. In our present approach, the oscillation frequencies are significantly slower and thus permit the neglect of the time delays, which simplifies the computational effort of the simulation and thus the computational fitting of the models against empirical data.

Our key finding is the demonstration that the optimal operating regime is at the edge of the local Hopf bifurcation, i.e. a balance of noisy excursions in the oscillatory state. We not only were able to demonstrate that previous findings on the optimal operating point based on grand-average FC hold true if we take into account the temporal dynamics of FC, i.e. FCD. We also demonstrated that a better way of constraining brain network models is by not only fitting the grand average FC, but by also fitting the temporal structure of the fluctuations using the FCD.

Another remarkable and important finding is that high metastability is only present in a narrow range of bifurcation parameter when a is close to the edge of the bifurcation. In other words, the FCD of the spontaneous resting state, in conjunction with brain network modelling provide evidence that the brain at rest is maximally metastable, refining and demonstrating the hypothesis of Tognoli and Kelso⁶. Note that there is also a region for very small G and positive a (oscillatory regime) where a relatively good fitting is obtained. This dynamic regime was previously observed with a pure oscillatory Kuramoto model of the BOLD signals at the mesoscopic level⁴⁷. Nevertheless, the level of fitting for the FC, metastability and even FCD is not as good as the one obtained in the region at the edge of the Hopf bifurcation. On the other hand, besides the extreme sensitivity of that working point (ultra-narrow regime of optimality) which means that the result is not so robust, the qualitative description of the BOLD signals is not realistic in the pure oscillatory regime in comparison with the noisy/oscillatory excursions evidenced in the regime of the bifurcation parameter a near zero.

For constructing a heterogeneous brain network model with different local parameter values, we took into account the spectral information of the BOLD data. We addressed the question if the oscillations at the individual nodes play a mechanistic role for the emergence of FC/FCD. In particular, we identified a cortical core of eight brain regions with the optimal fit of bifurcation parameter a close to the edge of bifurcation. We propose to call this the dynamical cortical core of the brain. Interestingly, three of these regions (the medial orbitofrontal cortex, posterior cingulate cortex and precuneus cortex) are part of the default mode network and thus re-experience past events and pre-experience possible future events^{48, 49}. In this vein other regions (parahippocampal and transverse temporal gyrus) have also been implicated in memory processing and may thus perhaps be helping integrate information over different timescales, binding fast and slow processes over time². This information is always contextual and in the noisy, unpredictable scanner it is perhaps not surprising that the brain is attending to the auditory signals (transverse temporal gyrus). As such this information processing is available for conflict monitoring and selection for action (rostral anterior cingulate cortex and caudal middle frontal gyrus) and motor execution (precentral gyrus)⁵⁰. Equally, the involvement of the cingulate cortex is interesting given that this region recently has been shown to be part of the common neurobiological substrate for mental illness across across six diverse diagnostic groups (schizophrenia, bipolar disorder, depression, addiction, obsessive-compulsive disorder, and anxiety) based on a meta-analysis of grey matter loss in 193 neuroimaging studies of 15892 individuals⁵¹. This reinforces the potential use of brain network computational modelling for understanding the underlying mechanisms of neuropsychiatric disorders⁵². The right-handed quality of Fig. 6, presumably arises from the specifics of the data used to fit the model and we will be exploring its biological validity in subsequent studies with larger group sizes.

Although the bifurcation parameter does not have a direct biophysical correlate, it seems to be involved in mediating biophysical effects. Another note of caution: we have presented a mesoscopic phenomenological model, i.e. the dynamical equation corresponds directly to the measured BOLD signals and not to the neural signals. It is a phenomenological model since the real coupling between regions does not, of course, occur between the hemodynamic signals, but between the underlying neural activity. Still, by using this simplified approach, we do not need to convolve modelled neural signals with the haemodynamic response function⁴⁷.

We therefore propose that in future both the global bifurcation parameter as well as the individual parameters could potentially serve as biomarkers for disease. It will be important to explore the changes for different brain diseases, e.g. within a standardized framework for connectome-based modelling such as The Virtual Brain (TVB)^{53, 54}, and applications such as fitting of TVB’s dynamic regime and TVB Processing pipeline¹⁷.

Overall, we have shown that neuroimaging data can be causally analysed by constructing a relatively simple brain network computational model using a Hopf bifurcation. This model was shown to be maximally metastable at the optimal fitting with the spatiotemporal dynamics of spontaneous brain activity. This dynamical regime may well allow for the optimal integration and segregation of fast and slow information over different time-scales, the “flights” and “perchings” of the stream of consciousness alluded to by William James over 100 years ago.

Methods

Ethics Statement

All participants of this study gave written informed consent before the study, which was performed in compliance with the relevant laws and institutional guidelines and approved by the ethics committee of the Charité University Berlin.

Empirical MRI Data Collection

Structural data from DTI and resting-state BOLD signal time series were acquired for 24 healthy participants (age between 18 and 33 years old, mean 25.7, 12 females, 12 males). A full description of the generation of SC and FC matrices from those data can be found in ref. 17. Here, we provide a quick overview of the employed methods. Empirical data were acquired at Berlin Center for Advanced Imaging, Charité University Medicine, Berlin, Germany. For simultaneous EEG-fMRI^{55, 56}, participants were asked to stay awake and keep their eyes closed. No other controlled task had to be performed. In addition, a localizer, DTI and T2 sequence were recorded for each participant. MRI was performed using a 3 Tesla Siemens Trim Trio MR scanner and a 12-channel Siemens head coil. Specifications for the employed sequences can be found in ref. 56. For each participant anatomical T1-weighted scans were acquired. DTI and GRE field mapping were measured directly after the anatomical scans. Next, functional MRI (BOLD-sensitive, T2*-weighted, TR 1940 ms, TE 30 ms, FA 78°, 32 transversal slices (3 mm), voxel size 3 × 3 × 3 mm, FoV 192 mm, 64 matrix) was recorded simultaneously to the EEG recording.

MRI Data Analysis

Processing steps executed by the public Berlin automatized processing pipeline⁵⁶ comprised 1) preprocessing of T1-weighted scans, cortical reconstruction, tessellation and parcellation, 2) transformation of anatomical masks to diffusion space, 3) processing of diffusion data, 4) transformation of anatomical masks to fMRI space, 5) Processing of fMRI data.

Anatomical MRI Data Analysis

The highly resolved anatomical images are important to create a precise parcellation of the brain. For each of those parcellated units, empirical functional data time series are spatially aggregated. T1-weighted images are pre-processed using FREESURFER including probabilistic atlas based cortical parcellation, here using Desikan-Killany (DK) atlas⁵⁷ (Table 1). This generates volumes that contain all cortical and subcortical parcellated regions with corresponding region labels used for fiber-tracking and BOLD time-series extraction.

Table 1 Anatomical labels for the 68 regions in the Desikan-Kahilly parcellation.

Empirical DTI Data Analysis and Tractography

Tractography requires binary WM masks to restrict tracking to WM voxels. Upon extraction of gradient vectors and values (known as b-table) using MRTrix, dw-MRI data are pre-processed using FREESURFER. Besides motion correction and eddy current correction (ECC) the b0 image is linearly registered (6 degrees of freedom, DOF) to the participant’s anatomical T1-weighted image and the resulting registration rule is stored for later use. We transformed the high-resolution mask volumes from the anatomical space to the participant’s diffusion space, to further use it for fiber tracking. The cortical and subcortical parcellations are resampled into diffusion space, one time using the original 1 mm isotropic voxel size (for subvoxel seeding) and one time matching that of our dw-MRI data, i.e., 2.3 mm isotropic voxel size. During MRTrix pre-processing diffusion tensor images that store the diffusion tensor (i.e., the diffusion ellipsoid) for each voxel location are computed. Based on that, a fractional anisotropy (FA) and an eigenvector map are computed and masked by the binary WM mask created previously. For subsequent fiber-response function estimation, a mask containing high-anisotropy voxels is computed. Fibre orientation distributions are estimated using constrained spherical deconvolution⁵⁸ based on a response function estimated in voxels that are expected to contain a single, coherently-oriented fibre bundle (commands dwi2response tournier and dwi2fod; see MRTrix Documentation: http://mrtrix.readthedocs.io/en/latest/). In order to resolve crossing pathways, fibers are prolonged by employing a probabilistic tracking approach as provided by MRTrix. In order to exclude spurious tracks, three types of masks are used to constrain tracking: seeding-, target- and stop-masks. In order to restrict track-prolongation to WM, a WM-mask that contains the union of GM-WM-interface and cortical WM voxels is defined as a global stop mask for tracking. To address several confounds in the estimation of connection strengths (information transmission capacities), a new seeding and fiber aggregation strategy was employed developed for this pipeline and described in detail in ref. 17. In combination with a new aggregation scheme, it is based on an appropriate selection of seed voxels and controlling for the number of generated tracks in each seed voxel. Instead of using every WM voxel, tracks are initiated from GM-WM-interface voxels and a fixed number of tracks are generated for each seed-voxel. Since a GM parcellation-based aggregation is performed, each seed-mask is associated with a ROI of the GM atlas. Along with seeding-masks complementary target-masks are defined specifying valid terminal regions for each track that was initiated in a specific seed voxel. The capacity measures that we derive between each pair of regions are intended to estimate the strength of the influence that one region exerts over another, i.e., their SC. In order to improve existing methods for capacities estimation the approach makes use of several assumptions with regard to seed-ROI selection, tracking and aggregation of generated tracks¹⁷. Upon tractography the pipeline aggregates generated tracks to structural connectome matrices. The weighted distinct connection count used here divides each distinct connection by the number of distinct connections leaving the seed-voxel (yielding asymmetric capacities matrix). Values have been normalized by the total surface area of the GWI of a participant.

Empirical fMRI Data Analysis

In order to generate the FC matrices, FSL’s FEAT pipeline is used to perform the following operations: deleting the first five images of the series to exclude possible saturation effects in the images, high-pass temporal filtering (100 seconds high-pass filter), motion correction, brain extraction and a 6 DOF linear registration to the MNI space. Functional data is registered to the participant’s T1-weighted images and parcellated according to FREESURFER’s cortical segmentation. By inverting the mapping rule found by registration, anatomical segmentations are mapped onto the functional space. Finally, average BOLD signal time series for each region are generated by computing the mean over all voxel time-series for each region. From the region wise aggregated BOLD data, FC matrices are computed within MATLAB using and Pearson’s linear correlation coefficient as FC metrics. We did not perform global signal regression on data.

Brain Network Model

The brain network model consists of 68 coupled brain areas (nodes) derived from the parcellation explained above. The global dynamics of the brain network model used here results from the mutual interactions of local node dynamics coupled through the underlying empirical anatomical structural connectivity matrix C _ij (see Fig. 2). The structural matrix C _i denotes the density of fibres between cortical area i and j as extracted from the DTI based tractography (scaled to a maximum value of 0.2). The local dynamics of each individual node is described by the normal form of a supercritical Hopf bifurcation, which is able to describe the transition from asynchronous noisy behavior to full oscillations. Thus, in complex coordinates, each node j is described by following equation:

$$\frac{d{z}_{j}}{dt}=z[{a}_{j}+i{\omega }_{j}-|{z}_{j}^{2}|]+\beta {\eta }_{j}(t)$$

(1)

where

$${z}_{j}={\rho }_{j}{e}^{i{\theta }_{j}}={x}_{j}+i{y}_{j}$$

(2)

and η _i (t) is additive Gaussian noise with standard deviation β=0.02. This normal form has a supercritical bifurcation at a _j  = 0, so that for a _j  < 0 the local dynamics has a stable fixed point at z _j  = 0 (which because of the additive noise corresponds to a low activity asynchronous state) and for a _j  > 0 there exists a stable limit cycle oscillation with frequency ${f}_{j}={\omega }_{j}/2\pi$. We insert equation 2 in equation 1 and separate real part in equation 3 and imaginary part in equation 4.

Thus, the whole-brain dynamics is defined by following set of coupled equations:

$$\frac{d{x}_{j}}{dt}=[{a}_{j}-{x}_{j}^{2}-{y}_{j}^{2}]{x}_{j}-{\omega }_{j}{y}_{j}+G\,\sum _{i}{C}_{ij}({x}_{i}-{x}_{j})\beta {\eta }_{j}(t)$$

(3)

$$\frac{d{y}_{j}}{dt}=[{a}_{j}-{x}_{j}^{2}-{y}_{j}^{2}]{y}_{j}+{\omega }_{j}{x}_{j}+G\,\sum _{i}{C}_{ij}({y}_{i}-{y}_{j})\beta {\eta }_{j}(t)$$

(4)

Please note that following the literature from physics, the equations are written in Cartesian, rather than polar coordinates^59,60,61,62. We couple the equations using the common difference coupling, which approximates the simplest (linear) part of a general coupling function. These equations are valid in the weakly coupled oscillator limit, in which the coupling preserves the periodic orbit of the uncoupled oscillators. If the linear coupling (following a Taylor expansion of the full coupling) does not exist, the next non-vanishing higher order term should be considered, which is a case we do not address here (please see Kuramoto⁵⁹ (see Eq 5.3.1) and Pikovsky, Arkady and Kurths⁶⁰ (see Eq. 8.12) for more detailed analytic treatments of the equations).

In the latter equations, G is a global scaling factor (global conductivity parameter scaling equally all synaptic connections). The global scaling factor G and the bifurcation parameters a _j are the control parameters with which we study the optimal dynamical working region where the simulations maximally fit the empirical FC and the FCD. We model with the variables x _j the BOLD signal of each node j. The empirical BOLD signals were band-pass filtered within the narrowband 0.04–0.07 Hz. This frequency band has been mapped to the gray matter and it has been shown to be more reliable and functionally relevant than other frequency bands^63,64,65,66. Within this model, the intrinsic frequency ω _j of each node is in the 0.04–0.07 Hz band (i = 1, …, n). The intrinsic frequencies were estimated from the data, as given by the averaged peak frequency of the narrowband BOLD signals of each brain region.

Grand average FC and FCD matrices

The grand average FC is defined as the matrix of correlations of the BOLD signals between two brain areas over the whole time window of acquisition. In order to characterize the time dependent structure of the resting fluctuations, we estimate the FCD matrix²⁹ (see Fig. 1). Each full-length BOLD signal of 22 min is split up into M=61 sliding windows of 60 sec, overlapping by 40 sec. For each sliding window, centered at time t, we calculated a separate FC matrix, FC(t). The FCD is a MxM symmetric matrix whose (t1, t2) entry is defined by the Pearson correlation between the upper triangular parts of the two matrices FC(t1) and FC(t2). Epochs of stable FC(t) configurations are reflected around the FCD diagonal in blocks of elevated inter-FC(t) correlations.

The grand average FC and the FCD matrices were estimated for the recordings of each of the 24 participants as well as for 24 simulations of 22 minutes of the computational model. We compared the FC matrices of the model (averaged Fisher’s z-transformed over the 24 sessions) and the empirical data (averaged Fisher’s z-transformed over the 24 participants), adopting as a measure of similarity between the two matrices the Pearson correlation coefficient between corresponding elements of the upper triangular part of the matrices. For comparing the FCD statistics, we collected the upper triangular elements of the matrices (over all participants or sessions) and generated the distribution of them. Then, we compared the simulated and empirical distribution by means of the Kolmogorov-Smirnov distance between them. The Kolmogorov–Smirnov distance quantifies the maximal difference between the cumulative distribution functions of the two samples.

Metastability

Here, we refer to metastability as a measure of how variable are the states of phase configurations as a function of time, i.e. how the synchronization between the different nodes fluctuates across time²⁴. Thus, we measure the metastability as the standard deviation of the Kuramoto order parameter across time. The Kuramoto order parameter is defined by following equation:

$$R(t)=|\sum _{k=1}^{n}{e}^{i{\phi }_{k}(t)}|/n$$

(5)

where φ _k (t) is the instantaneous phase of each narrowband BOLD signal at node k. The Kuramoto order parameter measures the global level of synchronization of the n oscillating signals. Under complete independence, the n phases are uniformly distributed and thus R is nearly zero, whereas R=1 if all phases are equal (full synchronization). Thus, for calculating the metastability of the empirical and simulated BOLD signals, we first band-pass filtered within the narrowband 0.04–0.07 Hz (as previously explained) and computed the instantaneous phase φ _k (t) of each narrowband signal k using the Hilbert transform. The Hilbert transform yields the associated analytical signals. The analytic signal represents a narrowband signal, s(t), in the time domain as a rotating vector with an instantaneous phase, φ _k (t), and an instantaneous amplitude, A(t), i.e., $s(t)=A(t)\cos (\phi (t))$. The phase and the amplitude are given by the argument and the modulus, respectively, of the complex signal z(t), given by $z(t)=s(t)+i.H[s(t)]$, where i is the imaginary unit and H[s(t)] is the Hilbert transform of s(t).

Local Optimization of Brain Nodes

The local optimization of each single bifurcation parameter a _j is based on the fitting of the spectral information of the empirical BOLD signals in each node. In particular, we aim to fit the proportion of power in the 0.04–0.07 Hz band with respect to the 0.04–0.25 Hz band (i.e. we remove the smallest frequencies below 0.04 Hz and consider the whole spectra until the Nyquist frequency which is 0.25 Hz)⁴⁷. For this, we filtered the BOLD signals in the 0.04–0.25 Hz band, and calculated the power spectrum P _j (f) for each node j. We define the proportion,

$${p}_{j}=\frac{{\int }_{0.04}^{0.07}{P}_{j}(f){\rm{d}}f}{{\int }_{0.04}^{0.25}{P}_{j}(f){\rm{d}}f}$$

(6)

and update the local bifurcation parameters by a gradient descendent strategy, i.e.:

$${a}_{j}={a}_{j}+\eta ({p}_{j}^{empirical}-{p}_{j}^{simulated})$$

(7)

until convergence. We used here η = 0.1. The updates of the a _j values are done in each optimization step in parallel.