Abstract
The ongoing activity of the brain at rest, i.e., under no stimulation and in absence of any task, is astonishingly highly structured into spatiotemporal patterns. These spatiotemporal patterns, called resting state networks, display low-frequency characteristics (<0.1 Hz) observed typically in the BOLD-fMRI signal of human subjects. We aim here to understand the origins of resting state activity through modeling via a global spiking attractor network of the brain. This approach offers a realistic mechanistic model at the level of each single brain area based on spiking neurons and realistic AMPA, NMDA, and GABA synapses. Integrating the biologically realistic diffusion tensor imaging/diffusion spectrum imaging-based neuroanatomical connectivity into the brain model, the resultant emerging resting state functional connectivity of the brain network fits quantitatively best the experimentally observed functional connectivity in humans when the brain network operates at the edge of instability. Under these conditions, the slow fluctuating (<0.1 Hz) resting state networks emerge as structured noise fluctuations around a stable low firing activity equilibrium state in the presence of latent “ghost” multistable attractors. The multistable attractor landscape defines a functionally meaningful dynamic repertoire of the brain network that is inherently present in the neuroanatomical connectivity.
Introduction
In the absence of external stimulation, any physical system is in its equilibrium state, which is often characterized by the lowest level of activity of that system. Eventually, if noise is inherent in the system, the fluctuations drive the system out of its equilibrium state resulting in low-amplitude random activity. The brain is a particular case of a noisy physical system composed of neurons interconnected through synapses in brain areas, whose activity is mainly characterized by the level of spiking activity in those areas. Nevertheless, the equilibrium state of the brain, i.e., the spontaneous, not stimuli- or task-evoked brain activity during rest, does not reflect just a trivial random activity as one may naively expect. During the last decade, numerous experimental investigations have shown that spontaneous brain activity during rest is highly structured into characteristic spatiotemporal patterns, the so-called resting-state networks (RSNs) (Biswal et al., 1995; Greicius et al., 2003; Fox et al., 2005, 2007; Fransson, 2005; Raichle and Mintun, 2006; Rogers et al., 2007; Vincent et al., 2007).
RSNs reflect the anatomical connectivity between brain areas in a network but cannot be understood in those terms alone (Bullmore and Sporns, 2009). The missing link for understanding the formation and dissolution of RSNs is the dynamics (Deco et al., 2011). Theoretical models allowed us to study the relation between anatomical structure and RSN. In particular, Ghosh et al. (2008) and Deco et al. (2009) were able to explain the formation and dissolution of slowly fluctuating RSNs by considering a simple local oscillatory dynamics at each node. Other RSN modeling efforts included more detailed physiological models for the dynamics of brain areas, but then again imposed specific constraints upon the network dynamics, mostly for reasons of computational efficiency [such as epileptiform node dynamics of a given brain area (Honey et al., 2007) or truncation of time delays via signal transmission and rescaling of the connectivity (Izhikevich and Edelman, 2008)]. A more detailed and complete physiological model for the dynamics of individual brain areas will allow making the link between neurophysiological parameters and RSN dynamics. To this aim, we formulate a detailed global attractor model of the brain network, which offers a realistic mechanistic model at the level of each single brain area based on spiking neurons and realistic AMPA, NMDA, and GABA synapses.
With this in mind, we here demonstrate the emergence of the slowly fluctuating (<0.1 Hz) RSNs as noise-driven transient fluctuations around the stable equilibrium state corresponding to low firing activity in all neurons in all areas. The best match between simulated and experimental resting state functional connectivity as measured in humans with fMRI occurs right below the threshold separating the trivial low firing activity equilibrium state from the multistable landscape of attractors with high firing activity in certain brain areas. The latter global attractor structure represents the brain's dynamic repertoire and shows the highest signal complexity as evidenced by entropy.
Materials and Methods
Empirical neuroanatomical and functional connectivity matrix.
Neuroanatomical connections in 5 healthy right-handed male human subjects were extracted by using diffusion spectrum imaging (DSI) white matter tractography (Hagmann et al., 2008; Honey et al., 2009). This neuroanatomical matrix expresses the density with which two different brain areas are connected through white matter fiber tracts, where we used a segmented gray matter parcellation into 66 areas. The neuroanatomical matrix was finally averaged across the 5 human subjects. Figure 1a shows graphically the structure of the connectivity matrix by encoding the strengths of the different connections in a color map. The connectivity matrix is symmetric at the voxel level, due to the fact that tractography cannot distinguish the directionality of the fibers. The internal local connectivity between neurons within a given brain area is considered explicitly in the spiking attractor model adopted (described below), so that at the global level described by the neuroanatomical connectivity the self-connection of a region to itself is set to 0 in the connectivity matrix. We order the different brain areas in the neuroanatomical connectivity matrix according to modules that have substantially denser connectivity within the module than with the complementary part of the network. Furthermore, homotopic regions in the two cerebral hemispheres were arranged symmetrically with respect to the center of the matrix. This reordering reveals graphically the small-world structure of brain networks through the presentation of clusters of varying size (Watts et al., 1998). In particular, the reordering of the connectivity matrix (Fig. 1a) shows the presence of clusters of nodes that are more connected inside than outside the cluster to which they belong, confirming previous observations (Bullmore and Sporns, 2009).
a, Neuroanatomical Connectivity matrix obtained by DSI and tractography after averaging across 5 human subjects. b, Empirical functional connectivity matrix of the BOLD activity averaged across the same 5 subjects during 20 min under resting state conditions (data from Hagmann et al., 2008).
Resting state activity was also obtained for the same 5 subjects by measuring the corresponding fMRI BOLD signal during 20 min in absence of stimulation or task. After regressing out the global signal (Fox et al., 2005) and averaging across subjects an empirical functional connectivity matrix was obtained (for details, see Honey et al., 2009). This empirical functional connectivity matrix reflects the correlation of the BOLD activity between different brain areas at rest. Figure 1b plots the resulting empirical functional connectivity matrix in a color map.
Global cortical model.
The global architecture of the model is shown in Figure 2. Each brain area is modeled by a local spiking attractor network consisting of mutually interconnected populations of excitatory pyramidal neurons and GABAergic inhibitory neurons. This population model serves as the representation of a node in the large-scale network, in which the structural connectivity matrix is composed of neuroanatomical connections between distinct brain areas in the human. The structural connectivity is obtained from DSI, which forces the connectivity matrix to be symmetric. This constraint will be only limiting, if the symmetry breaking of the real connectivity is large, since then it may introduce additional oscillatory behaviors of the network (Knock et al., 2009; Jirsa et al., 2010). At the current stage, there is no solution to overcome this constraint, though one could imagine that merging directed structural information from primate data (such as the Cocomac database; Kötter, 2004) with tractographic data from the human might bear promise.
Global cortical network. The network contains at each node (brain area) excitatory pyramidal cells (red triangles) and inhibitory interneurons (blue circles). Neurons are fully connected and clustered into excitatory and inhibitory populations (large circles). Black and yellow arrows indicate excitatory and inhibitory recurrent connections between neurons in a local brain area respectively, and red arrows indicate excitatory connections between neurons in different brain areas.
Here, the interconnection between different brain areas is specified by the neuroanatomical human matrix described above. We consider here that the white matter tract connections between two distinct brain areas describe synaptic connections between pyramidal neurons in those areas. We weight those inter-areal connections by the strength specified in the neuroanatomical matrix (numbers of fibers connecting those regions) and by a global factor denoted by W, which we will vary systematically to study the dynamics and fix points of the global cortical system (attractors) as a function of it.
Local brain area model: spiking attractor network.
For modeling a local brain area (i.e., a node in the global network) we have taken a two-tiered approach, that is modeling the population of single neurons, as well as a mean-field approximation (after Brunel and Wang, 2001).
For the microscopic level of description, we use a biophysically realistic attractor network model consisting of integrate-and-fire spiking neurons with excitatory (AMPA and NMDA) and inhibitory (GABA-A) synaptic receptor types (Brunel and Wang, 2001). This type of attractor network of spiking neurons is a dynamical system that in general has the tendency to settle in stationary states, fix points called “attractors,” typically characterized by a stable pattern of firing activity. External or even intrinsic noise that appears in the form of finite size effects could provoke destabilization of an attractor inducing therefore transitions between different stable attractors. The dynamics of the network can be described by coupling the dynamical equations describing each neuron and the synaptic variables associated with their mutual coupling. We use here an Integrate-and-Fire (IF) model for describing the spiking activity of neurons, which is characterized by the dynamics of the membrane potential V(t). An IF neuron consists of a basic RC-circuit with the cell membrane capacitance Cm in parallel with a membrane resistance Rm. If the membrane potential is below a given threshold Vthr (subthreshold dynamics), then the membrane potential of each neuron in the network is given by the following equation:
where gm = 1/Rm is the membrane leak conductance, VL is the resting potential, and Isyn is the synaptic current. The membrane time constant is defined by τm = Cm/gm. When the voltage across the membrane reaches the threshold Vthr, the neuron generates a spike, which is then transmitted to other neurons, and then the membrane potential is instantaneously reset to Vreset and maintained there for a refractory time τref, during which the neuron is unable to produce further spikes.
Incoming input currents coming from connected neurons or from external inputs drive the membrane potential. Indeed, the spikes arriving at a given neural synapse produce an input to the neuron, which induces postsynaptic excitatory or inhibitory potentials essentially given by a low-pass filtering formed through the synaptic receptors. In our case, the total synaptic current is given by the sum of glutamatergic AMPA external excitatory currents (IAMPA,ext), AMPA (IAMPA,rec) and NMDA (INMDA,rec) recurrent excitatory currents, and GABAergic recurrent inhibitory currents (IGABA):
where
Here gAMPA,ext, gAMPA,rec, gNMDA,rec, and gGABA are the synaptic conductances, and VE, VI the excitatory and inhibitory reversal potentials, respectively. The dimensionless parameters wj of the connections are the synaptic weights. The NMDA currents are voltage dependent and they are modulated by intracellular magnesium concentrations. The gating variables sji(t) are the fractions of open channels of neurons and are given by:
The sums over the index k represent all the spikes emitted by the presynaptic neuron j (at times tjk). In Equations 7–11, τNMDA,rise and τNMDA,decay are the rise and decay times for the NMDA synapses, and τAMPA and τGABA the decay times for AMPA and GABA synapses. The rise times of both AMPA and GABA synaptic currents are neglected because they are short (<1 ms). The values of the constant parameters and default values of the free parameters used in the simulations are displayed in Table 1.
Neural and synaptic parameters for each local brain area
Each local attractor network contains 100 excitatory pyramidal neurons and 100 inhibitory neurons. The total cortical network has, therefore, 66 × 200 = 13,200 neurons. We use local attractor networks where neurons are organized into two sets of populations (Fig. 2), namely: an inhibitory population, and an excitatory population. The network is fully connected, meaning that each neuron in the network receives NE excitatory and NI inhibitory synaptic contacts. The connection strengths between and within the populations are determined by dimensionless weights wj. The recurrent self-excitation within each excitatory population is given by the weight w+, and within each inhibitory population is given by the weight w = 1. The connections between excitatory and inhibitory neurons have the weight w = 1.
All neurons in the network always receive an external background input from Next = 800 external neurons emitting uncorrelated Poisson spike trains. More specifically, and for all neurons inside a given population p, the resulting global spike train, which is still Poissonian, has a time-varying rate νextp(t), governed by
where tn = 300 ms, ν0 = 2.4 kHz, σν is the SD of νextp(t), and np(t) is normalized Gaussian white noise. Negative values of νextp(t) that could arise due to the noise term are rectified to zero. These input rate fluctuations represent the noisy fluctuations that are typically observed in vivo.
Mean-field reduction.
For the mesoscopic level of description, the mean-field approximation reduces the number of integration variables to one for each neural population (Brunel and Wang, 2001). Solving the mean-field equations provides the fixed points of the population firing rates, i.e., the stationary states of the populations after the period of dynamical transients. As this can be done much more quickly than integrating the full spiking model, scanning the parameter space to find a parameter set matching the experimental findings becomes feasible. Other mean-field approaches consider parameter dispersion (Assisi et al., 2005; Stefanescu and Jirsa, 2008) and different types of coupling including gap junctions and chemical synapses (Jirsa and Stefanescu, 2011; Stefanescu and Jirsa, 2011), which typically results in a higher-dimensional state vector and hence is computationally more expensive.
In the mean-field formulation the potential of a neuron is calculated according to:
where V(t) is the membrane potential, x labels the populations, τx is the effective membrane time constant, μx is the mean value of the membrane potential in absence of spiking and fluctuations, σx measures the magnitude of the fluctuations and η is a Gaussian process with absolute exponentially decaying correlation function and time constant τAMPA. The quantities μx and σx are given by:
where νext is the external incoming spiking rate, νI is the spiking rate of the inhibitory population, τm= Cm/gm with the values for the excitatory or inhibitory neurons depending of the population considered. The other quantities are given by:
where p is the number of excitatory populations, fx is the fraction of neurons in the excitatory population x, ωj,x the weight of the connections from population x to population j, νx is the spiking rate of the excitatory population x, γ = 0.28, β = 0.062, and the average membrane potential <Vx> has a value between −55 mV and −50 mV.
The mean-field approximation finally yields a set of n nonlinear equations describing the average firing rates of the different populations in the network as a function of the defined quantities μx and σx:
where φ is the transduction function of population x, which gives the output rate of a population x in terms of the inputs, which in turn depend on the rates of all the populations.
with erf(u) the error function and τrp the refractory period which is considered to be 2 ms for excitatory neurons and 1 ms for inhibitory neurons. To solve the equations defined by Equation 31 for all x, we numerically integrate Equations 32–34 and the differential equation below, whose fixed-point solutions correspond to solutions to Equation 31:
To find the possible fixed points that coexist for a given parameter set, Equation 35 has to be integrated for different initial conditions of population firing rates over a range of external inputs. Generally, the firing rates obtained by the mean-field approximation would be exact if the number of neurons was infinitely large and the unitary postsynaptic potentials elicited by presynaptic spikes were infinitesimally small.
With the conductivity values specified in Table 1, the values of w+ studied ranged between 1.1 and 1.7, such that the local brain area network, when isolated from the rest of the cortex, has a unique stable attractor, namely the spontaneous state attractor, where all excitatory neurons spike asynchronously at a rate of 3 Hz and the inhibitory neurons at a rate of 9 Hz (Brunel and Wang, 2001). We also studied the effect of different levels of noise from σν = 0.1–0.4. The results are robust and not dependent of the values of w+and σν. In the Results section we show the simulations by using always w+ = 1.5 and σν = 0.14.
BOLD-fMRI signal.
The simulation of the fMRI BOLD signal in the global cortical model is computed by means of the Balloon-Windkessel hemodynamic model of (Friston et al., 2003). The Ballon-Windkessel model describes the coupling of perfusion to BOLD signal, with a dynamical model of the transduction of neural activity into perfusion changes. The model assumes that the BOLD signal is a static nonlinear function of the normalized total deoxyhemoglobin voxel content, normalized venous volume, resting net oxygen extraction fraction by the capillary bed, and resting blood volume fraction. The BOLD-signal estimation for each brain area is computed by the level of neuronal activity summed over all neurons in both populations (excitatory and inhibitory populations) in that particular area. In all our simulations shown here this level of neuronal activity is given by the rate of spiking activity in windows of 1 ms. In brief, for the ith region, neuronal activity zi causes an increase in a vasodilatory signal si that is subject to autoregulatory feedback. Inflow fi responds in proportion to this signal with concomitant changes in blood volume vi and deoxyhemoglobin content qi. The equations relating these biophysical variables are:
where ρ is the resting oxygen extraction fraction. The BOLD signal is taken to be a static nonlinear function of volume and deoxyhemoglobin that comprises a volume-weighted sum of extra- and intravascular signals:
where V0 = 0.02 is the resting blood volume fraction. The biophysical parameters were taken as in (Friston et al., 2003). We also tested simulations in which the activity zi is given by the glutamate turnover (Logothetis et al., 2001; Attwell et al., 2010) and obtained equivalent results, since the synaptic activity typically strongly correlates with the spiking activity in our network.
Entropy.
To explore how the dynamics depend on W (scaling of the global inter-areal coupling strength) the mean-field equations are solved for particular values of W starting from different initial conditions. These initial conditions are chosen at random for each node by initially setting the average firing rate of the selective excitatory pool randomly to 10 Hz (active) or 3 Hz (inactive) and the average firing rate of the inhibitory pool to 9 Hz. To determine at which value for W a bifurcation has occurred the number of attractors and the entropy of the system is derived from the final average firing rates of the selective pools after 4000 steps found for each initial condition. If the final average firing rate of each selective pool is self-similar across initial conditions; i.e., the final average firing rate of each selective pool after initial condition X does not differ from any of the final average firing rates of this same pool by >10 Hz, there is only a single attractor. If this is not the case, the number of attractors can be established by counting the number of initial conditions in which at least one selective pool has a different final average firing rate from all previously found final average firing rates of that same pool. In the case of only a single attractor the entropy is 0 as the system will invariably settle within it. If, however, the number of attractors is larger than 1, the entropy is given by
where p(i) is the probability that the system settles in attractor i.
Results
We consider the dynamics of a realistic global attractor spiking network structured in brain areas, which are interconnected according the diffusion tensor imaging (DTI)/DSI-based neuroanatomical connectivity matrix of the brain of five human subjects (Human Connectome; Hagmann et al., 2008; see Materials and Methods). Precisely, the interplay between the particular underlying structure of the neuroanatomical inter-areal connections, the local attractor spiking dynamics of each area and the noise is able to explain the generation of the spatiotemporal structured resting state networks evidenced by the fMRI functional connectivity gained from the same subjects. In the following we will demonstrate that resting state networks in fMRI result from structured noise fluctuations around the trivial low firing equilibrium state induced at the edge of a bifurcation by the presence of latent “ghost” multistable attractors corresponding to distinct foci of high firing activity in particular brain areas.
We first study the stationary fix points (attractor landscape) of the cortical spiking network described in the Materials and Methods. We reduce the spiking dynamics to a set of mean-field equations that describes the stationary states (see Materials and Methods). In particular, we study the attractors of the global cortical system as a function of the parameter W scaling the global inter-areal coupling strength. Figure 3 plots the number of attractors found in the system as a function of W. We find these attractors by iterating as usual the mean-field reduced rate equations initialized with 1000 different initial conditions. Furthermore, Figure 3 also shows the entropy of the attractors calculated by considering the probability of obtaining the different attractors after departing from random initial conditions (see Materials and Methods).
Mean-field analyses of the attractor landscape of the cortical spiking network as a function of the global inter-areal coupling strength. The dashed line plots the number of stable attractors, whereas the continued line shows the entropy of the attractors.
This entropy characterizes the effective expected variability of cortical activity due to nonstationary noise-driven transitions between multistable attractors (which will be studied later with the full spiking implementation). For very small values of W only one attractor is stable and therefore the entropy is zero. That attractor corresponds to the trivial spontaneous ground state of the system where all neurons in all brain areas are firing at a low level of activity (excitatory neuron at 3 Hz and inhibitory neurons at 9 Hz). For very large values of W also only one attractor is stable and therefore the entropy again turns to zero. That attractor corresponds to the “epileptiform” case where all excitatory neurons are highly activated in all brain areas. In intermediate regions of W multistability of many attractors corresponding to distinct foci of high firing activity in particular brain areas emerges, causing consequently an increase of the entropy. For our further analysis, the bifurcation separating the trivial spontaneous state and the multistability region at W = Wc = 1.6 will be particularly relevant.
To identify the mechanisms underlying resting state generation, we study next how the dynamics of the model depends on the global coupling strength W. The dynamics of the cortical model is simulated by using the spiking model. We first calculate the firing activity in all brain areas, and then simulate the BOLD-fMRI signal by using the Balloon-Windkessel model (Friston et al., 2003; see Materials and Methods). Then the simulated BOLD signal was downsampled at 2 s to have the same resolution as in (Honey et al., 2009) and the global signal (average over all regions) (Fox et al., 2005, 2007) was regressed out of the BOLD time series. Finally, we computed the simulated functional connectivity by calculating the correlation matrix of the BOLD activity between all brain areas. To identify the region of the parameter W where the model best reproduces the empirical functional connectivity, we computed the Pearson correlation between both the empirical and the simulated functional connectivity matrix. Figure 4 shows how the fitting of empirical data as measured by the correlation between simulated and empirical functional connectivity is maximal for values of W at the edge of the bifurcation found in Figure 3. Interhemispherical correlations are much weaker and sometimes missed in the model, because the DTI/DSI tractography also missed those connections in the neuroanatomical matrix that we used (see Fig. 1a). Figure 4a plots the fitting of the functional connectivity at rest between theory and empirical data taking into account all pairs, i.e., including also the interhemispherical, whereas Figure 4b plots the results for single hemispheres, i.e., excluding interhemispherical pairs.
Fitting of empirical data as measured by the correlation between simulated and empirical functional connectivity as a function of the global coupling parameter W. The best fit is achieved at the edge of the bifurcation.
For the particular value W = Wc, Figure 5 shows a detailed comparison between the simulated and empirical functional connectivity. Figure 5a shows the simulated (left and middle subpanel) and empirical (right subpanel) functional connectivity matrices. In particular, for the simulated model we present the functional connectivity matrix based on the firing rate dynamics (left subpanel) and on the simulated BOLD signal (middle subpanel). The highest empirical functional connectivity correlations correspond mainly to pairs linked by the neuroanatomical connectivity matrix (Honey et al., 2009). The best agreement occurs when the simulated correlation is also high for these pairs. Figure 5b presents the neuroanatomical connectivity matrix (left subpanel), the empirical (middle subpanel) and the simulated (right subpanel) functional connectivity between one seed region and all other brain regions. We took the left posterior cingulated (lPC) as a seed (red bars in the plots). The simulated functional connectivity largely reproduces a lot of details of the empirical functional connectivity between the Default Mode Network regions (Greicius et al., 2003) along the medial axis. Figure 5c shows that Pearson correlation between both functional connectivity matrices for a given seed is significantly positive for all seeds.
Detailed comparison between the neuroanatomical connectivity matrix, the empirical and the simulated functional connectivity for the working point W = Wc = 1.6 at the edge of the bifurcation. a, Left, Functional connectivity matrix based on the firing rates dynamics. Middle, Functional connectivity matrix based on the simulated BOLD signal. Right, Empirical functional connectivity. b, From left to right: neuroanatomical connections, empirical and simulated functional connectivity for the seed lPC. c, Pearson correlation between empirical and simulated functional connectivity for each individual seed (blue), and the corresponding p-value (red).
For the same working point, at the edge of the bifurcation, Figure 6 shows typical temporal evolution of the simulated BOLD signal for three particular brain regions. We took two regions (lPC and left superior frontal (lSF)) in the Default Mode Network which are correlated and a third region in temporal (lST) brain areas (Attentional Network) which is anticorrelated with the first two regions. The model is able to capture the correlation of BOLD activity between brain regions in the same network and the anticorrelation between brain region in different resting state networks as evidenced previously by Fox et al. (2005). Further more, these patterns of spatiotemporal correlations defining the concept of resting state networks emerge at the level of slow fluctuations ∼0.1 Hz.
Simulated BOLD signal for 3 brain regions (lPC, green; lSF, blue and lST, red). lPC and lSF are part of the Default Mode network, whereas lST is part of the Attentional Network. Both in the empirical and in the simulated data, the functional connectivity between lPC and lSF is positive, while negative with lST. The temporal evolution of the signal shows this pattern of correlations.
At last, we show that at the edge of the bifurcation the tail of the distribution of the pair correlations between brain areas is maximally consistent with a power law distribution. To test the hypothesis that a given dataset is consistent with the assumption that they are drawn from a power-law distribution, we follow the standard goodness-of-fit test, which generates a p-value that quantifies the plausibility of the hypothesis. This test basically compares the “distance” between the distribution of the given dataset and the hypothesized power law distribution. Concretely we use here the method of Clauset et al. (2009) which first fit the given dataset to the power-law and calculate the Kolmogorov–Smirnov (KS) statistic for this fit. Next, we generate a large number of power-law distributed synthetic datasets with parameters equal to those of the distribution that best fit the given dataset. Each synthetic dataset is fit individually to its own power-law model and a KS statistic is obtained for each one relative to its own model. At last we count what fraction of the time the resulting statistic is larger than the value for the given dataset. This fraction establishes the p-value. If p is large, then the difference between the given dataset and the assumed underlying power law model can be attributed to statistical fluctuations alone; if it is small (<0.1), the power law model is not a plausible fit to the data. Figure 7 plots the p-values of the Kolmogorov–Smirnov test, the goodness of fit (gof), and the estimated maximum-likelihood estimate of the scaling exponent (alpha), as a function of the global coupling parameter W. As before Figure 7a plots the results taking into account all pairs, i.e., including also the interhemispherical, whereas Figure 7b plots the results for single hemispheres, i.e., excluding interhemispherical pairs. The horizontal line in each subpanel corresponds to the value obtained for the empirical data. A maximum of the p-value at the edge of the bifurcation is observed, indicating that at this point characteristic criticality behavior, like power law distributions, is observed. This is consistent with the recent analyses of Kitzbichler et al. (2009), who demonstrated criticality in empirical resting state data. Further more, the p-value, gof, and alpha estimated scaling exponent of the Kolmogorov–Smirnov test on the empirical functional connectivity data also gives a significant value and are very similar to the ones obtained in the simulations precisely at the brink of the bifurcation.
p-values of the Kolmogorov–Smirnov test, the gof, and the estimated maximum-likelihood estimate of the scaling exponent (alpha), as a function of the global coupling parameter W. a plots the results taking into account all pairs, i.e., including also the interhemispherical, whereas b plots the results for single hemispheres, i.e., excluding interhemispherical pairs. The horizontal line in each subpanel corresponds to the value obtained for the empirical data.
Discussion
In this paper, we investigated the neuronal and cortical mechanisms underlying the generation of resting state networks in the brain in absence of any stimulation and task. We formulated and studied a detailed and realistic spiking cortical network that is microscopically organized as standard attractor models (known from neural models of memory, attention, decision making, etc.) and mesoscopically organized through a neuroanatomical large-scale connectivity matrix obtained from human subjects via DTI/DSI tractography. More specifically, the model of a local brain area consists of integrate-and-fire spiking neurons with excitatory (AMPA and NMDA) and inhibitory (GABA-A) synaptic receptor types. We demonstrated that the best fit of the empirical data as characterized by the fMRI BOLD signal-based functional connectivity is obtained when the brain network is critical. In other words, the brain network operates at the brink of a bifurcation that separates the stable equilibrium low activity state from the multistable state region where many attractors corresponding to high activity in different brain areas coexist. Under these conditions the slowly fluctuating (<0.1 Hz) resting state networks emerge as structured noise fluctuations around the stable low activity state induced by the presence of latent “ghost” multistable attractors at the edge of the bifurcation.
Conceptually the notion of a global attractor model of the brain network is not new (Grossberg, 1980; Hopfield, 1984; Carpenter and Grossberg, 2003; Haken, 2004). In these instances the equilibrium of a network is typically multistable, in other words different equilibrium states coexist and one state is finally realized. Ghosh et al. (2008) proposed originally that the RSN settles in one of these states and then explores in the presence of noise the neighborhood of this particular state. They have clearly demonstrated that the brain network must operate close to instability of the equilibrium to allow for such exploration. This phenomenon is called criticality. However, Ghosh et al. (2008) have not established the coexistence of multiple states, which will be eventually visited during the noise-driven exploration. Furthermore, if the brain network dynamics is indeed critical, then the multiple states may be stable, but they also may represent latent “ghost” attractors, which are regions of the state space just at the edge of the bifurcation. Answers to these questions will shed light not only on the mechanisms giving rise to the emergence of RSNs, but also have implications regarding our understanding of brain function.
In contrast to previous models (Honey et al., 2007; Ghosh et al., 2008; Deco et al., 2009) the spiking attractor model used here does not exhibit oscillations explicitly, but describes the firing rate of pools of neurons. Modulations of the neuronal subthreshold activity giving rise to rhythms as observed in local field potentials can be studied in principle, but have not been considered here, since we focused upon the emergent fix point attractors and their role in shaping the brain network's dynamic repertoire. This focus upon fix point attractors does not imply a static nature of the resting state network dynamics. On the contrary, the dynamics of the asynchronously firing neurons within a brain area is complex, but the fix point attractors prescribe its firing rate. The advantage of neglecting the oscillatory physiological signals lies in the fact that we can neglect the time delays via signal transmission between brain areas. This is justified in this particular case, because time delays do not alter the stationary attractor states of a network, however, they may alter the stability boundaries of the oscillatory network states (for review, see Jirsa, 2004; Jirsa and Ding, 2004; Jirsa, 2009). As a consequence, the resulting dynamic repertoire will be mostly determined by the properties of the neuroanatomical connectivity matrix. Through this approach we demonstrated that a large number of resting state networks, and therefore a large number of genuine attractors beyond the bifurcation threshold, is crucial for the size and richness of a network's dynamic repertoire.
It is intellectually tempting (and not too far-fetched, we believe) to interpret the multistable attractor landscape beyond the bifurcation point as functionally meaningful. When choosing the working point of the brain network at rest just below the bifurcation, then the multistable attractors express themselves as latent ghost attractors, i.e., the stable fixed points do not exist yet, but are either saddle points with at least one unstable direction, or regimes in the state space with slow flow. Criticality here means that the magnitude of the flow is close to zero at these points. As a consequence, these states can be easily stabilized when needed in a given task context or for a given function. This is in contrast to existing approaches in cognitive neurosciences, in which cognitive/behavioral states coexist and (task-specific) stimuli drive the system from one state to the other, or alternatively, state transitions are considered as mechanisms for goal-directed behavior (Kelso, 1995; Calvert et al., 2004; Haken, 2004; Kelso and Jirsa, 2004; Rolls and Deco, 2010). Our here proposed scenario is different: there exists a set of available brain states in a repertoire that can be rapidly activated. Even in absence of any task, the repertoire exists and expresses itself by shaping the dynamic flow in the neighborhood of the equilibrium point of the network. This subtle but crucial element establishes also the main difference of the here presented work and the model of Ghosh et al. (2008) with other RSN models. It is the fluctuations around the equilibrium point of the network that are structured, rather than the network state. In resting conditions, the system does not hop from a multistable state to another multistable state like in former models (Deco et al., 2009) but the fluctuations around the unique equilibrium (the low activity ground state) become structured induced by the attractor structure of the cortical model above the bifurcation point. The mechanism as proposed by Ghosh et al. (2008) is identical to the here proposed structured fluctuations, but it has never been established that the attractor landscape of Ghosh et al. (2008) is indeed multistable. This crucial multistability, however, we have demonstrated here explicitly. At the working point just below the critical threshold value, the brain network is not particularly prone to be in any other state than the equilibrium state, but with gentle external stimulation the system could be rapidly stabilized in one of the attractor states. The stabilization of a particular attractor will require well controlled, but minor reconfigurations of the network parameters rather than creating an attractor entirely from scratch. As we demonstrated here, the attractors are encoded in the connectivity matrix and are presumably associated with the computation of a specific brain function, which is evidenced by the fact that the resting state networks resemble intermittent activations of spatial network patterns that are consistently known to be activated under cognitive and behavioral task conditions (Fox et al., 2005). There have been several other studies, which demonstrated that in the absence of an overt task, spontaneous fluctuations in the BOLD signal correlate across functionally related brain regions (Lowe et al., 1998; Gusnard and Raichle, 2001; Greicius et al., 2003; Rogers et al., 2007). Some of these regions are also part of the Default Mode Network (Greicius et al., 2003), which identifies regions showing the greatest deactivation during externally imposed cognitive challenges. Interestingly multistable attractors (more accurately, bistable limit cycles) have been also postulated to be involved in the generation of the alpha rhythm (Freyer et al., 2011), which is a prominent mostly posterior rhythm of 8–12 Hz in the human. Alpha activity arises in absence of any stimulation or task context, in particular with eyes closed. Under conditions of noise, the corticothalamic model of Freyer et al. (2011) operates close to a subcritical Hopf bifurcation and captures various detailed temporal characteristics of alpha activity, including the bimodal distribution of the power spectrum and scale-invariant fluctuations. In its spirit, this modeling approach is closer to Deco et al. (2009) since it places the ongoing cortical oscillations beyond the threshold onto self-sustained oscillatory attractors (limit cycles). However, no spatiotemporal extension of this effort exists so far. Nevertheless, the detailed elaboration of the temporal characteristics of the alpha spectrum and the authors' conclusion that the brain operates close to instability is intriguing. In other words, the brain at rest appears to be critical.
We conclude that the multistable attractor landscape, as expressed by the latent ghost attractors during resting state conditions, defines a functionally meaningful dynamic repertoire of the brain network, that is inherently present in the neuroanatomical connectivity. This repertoire is functionally relevant, since it allows the system to rapidly compute a specific brain function through stabilization of one of its attractors. Consequently, the more entropy of attractors exists, the richer is the dynamical repertoire and therefore the brain network displays more capabilities of computation. We speculate therefore that the neuroanatomical connections developed a scale-free type of architecture to be able to store a large number of different and flexibly accessible brain functions. The numerous brain functions are evidenced indirectly under resting state conditions by the generation of a large diversity of networks reflecting different ways of structured fluctuations, i.e., by the resting state networks. These mechanisms will be studied in future works. An advantage of spiking models is precisely that we can incorporate effects like oscillations, neuromodulations, etc., in a simple way and without destroying the structure of attractors.
Footnotes
G.D. was supported by the European Union Grant EC005-024, by SAF2010-16085 and the “La Marato” Foundation, and by the CONSOLIDER-INGENIO 2010 Programme CSD2007-00012. The research reported herein was supported by the Brain Network Recovery Group through the James S. McDonnell Foundation and the FP7-ICT BrainScales.
The authors declare no competing financial interests.
- Correspondence should be addressed to Gustavo Deco, DTIC-Universitat Pompeu Fabra, C/Roc Boronat, 138, 08018 Barcelona, Spain. gustavo.deco{at}upf.edu
