Abstract
Brain fluctuations at rest are not random but are structured in spatial patterns of correlated activity across different brain areas. The question of how resting-state functional connectivity (FC) emerges from the brain's anatomical connections has motivated several experimental and computational studies to understand structure–function relationships. However, the mechanistic origin of resting state is obscured by large-scale models' complexity, and a close structure–function relation is still an open problem. Thus, a realistic but simple enough description of relevant brain dynamics is needed. Here, we derived a dynamic mean field model that consistently summarizes the realistic dynamics of a detailed spiking and conductance-based synaptic large-scale network, in which connectivity is constrained by diffusion imaging data from human subjects. The dynamic mean field approximates the ensemble dynamics, whose temporal evolution is dominated by the longest time scale of the system. With this reduction, we demonstrated that FC emerges as structured linear fluctuations around a stable low firing activity state close to destabilization. Moreover, the model can be further and crucially simplified into a set of motion equations for statistical moments, providing a direct analytical link between anatomical structure, neural network dynamics, and FC. Our study suggests that FC arises from noise propagation and dynamical slowing down of fluctuations in an anatomically constrained dynamical system. Altogether, the reduction from spiking models to statistical moments presented here provides a new framework to explicitly understand the building up of FC through neuronal dynamics underpinned by anatomical connections and to drive hypotheses in task-evoked studies and for clinical applications.
Introduction
Ongoing brain fluctuations at rest (i.e., in the absence of external stimulation and response demands) are not random but are highly structured. Evidence from numerous neuroimaging and electrophysiological experiments demonstrates that the brain at rest displays spatial patterns of correlated activity across different brain areas known as resting-state networks (RSNs) (Biswal et al., 1995; Greicius et al., 2003; Fox et al., 2005; Fox and Raichle, 2007; Fransson, 2005; Deco and Corbetta, 2011). It has been demonstrated that the correlation structure of spontaneous BOLD fluctuations (i.e., the resting functional connectivity [FC]) relates to the underlying anatomical circuitry as obtained by diffusion tensor spectrum imaging (DTI/DSI) (Honey et al., 2007; Hagmann et al., 2008). Similar results were obtained in the monkey by combining fMRI and retrograde tracer maps (Vincent et al., 2007). This suggests that RSNs arise from correlations of neuronal noise between brain areas that are coupled by the underlying anatomical connectivity. Nevertheless, how FC relates to the anatomical connectivity and brain dynamics remains an open question and has motivated several experimental and network-modeling studies to understand structure–function relationships (Honey et al., 2007, 2009; Ghosh et al., 2008; Bullmore and Sporns, 2009; Deco et al., 2009, 2011; Cabral et al., 2011; Freyer et al., 2011; Pernice et al., 2011; Deco and Jirsa, 2012; Robinson, 2012).
One important conclusion of previous computational studies is that the relation between anatomical structure and functional correlations is strongly dependent on the local dynamics and the global dynamical state of the network. Thus, it is fundamental to describe the local dynamics of the different brain areas as realistically as possible. In particular, a recent study formulated a detailed and realistic large-scale model of global resting brain activity considering local networks of excitatory and inhibitory populations of spiking neurons coupled through NMDA, AMPA, and GABA synaptic dynamics (Deco and Jirsa, 2012). Within this model, the emerging RSNs best account for the FC measured experimentally using fMRI when the dynamical state of the network is at the edge of a critical point. Nevertheless, the spiking model, although biophysically realistic, is very difficult to treat analytically because of its complexity and its computational cost. Hence, there is a need for reducing the spiking model into an efficient simplified model that, while keeping the fair description of dynamics, allows for analytical insights.
In this study, we used a dynamic mean field technique (Wong and Wang, 2006) that approximates the average ensemble behavior, instead of considering the detailed interactions between individual neurons. Using this reduction, we studied the dynamics around the ground state, in which the firing rate in all areas is low, and the contribution of linear terms. Linearization enables further reduction of the originally stochastic dynamical system to its first- and second-order statistics, by deriving deterministic motion equations for means and covariances. This allows for parameter explorations, which are otherwise impractical because precise estimation of the different statistics requires extensive simulations of the entire stochastic system. More importantly, the reduction from the spiking large-scale model to statistical moments' equations leads to a compact relationship between correlation structure, underlying connectivity, and dynamics.
Materials and Methods
Anatomical connectivity.
The anatomical structure of the global brain model is defined by the connections between different cortical areas extracted from healthy humans using diffusion spectrum imaging (DSI) techniques. Here, we used the DSI data of five healthy right-handed male human subjects (Hagmann et al., 2008; Honey et al., 2009). The gray matter was first parcellated in 33 cortical areas (Table 1) per hemisphere (66 areas in total) according to anatomical landmarks. This was followed by a further parcellation into 998 ROIs. White matter tractography was used to determine which voxel pairs were connected by putative white matter fiber tracts and to estimate their density and corresponding lengths, from which the anatomical connectivity matrix was obtained at the ROI level and then averaged across five subjects. Because tractography does not give fiber directionality, the connectivity matrix is symmetric at the voxel level. To downsample the anatomical connectivity to 66 regions, the connection strength between two regions was calculated by summing all incoming fiber strengths to the target region, and dividing it by its region-dependent number of ROIs, resulting in a nonsymmetric downsampled connectivity matrix. This normalization by the number of ROIs, which have approximately the same surface on the cortex (i.e., the same number of neurons) is required because neuronal activity is sensitive to the number of incoming fibers per neuron in the target region. As the dynamical model of one region already takes into account the effect of its internal connectivity (see below), the connection of a region to itself was set to 0 in the connectivity matrix for the simulations.
Names and abbreviations of the brain regions considered in the human connectome from Hagmann et al. (2008)a
fMRI FC.
The global dynamics of the model was studied by contrasting the simulated emerged resting FC with the empirical one. The empirical resting FC was measured in 24 right-handed healthy young volunteers (15 females, age range 20–31 years). Subjects lay in a supine position and viewed a black screen with a centered red fixation point of 0.3 visual degrees, through a mirror tilted by 45 degrees. Each volunteer underwent two scanning runs of 10 min in a resting-state condition. Specifically, they were instructed to relax but to maintain fixation during scanning. The eye position was monitored at 120 Hz during scanning using an ISCAN eye tracker system. Scanning was performed with a 3T MR scanner (Achieva; Philips Medical Systems) at the Institute for Advanced Biomedical Technologies in Chieti, Italy. Functional images were obtained using T2-weighted echo-planar imaging with BOLD contrast using SENSE imaging. Echo-planar imaging (TR/TE = 2000/35 ms) comprised 32 axial slices acquired continuously in ascending order covering the entire brain (voxel size = 3 × 3 × 3.5 mm3). For each scanning run, initial 5 dummy volumes allowing the MRI signal to reach steady state were discarded. The next 300 functional volumes were used for the analysis. A three-dimensional high-resolution T1-weighted image (TR/TE = 9.6/4.6 ms, voxel size = 0.98 × 0.98 × 1.2 mm3) covering the entire brain was acquired at the end of the scanning session and used for anatomical reference. Initial data preprocessing was performed using the SPM5 software package (Wellcome Department of Cognitive Neurology, London) running under MATLAB (MathWorks). The preprocessing steps involved the following (Mantini et al., 2007): (1) correction for slice-timing differences, (2) correction of head motion across functional images, (3) coregistration of the anatomical image and the mean functional image, and (4) spatial normalization of all images to a standard stereotaxic space (MNI) with a voxel size of 3 × 3×3 mm3. Furthermore, we used the GIFT toolbox (http://mialab.mrn.org/software/gift/index.html) to perform to a spatial-independent component analysis of the BOLD time series in MNI space. After data decomposition by independent component analysis, we detected artifactual independent components (ICs) on the basis of the following: (1) the ratio between the spatial correlations of the IC map with a ventricular CSF template and a gray matter template, respectively; and (2) the sparseness of the regions belonging to the IC map, assessed by combining spatial entropy and clustering degree measures (Sui et al., 2009). The removal of the artifactual ICs was performed using the GIFT toolbox. No global signal regression was performed. For each recording session (subject and run), we extracted the mean BOLD time series from the 998 ROIs of the brain atlas used by Hagmann et al. (2008) and calculated the 66 × 66 correlation matrix representing the mutual FC between each pair of cortical areas. To further select the sessions, we discarded those for which the correlation matrix showed a pattern of global correlation. Finally, we concatenated in time the remaining sessions for each parcel and calculated the correlation matrix.
Spiking brain model.
The global brain large-scale model is a network of networks, used previously (Deco and Jirsa, 2012). It is composed of local networks, corresponding to different cortical brain areas, interconnected in a large-scale network according to a structural connectivity matrix, noted C and defined by the neuroanatomical connections between those brain areas in the human as obtained by DSI and tractography (see above). Local networks consist of mutually interconnected populations of excitatory (pyramidal) and inhibitory spiking neurons. The interarea connections are established as long-range synaptic connections between pyramidal neurons in those areas and are weighted by the strength specified in the matrix C and by a global scaling factor G (free control parameter).
Local networks describing the dynamics of different local cortical areas consist of spiking neurons with excitatory (AMPA and NMDA) and inhibitory (GABA-A) synaptic receptor types (Brunel and Wang, 2001; Deco and Jirsa, 2012). Each of these local networks is a dynamical system specified by a large set of coupled equations describing each neuron and synapse. The stationary fixed points (stable patterns of firing activity) of this type of networks are called “attractors.” Each local cortical area is modeled by a fully connected recurrent network of a population of NE excitatory pyramidal neurons and a population of NI inhibitory neurons. Each neuron in the network is described by the classical Integrate-and-Fire model, which expresses the dynamical evolution of the neuron's membrane potential V(t) that is driven by incoming local excitatory and inhibitory inputs within the same cortical area, long range excitatory inputs from all other cortical areas, and external inputs. The time evolution of the membrane potential of a given neuron i in the network obeys the following differential equation:
if the membrane potential is below a given threshold Vthr. The neuron generates a spike, when the membrane potential reaches the threshold Vthr. The spike is transmitted to other neurons, and 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. In Equation 1, gm is the membrane leak conductance, Cm is the capacity of the membrane, and VL is the resting potential. The membrane time constant is defined by τm = Cm/gm. The synaptic input current is given by the last four terms on the right-hand side of Equation 1. The spikes arriving at the synapse produce a postsynaptic excitatory or inhibitory potential given by a conductance-based model specified by the synaptic receptors, modeled by the four last terms of the right hand side of Equation 1, corresponding to the following: glutamatergic AMPA external excitatory currents, AMPA and NMDA recurrent excitatory currents, and GABAergic recurrent inhibitory currents. The respective synaptic conductances are gAMPA,ext, gAMPA,rec, gNMDA,rec, and gGABA, and VE and VI are the excitatory and inhibitory reversal potentials, respectively. The dimensionless parameters wij of the connections are the synaptic weights defined as following: the recurrent self-excitation within each excitatory population is given by the weight w+ = 1.55 (for both AMPA and NMDA recurrent synapses); recurrent inhibition and connections between excitatory and inhibitory neurons have the weight wII = wEI = wIE = 1; the weight of connections between local cortical area interareal connections are proportional to the neuroanatomical matrix and given by wij = G.C. The global scaling factor G is a free control parameter that we will vary systematically to study the dynamics of the global cortical system. The gating variables sjI(t) represent the fractions of open channels for AMPA and GABA receptor-mediated currents is described as follows:
NMDA channels are represented by second-order kinetics as follows:
Where the sums over the index k represent all the spikes emitted by the presynaptic neuron j (at times tjk); τAMPA and τGABA represent the decay times for AMPA and GABA synapses, respectively, and τNMDA,rise and τNMDA,decay are the rise and decay times for the NMDA synapses, respectively. All neurons in the network receive an AMPA-mediated external background input of uncorrelated Poisson spike trains with a time-varying rate νextp(t) governed by the following:
where τn = 30 ms, ν0 = 2.18 kHz (equivalent to 800 excitatory connections from external neurons firing at a rate of 3 Hz), σν = 0.01 kHz is the SD of νep(t) and dW is a Wiener process. These Poisson spikes drive the evolution of sAMPA,ext through Equation 2. Negative values of νep(t) that could arise because of the noise term are rectified to zero. These input rate fluctuations represent the noisy fluctuations that are typically observed in vivo, arising from finite size noise effects of the spiking dynamics of the individual neurons with their Poisson-like spike trains in a system of limited size. Table 2 specifies all the values used in the simulations. The parameters of the spiking model are the same as in Deco and Jirsa (2012), except for three of them: in the present study, we slightly increased the recurrent self-excitation within each excitatory population (w+ = 1.55, instead of 1.5) and decreased both the intensity (ν0 = 2.18 kHz, instead of 2.4 kHz) and the noise amplitude (σν = 0.01 kHz, instead of 0.14 kHz) of the external background input. These modifications change the bifurcation diagram of the large-scale network: with the new parameters, we obtained a network for which the first (i.e., coexistence of multistable states of high activity and a state of low activity) and the second bifurcation (destabilization of the state of low activity) are separated (see Fig. 2A), whereas in Deco and Jirsa (2012) both bifurcations overlapped. This allows us to assess which of the two bifurcations is the relevant one for studying the emergence of FC.
Neural and synaptic parameters for each local brain area
Reduced brain model: dynamic mean field.
To simplify the mathematical and numerical description of the model, we derived a dynamical mean field reduction of each local spiking attractor (Wong and Wang, 2006). The dynamic mean field reduction is based on the classical mean field model (Brunel and Wang, 2001) but, importantly, it fairly approximates the temporal dynamics of the spiking network, contrary to classical mean field techniques, which are used to calculate the steady states only. In brief, the mean field approach considers the diffusion approximation according to which sums of synaptic gating variables (Eqs. 2 and 4) are replaced by a DC component and a Gaussian fluctuation term. Moreover, the first passage equation for calculating the firing rate is approximated by a simple sigmoidal input–output function (Abbott and Chance, 2005). Further reduction is made by, first, linearizing the input–output relation of the inhibitory interneurons, given that their mean firing rate is higher than the one of the excitatory population and typically falls between the range of 8–15 Hz, within which the single-cell input–output relation is almost linear (Wong and Wang, 2006). Furthermore, the synaptic gating variable of NMDA receptors has a much longer decay time constant (100 ms) than AMPA and GABA receptors. Therefore, the system can be further reduced by using an adiabatic approximation and assuming that the dynamics of the NMDA gating variable dominates the time evolution of the system, whereas the other synaptic variables instantaneously reach their steady state (Wong and Wang, 2006). Moreover, for simplicity here, we have neglected contributions by the AMPA receptors to local recurrent excitation. In conclusion, the whole dynamics of each local network of excitatory and inhibitory populations of spiking neurons interconnected via NMDA synapses can be expressed by a single one dimensional equation. The global brain dynamics of the network of interconnected local networks can be simply and consistently described by the following set of coupled nonlinear stochastic differential equations:
Where H(xi) and Si denote the population firing rate and the average synaptic gating variable at the local cortical area i (from 1 to N = 66 cortical areas in our case), respectively, w = 0.9 is the local excitatory recurrence, and Cij is the structural connectivity matrix expressing the neuroanatomical links between the areas i and j. Parameter values for the input–output function are a = 270 (n/C), b = 108 (Hz), and d = 0.154 (s). The kinetic parameters are γ = 0.641/1000 (the factor 1000 is for expressing everything in ms), and τS = 100 ms. The synaptic couplings are JN = 0.2609 (nA) and the overall effective external input is I0 = 0.3 (nA). The value of the parameters I0 was tuned to obtain a bifurcation diagram of the dynamic mean field model similar to the one obtained for the spiking model (see below). This tuning is necessary to get the system of Equations 6–8 to result from a succession of approximations, among them the neglect of AMPA receptors, which inclusion would otherwise highly complicate the reduced model (Wong and Wang, 2006). In Equation 6, υi(t) is uncorrelated standard Gaussian noise and the noise amplitude at each node is σ. Throughout this study, unless otherwise specified, we fixed the noise amplitude to σ = 0.001 (nA).
Linearization of the dynamical mean field model.
To test whether the resting-state fluctuations could be appropriately described by a linearization of the dynamics, we formulated and studied also a linear model derived from the first-order Taylor expansion of the dynamics mean field model described by Equations 6–8 around the stable stationary spontaneous state. The resulting model is defined by following system of linear equations:
where the supra-index (0) denotes the values of the corresponding variables for the stationary stable spontaneous fix point attractor (i.e., all values of the firing rate activity are low, ∼2–3 Hz). Of course, the linear Equation 9 is only valid if the spontaneous state is stable (i.e., until the bifurcation that describes the instability of the spontaneous attractor).
A semianalytical treatment of the global brain dynamics: the moments' method.
To estimate the network's statistics, we approximate deterministic dynamical equations for statistical moments of network's gating variables (Rodriguez and Tuckwell, 1996). This “moments' method” reduces the number of variables and allows for parameter exploration that, otherwise, is impractical because estimation of the moments would require extensive simulations of the entire stochastic system. For this, we express the system of stochastic differential Equations 6–8 in terms of the first- and second-order moments of the distribution of synaptic gating variables: μi, the expected mean gating variable of a given local cortical area i, and ρij, the covariance between gating variables of local cortical areas i and j. The moments are defined as follows:
where the angular brackets denote the average over realizations. To derive equations of motion of the moments, we used the Fokker-Plank equation for the distribution of gating variables, p(S⃗,t) where S⃗ = {S1,…,SN} which is given by the following:
where, Ak = −Sk/τs and Bk = (1 − Sk)γH(xk). Using Equation 13 and noting that 
And similarly:
Taylor expanding Si around μi, i.e. Si = μi + δSi, and keeping the terms up to <δSiδSj>, the equations of motion of the moments are derived:
where the matrix W is given by the following: Wij = wJNδij + GCij. In matrix form Equation 17 is as follows:
where the superscript T is the transpose, Qn is the covariance matrix of the noise, Pij = ρij, and J is the Jacobian matrix, given by the following:
Finally, using this method, we estimated the correlation structure between gating variables, defined as follows:
Stationary covariance and the link between structure, dynamics, and FC.
Equations 17 and 18 indicate that the covariance of gating variables is determined by both the underlying connectivity and the dynamics. To get further theoretical insights, in this section we show how the stationary covariance of spontaneous fluctuations explicitly relates to the underlying connectivity and dynamics. The covariance matrix of fluctuations around the stationary spontaneous state is given by the algebraic equation:
which can be solved using the eigen decomposition of the Jacobian matrix evaluated at the fixed points: J = LDL−1 where D is a diagonal matrix containing the eigenvalues of J, noted λi, and the columns of matrix L are the eigenvectors of J. For uncorrelated white noise, Qn is diagonal. Multiplying Equation 20 by L−1 from the left and by L−† from the right (the “†” being the conjugate transpose), we get the following:
where M is given by the following: Mij = − Q̃ïj/(λi + λ*j), and Q̃ = L−1QnL−†.
Thus, the covariance matrix of spontaneous fluctuations is determined by the eigenvalues of the Jacobian matrix, which, in turn, is related to the connectivity matrix and the dynamics (Eq. 19). Hence, Equation 21 provides a direct link between the correlation structure, the underlying connectivity, and dynamics. Indeed, the interpretation of Equation 21 is that the input covariance (M) is propagated through the dynamical system and is mapped to its approximated output (P).
Power spectrum of fluctuations.
The power spectrum of fluctuations around the fixed points, Φ(ω), can be obtained by, first, writing the equation of the fluctuations around the fixed points, given by the following:
where J is the Jacobian matrix calculated at the fixed points.
Then, tacking the Fourier transform of Equation 22:
Hence, δS̃(ω) = − (J + iω)−1συ̃(ω). The power spectrum being the expectation of the Fourier transformed autocorrelation function, we get:
Thus, the power spectrum of δSk is as follows:
Simulation of BOLD fMRI signals.
To constraint and to study the dynamical working point of the models described above, we use the empirical FC, defined as the correlation matrix of BOLD fMRI signals obtained in healthy humans during rest. To simulate BOLD fMRI signal, we used the Balloon-Windkessel hemodynamic model of Friston et al. (2003). The Balloon-Windkessel model describes the coupling of perfusion to BOLD signal, with a dynamical model of the transduction of neuronal 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 synaptic activity in that particular cortical area, noted zi for a given cortical are i. In the detailed spiking model, we considered the synaptic activity summed over the excitatory and inhibitory populations in each cortical area; and for the reduced model, we used the synaptic variable S. Briefly, for the i-th region, synaptic activity zi causes an increase in a vasodilatory signal xi 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 processes are as follows:
where ρ is the resting oxygen extraction fraction. The BOLD signal is given by the following:
where V0 = 0.02, k1 = 7ρ, k2 = 2, and k3 = 2ρ − 0.2. All biophysical parameters were taken as in Friston et al. (2003). The BOLD model converts the local synaptic activity of a given cortical area into an observable BOLD signal and does not actively couple the signals from other cortical areas.
The simulated BOLD signal was down-sampled at 2 s to have the same temporal resolution as in the empirically measured BOLD signal. Simulation length for computing the model FC was equivalent to 20 min.
Goodness of fit of model FC.
We measured the agreement between the empirical FC and the model FC by calculating the Pearson correlation coefficient between all corresponding pairwise correlation coefficients of both matrices. Because the sampling distribution of Pearson correlation rc coefficient is not normally distributed, we used the Fisher's z-transformation to convert rc to the normally distributed variable z, given by z = 0.5ln[(1 + rc)/(1 − rc)], before applying any subsequent test on correlation coefficients.
Results
Linking structure and dynamics through large-scale network modeling
Spatiotemporal patterns across the brain at rest are typically characterized by the resting FC in neuroimaging experiments. Theoretical models allow us to investigate explicitly the link between FC and the underlying anatomical connectivity (Fig. 1). The first component of resting-state models is the anatomical matrix that expresses the density of white matter fiber tracts connecting two different brain areas. Here, we used a 66 × 66 cortical human matrix extracted from the average of DSI-based tractography (Hagmann et al., 2008) obtained in five healthy human subjects (see Materials and Methods). The second component of the models is the type of dynamics assumed for the local nodes. Here, we assumed that a faithful description of the local brain area dynamics is indispensable for discovering the mechanisms underlying the experimentally observed resting global brain dynamics. Therefore, we described the local brain area dynamics by the detailed spiking and conductance-based synaptic model previously used by Deco and Jirsa (2012). This model is able to reproduce properly the typical asynchronous spontaneous low level of spiking activity at ∼3−10 Hz, as evidenced experimentally with neuronal recordings (Burns and Webb, 1976; Koch and Fuster, 1989; Wilson et al., 1994). The obtained spiking activity was transformed to BOLD fMRI signals using the Balloon-Windkessel hemodynamic model (Friston et al., 2003). Combining both structure and dynamics in the models, the next step is to evaluate the model predictions and compare them with empirical data (i.e., to identify the model parameters for which the model best fits the empirical resting FC). In this paper, the empirical FC was calculated by averaging the results obtained on 24 healthy human subjects (Hlinka et al., 2011) and projected to the same parcellation adopted for the anatomical structural matrix.
Linking anatomical connections and FC. Neuroanatomical connectivity data were obtained by DSI and tractography after averaging across subjects. Parcellation provides a connectivity matrix C linking the N cortical areas with clear anatomical landmarks. A neurodynamical model was then constructed using a set of stochastic differential equations coupled according to the connectivity matrix C. To validate the model, we compared the model spatiotemporal patterns to the ones observed in empirical data. In our case, the empirical FC was measured using fMRI BOLD activity. This framework enables us to study the link between anatomical structure and resting-state dynamics.
Dynamic mean field approximation
It has been shown that the detailed spiking model the best fit of empirical resting fMRI data is obtained when the brain network is critical (i.e., at the border of a dynamical bifurcation point) (Deco and Jirsa, 2012). Nevertheless, the model suffers from the disadvantage that its computational study involves the solution of typically many hundred thousands of coupled nonlinear differential equations corresponding to each single neuron and each single synapse. Therefore, for studying the global brain model in more detail, we reduced consistently the model by using dynamic mean field techniques (see Materials and Methods).
The reduced dynamic mean field model (DMF) ignores the interaction between single neurons within a cortical area and instead considers the ensemble dynamics, for which the temporal evolution is dominated by the longest time scale of the system. Specifically, because the synaptic gating variable (i.e., the fraction of open channels) of NMDA receptors has a much longer decay time constant (100 ms) than AMPA and GABA receptors, the system can be reduced by assuming that the dynamics of the NMDA gating variable dominates the time evolution of the system, whereas the other synaptic variables instantaneously reach their steady state. This leads to a reduced system of N coupled nonlinear differential equations (Eqs. 6–8), where N is the number of brain areas considered (in our case N = 66), expressing the dynamics of NMDA gating variables (see Materials and Methods).
We first calculated the bifurcation diagrams characterizing the stationary states of the brain system for both the spiking (Fig. 2A,C) and the reduced DMF model (Fig. 2B,D) as a function of the control parameter G, describing the scaling or global strength of the intercortical coupling (see Materials and Methods). We studied 1000 different random initial conditions to identify all possible stationary states. In the case of the spiking model, we used the classical mean field (Brunel and Wang, 2001; Deco and Jirsa, 2012) for studying the fixed points of the system. Figure 2A shows for the different possible stable states the maximal firing rate activity over all cortical areas. For small values of the global coupling G, only one stable state exists, characterized by a low firing activity in all cortical areas. For a critical value of G, a first bifurcation emerges where new multistable states of high activity appear while the state of low activity is still stable (Fig. 2C). For even larger values of G, a second bifurcation appears where the state of low activity becomes unstable. In this study, we slightly changed the external background input and the self-recurrence with respect to Deco and Jirsa (2012), to obtain a system with the first and the second bifurcations separated (see Materials and Methods). In the case of the reduced DMF model, we studied the fixed points of the system of equations (Eqs. 6–8), in the absence of noise. Here, we followed the reduction of Wong and Wang (2006) and used their proposed simplified version, which is obtained by neglecting the effect of AMPA receptors. This highly simplifies the mean field equations but also comes at the cost of departure of the mean field model from the spiking network model. For this reason, we tuned one parameter of the DMF model, the constant external input, to get a structure of the stationary attractors of the global system similar to the full spiking model (see Materials and Methods). The bifurcation diagram of the DMF model is shown in Figure 2, B and D. In this study, we assumed that the state of low activity represents the spontaneous state of the brain, which is assumed to be the resting state of the brain. In the following, we seek to identify the parameter regimen in which the emergent model FC matches the empirical one.
Bifurcation diagrams of full and reduced models. A, B, Mean-field analyses of the attractor landscape of the cortical spiking network (A) and the reduced DMF model (B) as a function of the global coupling strength. Each point represents the maximum firing rate activity among all nodes. In both cases, depending on the coupling strength (G) and the initial conditions, the network converges to one of three possible dynamical regimens: for low values of G, the network converges to a single stable state of low firing activity (i.e., the spontaneous state); for increasing values of G, new stable states of high activity together with the spontaneous state coexist; for even larger values of G, the spontaneous state becomes unstable. Insets, For a particular set of parameters in the multistability regimen (with spontaneous state), the set of possible stable states (colored dots). For visualization, the states (firing rate vectors of dimension, N = 66) are projected in the two first principal components of the state space using the PCA method. C, D, These states are represented by the specific configuration of firing rates shown for the spiking model (C) and its dynamic mean field reduction (D). These are as follows: the spontaneous state of low activity, plus three states for which some nodes have a high firing rate. spont., Spontaneous.
Empirical FC versus model FC
We calculated the FC for both the spiking model and the reduced DMF model. In the detailed spiking model, simulated BOLD signals were obtained by transforming the synaptic activity summed over the excitatory and inhibitory populations in each cortical area by the Balloon-Windkessel hemodynamic model (see Materials and Methods), whereas for the DMF the synaptic gating variables of each cortical area were used. The resulting correlation matrix of the BOLD activity between all brain areas was calculated to characterize the simulated resting FC. Here, we used as a measure of fitting the Pearson correlation coefficient between the Fisher z-transformed pairwise correlation coefficients of both the empirical and the simulated FC matrix (see Materials and Methods). Figure 3A plots this fit as a function of the scaling control parameter G. In the same plot, the second bifurcation is also shown (vertical line). The best fit of the empirical data is observed at the brink of the second bifurcation in both models. Increasing the noise amplitude in the DMF model shifts the bifurcation upon which the spontaneous state is destabilized toward lower values of G, but the best fit of the empirical data is still obtained near the instability. For strong noise, the agreement between the model FC and empirical FC is reduced. Figure 3B, D shows the comparison between pairwise correlations of both the optimal FC of the spiking model and DMF model and the empirical FC. The corresponding FC matrices are shown in Figure 3C–E. The model FC particularly missed the interactions between homotopic cortical areas from different brain hemispheres, which are prominent in the empirical FC (Fig. 3C,E, second diagonal of matrices). This is the result of the low quality of structural connectivity for long interhemispheric fiber tracts through the corpus callosum, which are not detected because of the limited resolution of current imaging/tractography techniques (Hagmann et al., 2008). These missed pathways have direct repercussions on the model results. Finally, the FCs obtained using the spiking and the DMF models are extremely similar. Indeed, the correlation coefficient between the elements of the FC matrices corresponding to the DMF model and the spiking model, shown in Figure 3C and Figure 3E, respectively, is equal to 0.80, meaning that the reduced DMF captures properly the relevant dynamical entrainments. This is crucial because the dimensionality reduction achieved using the DMF allows us to study much more exhaustively and in much shorter computational time the dynamics of the system. For example, the 28 points shown in Figure 3A for the spiking model took 1 month CPU on a powerful workstation, whereas the 270 points calculated with the reduced DMF took just a couple of hours.
Agreement between empirical and model FC. A, The fitting of empirical data as measured by the correlation between simulated and empirical FC, as a function of the global coupling parameter G, was calculated for both the spiking (red stars) and the reduced DMF models for two levels of noise (blue points: σ = 0.001, light blue trace: σ = 0.010). In the two models, the best fit is achieved at the edge of the second bifurcation (black vertical line) For σ = 0.010, the gray vertical indicates the critical value of G, for which the spontaneous state is destabilized. B, Comparison of pairwise correlation coefficients obtained using empirical BOLD activity and simulated BOLD activity from the spiking network. C, The subdiagonal elements represent the empirical FC, in warm color code, whereas upper-diagonal elements represent the FC generated by the spiking model (cold colors). D, Comparison of pairwise correlation coefficients obtained using empirical BOLD activity and simulated BOLD activity from the DMF model. E, The subdiagonal elements represent the empirical FC, in warm color code, whereas upper-diagonal elements represent the FC generated by the DMF model (cold colors). dyn., Dynamics.
Fluctuations around the spontaneous state
Given the fact that we were able to reduce significantly the model without losing the realism in the description of the dynamics, we can further investigate the nature of fluctuations and ask what is the contribution of linear terms in the exploration of the dynamics. The reduced DMF is simple enough (just one nonlinear equation per node) to linearize the equation by performing a first-order Taylor expansion around the global spontaneous fixed point (see Materials and Methods). Figure 4A compares the results of the nonlinear and the linearized DMF model by comparing again the fit of the BOLD simulations with the empirical resting fMRI data results. The linearized DMF model is only properly defined in the regions where the global spontaneous state of the full model is stable. The linearized version of the model fairly reproduces the resting global dynamics, meaning that a linearization is an excellent approximation.
Structured linear fluctuations around the spontaneous state. A, Fitting of empirical FC using simulated BOLD signals with the DMF model (black points) and its linearization (gray squares). B, The black trace represents the fitting between the empirical FC and the correlation structure of gating variables obtained using the derived equations for statistical moments (Eqs. 16 and 17). Gray points indicate the eigenvalues of the Jacobian matrix of the model for each value of G. C, Transfer function of the Balloon-Windkessel model, used for simulating the BOLD signals. The transfer function was computed as the amplitude of the obtained BOLD signal in response to cosine inputs of different frequencies, i.e., zi = cos(ωt). The BOLD model attenuates the high frequencies and allows low frequencies (<0.5 Hz) to pass. D, Power spectrum of gating variable fluctuations around the spontaneous state, Φ(ω), for three different cortical area (rBSTS, rCAC, and rCMF). The power spectrum Φ(ω) was obtained using Equation 25. Power of slow fluctuations increases with G.
Structure–function relation
Hence, we can further simplify the system of stochastic differential equations by expressing it in terms of the first- and second-order statistical moments (means and covariances) of network activity. In general, the nonlinear and stochastic nature of the original equations hinders analytical progress, and the main analysis of such stochastic differential equations remains based on numerical investigations, which are time consuming because of the need for sufficiently many trials to generate statistically meaningful data. But in the case of weak noise, one can derive deterministic differential equations for the first- and second-order moments of the distribution of the gating variables of the neuronal populations (see Materials and Methods). This is often referred to as the moments' method (Rodriguez and Tuckwell, 1996). The moments' equations for the reduced nonlinear DMF are derived in Materials and Methods. Briefly, this method uses the Fokker-Plank equation and Taylor expansion to derive motion equations for the means and covariances. Alternatively, equations for moments can be obtained using the Laplace transform (Marreiros et al., 2009).
Figure 4B shows the fitting between the empirical FC and the stationary correlation structure of gating variables obtained using the moments' method. In this case, the resting FC is not based on the BOLD signal, but on the correlation matrix between the gating variables of each cortical area, noted Q. The correlation structure Q is calculated semianalytically by solving the set of deterministic Equations 16 and 17. We found that the fitting between the model correlation structure and the empirical FC presents a maximum just before the (second) bifurcation, as for the correlation structure of BOLD from the DMF. However, the increase of the fitting in the former case is not as substantial as in the latter. This apparent discrepancy arises from the filtering effect of the BOLD model, which is not considered in the moment's equations. The Balloon-Windkessel model used to simulate the BOLD fMRI signal does not couple the signals of different brain areas but has an impact in the frequency domain. Indeed, the Balloon-Windkessel model can be seen as a low-pass filter: its transfer function passes frequencies ∼<0.5 Hz but attenuates higher frequencies (Fig. 4C) (Robinson et al., 2006). Thus, only gating variable correlations at slow time scales are transmitted through the BOLD model. Slow fluctuations are expected close to the edge of stability because the time scale of dynamics increases when the spontaneous state approaches the instability, for which the real part of some eigenvalues is close to zero (Fig. 4B, gray points). Indeed, far from the bifurcation, the spontaneous state is strongly attractive and thus fluctuations relax rapidly, with an exponential decay determined by the associated real part (negative) of the eigenvalues; as the global coupling increases, the attraction toward the spontaneous state weakens and the real part of some eigenvalues approaches zero, making the decay long, thus slowing down the time scale of dynamics; beyond the bifurcation, fluctuations grow exponentially, making the spontaneous state unstable. To quantify this effect, we considered the power spectrum of fluctuations around the spontaneous state, which can be derived by taking the Fourier transform of the equation governing the time evolution of fluctuations (see Materials and Methods). As shown in Figure 4C, the power of low frequencies (<1 Hz) increases when the global coupling increases. In conclusion, slow and structured fluctuations emerge when the spontaneous state starts losing stability, so that the correlation structure of gating variables is seen through the hemodynamically filtered response.
Furthermore, the deterministic equation for the covariance derived using the moments' method gives further theoretical insights about the relationship between anatomical structure, dynamics, and FC. Indeed, as shown in Materials and Methods, the propagation of noise through the network is determined both by the anatomical connectivity and dynamics, through the Jacobian matrix of the linearized system. More explicitly, Equation 21 expresses the contribution of both the anatomical structure and dynamics in the mapping of the intrinsic noise into the resulting functional correlations. Indeed, this mapping is expressed by the eigen decomposition of that Jacobian matrix, which, in turn, is explicitly determined by the connectivity matrix and the network dynamical state at which it is evaluated. Altogether, these results indicate that FC arises from reverberation of noise in the linearized dynamical system, which is constrained by the anatomical connections and that is subject to a dynamical slowdown because of the loss of stability of the spontaneous state.
Seed-based correlations and spatial modes
We next studied in more detail which features of the FC are explained by the model. For this, we first concentrated on seed-based correlations and principal components of the FC. Figure 5 plots the neuroanatomical connectivity matrix, the empirical FC, and the model FC between one seed region and all other brain regions at the best operating point (i.e., at the edge of the second bifurcation) for the spiking, the reduced DMF, and the moments-based models. We took the left posterior cingulate (lPC) as a seed, which is part of the well-known “default-mode network” (Shulman et al., 1997; Raichle et al., 2001; Snyder and Raichle, 2012). The simulated FC of all three models largely reproduces the empirical FC. In conclusion, the linearized reduction in form of the moments' method is an analytical model that captures excellently the dynamics of resting-state activity in an extremely simple way. In particular, note the excellent agreement between the correlation structure obtained using the moments' method and the one using both DMF and spiking models.
Comparison between empirical and model seed-based correlations. Comparison between the neuroanatomical connectivity matrix, the empirical FC, and the simulated FC. The model FC was calculated for the working point for which the model FC best fits the empirical FC (at the edge of the bifurcation). The first column represents the neuroanatomical connections linking the cortical area lPC. From the second to the fifth columns: FC for the seed lPC obtained from the empirical data, the spiking model, the DMF, the and moments' equations, respectively. Gray bars indicate the seed.
A further detailed comparison between the empirical and model FC can be obtained using the principal components, or “spatial modes,” of empirical and model spatial modes. The spatial modes are the eigenvectors of the covariance matrix. Figure 6A shows the first four dominant spatial modes of both the empirical covariance matrix and the covariance obtained using the moments' method, calculated for the working point for which the FC of the model best fits the empirical one, as previously done. The first two empirical and model principal components are similar (vector projection > 0.7). These two first modes contribute together to 27.8% of the total variance (Fig. 6B). Figure 6C shows the projection of the spatial modes of the different models into the empirical spatial modes. All models explain the emergence of the two first principal components, whereas the empirical structural connectivity alone matches the first principal component of the empirical FC. As shown in Materials and Methods, spatial modes of spontaneous fluctuations are related both to the dynamics and the structural connectivity. Indeed, how fluctuations around the spontaneous state are structured by the dynamical model (i.e., connectivity and dynamics) and how the linear approximation of the covariance of the fluctuations can be mapped are indicated by Equation 21. Finally, there is high projection along the diagonal between the spatial modes of spontaneous activity obtained using the moments' method and the one using both DMF and spiking models.
Comparison between empirical and model principal modes. A, The first four principal components (PCs) of both the empirical covariance matrix of BOLD signals (black) and the covariance matrix of gating variables of local cortical areas given by the moments' method (gray). B, Percentage of variance of the empirical PCs (ordered). C, Left, Projection (scalar product) between empirical PCs and the PCs of the different models. Right top, Projection between empirical PCs and the PCs of the anatomical connectivity. Right middle, Projection between PCs of the covariance of BOLD signals obtained using the dynamical mean field model and the PCs of the covariance matrix of gating variables of local cortical areas given by the moments' method. Right bottom, Projection between PCs of the covariance of BOLD signals obtained using the spiking model and the PCs of the covariance matrix of gating variables of local cortical areas given by the moments' method. Only the first 15 PCs are considered.
Exhaustive parameter exploration using the moments' method
The above results show that the moments' equations are an excellent approximation. We can thus benefit from the simplicity of the moments' method to explore the parameter space of the model. Using the moments' method, we calculated the correlation structure of the model and compared it with the empirical FC, for many different values of the global coupling (G) and the local excitatory recurrence (w) (Fig. 7). We found that the best match between the correlation structure of the model and the empirical FC is obtained near the bifurcation, for a wide range of parameters, indicating that fluctuations are relevantly structured near the critical point, independently of other parameters.
Landscape of fitting values in the G-w parameter space. The correlation structure of gating variables, Q, was calculated using the moments' method for many values of the global coupling (G) and the local excitatory recurrence (w). The resulting fitting (Pearson correlation) between Q and the empirical FC is shown in color code in the G-w plane. The green line indicates the second bifurcation, for which the spontaneous state becomes unstable. The red line indicates the maximum fit between Q and the empirical FC for a given w.
Dynamical changes of effective connectivity
Finally, the fact that the network's covariance arises from noise propagation through the (eigenvectors of) the Jacobian matrix, evaluated at the fixed point of the global network (Eq. 21), indicates that the covariance is essentially state-dependent. For instance, the application of external inputs/stimuli to a set of cortical areas changes the dynamical state of the brain network, to a new created fixed point, for which the Jacobian matrix is different; consequently, the correlations between areas are changed. To illustrate this and as a proof of principle, we applied an arbitrary input to the system and compared the resulting correlation structure of gating variables in absence (spontaneous or “rest”) and presence (“evoked”) of stimulation, respectively, obtained using the moments' method. Correlation matrices in spontaneous and evoked conditions, together with their corresponding firing rate state, are presented in Figure 8A. The difference between the two matrices (Δrc) is presented in Figure 8B. We found that application of an external input can substantially increase (Δrc > 0) or decrease (Δrc < 0) the pairwise correlations between cortical areas with respect to the spontaneous condition and therefore generates a different “effective” connectivity by changing the dynamics built on the anatomical connectivity.
Spontaneous FC versus evoked FC. A, The correlation matrix evoked by an external stimulus (upper triangular part of the matrix) was compared with the spontaneous correlation matrix (lower triangular part of the matrix); both matrices were calculated using the moment's method. In both conditions, the global coupling was G = 2.2, for which the agreement between the spontaneous correlation matrix and the empirical FC was maximum. In the evoked condition, the right and left lateral occipital cortex (LOCC) received an augmented constant input, equal to I0 + dI, with I0 = 0.3 (nA) and dI = 0.05 (nA), whereas all other cortical areas received the baseline input I0 (Eq. 8). The stationary states, at which correlations are calculated, in both spontaneous and evoked conditions, are represented by the histograms. Red bars indicated the stimulated nodes (rLOCC and lLOCC). B, Element-by-element difference (Δrc) between evoked correlation matrix and spontaneous (rest) correlation matrix.
Discussion
In this article, we demonstrated that resting FC emerges from slow linear explorations of the dynamical capabilities of the system at the brink of instability of the underlying global spontaneous state of the brain in which all brain areas spike at low level of activity. In other words, the resting state results from linear excursions at the edge of a nonlinear cliff (i.e., at the brink of a bifurcation), which are dynamically slowed down, resulting from loss of attraction, and shaped by propagation through the network architecture. We were able to show this by reducing the realistic coupled system of excitatory and inhibitory populations connected through conductance-based synapses tuned to express asynchronous low level spontaneous spiking activity. The reduction was made using a dynamic mean field approximation, describing in a simple way the ensemble dynamics, for which the time evolution is dominated by the dynamics of the slow NMDA gating variable. The reduced model fairly approximates the emergent FC of the full spiking model and is amenable to further reductions by linearization. Indeed, there exists evidence that the resting state of the brain is properly characterized by linear statistics, suggesting that a linear underlying dynamics may be sufficient (Hlinka et al., 2011). A further reduction based on the moments' method allowed us to derive a deterministic equation for the first- and second-order moments of the neuronal activity, such that the resting FC could be obtained very efficiently semianalytically and thereby avoiding long-lasting and multiple explicit simulations trials. This analytical simplification is essential because it allows a thorough investigation of the brain system and a more exhaustive optimization of multiple parameters. Moreover, the reduction into statistical moments provides an explicit link between structure, dynamics, and FC. We demonstrated analytically the contribution of both the anatomical structure and relevant dynamics in the propagation of the intrinsic noise into the building up of the FC.
Our results show that the different models (the spiking model, the dynamical mean field, and its linear approximations) fit well the first two principal components of the FC. This incomplete fitting is principally the result of the lack of interhemispherical correlations in the model because the DTI/DSI tractography also missed those connections in the neuroanatomical matrix that we used. Moreover, the anatomical matrix used here did not include subcortical routes that are known to play an important role in shaping the spontaneous activity of the brain (e.g., Robinson et al., 2001; Freyer et al., 2011). In addition, we have made several simplifying assumptions that also may impose an upper bound to the agreement between the model predictions and the empirical observations. Indeed, in this study, we consider that all connections between brain areas are excitatory and instantaneous, thus neglecting the effects of feedforward inhibition and conduction delays that are likely to shape the spatiotemporal patterns of brain dynamics. Indeed, it is known that time delays give rise to complex spatiotemporal patterns, oscillations, multistability, and chaos (Roxin et al., 2005; Ghosh et al., 2008; Deco et al., 2009).
A common characteristic of the present and previous models is that the optimal operating point for explaining the emergence of RSN is always near a critical point (i.e., a bifurcation). For example, Ghosh et al. (2008) showed that the resting state emerges from noise-induced transient fluctuations around the stable equilibrium state of a network of coupled FitzHugh-Nagumo oscillators at the edge of instability, whereas in the model proposed by Deco et al. (2009), the noise induces fluctuations between different multistable oscillatory brain states in a network of Wilson-Cowan oscillators, again at the edge of a critical point. In the study of Honey et al. (2007), RSNs also result from instability due to chaotic fluctuations. Deco and Jirsa (2012) and Robinson (2012) showed that the best agreement between model FC and empirical FC arises close to criticality. Moreover, the low-frequency 1/f behavior of EEG spectrum has been shown to be a signature of edge-of-stability dynamics (Robinson et al., 1997; Aburn et al., 2012). Furthermore, bistability of alpha rhythms is explained in a thalamocortical model also at the edge of a dynamical instability (Freyer et al., 2011). In other words, the type of local dynamics is relevant for determining the working point that generates resting FC. Nevertheless, in all models, the underlying anatomical structure shapes indeed in the same way as the dynamical landscape that is explored by the noisy fluctuations at its critical working point. This is the reason why, with more or less detail, all these models could explain the spatial functional correlations, which defines the different RSNs. In other words, at the edge of the critical instability of any model, the spatial correlations of the noisy excursions are mainly shaped by structure. Critical dynamics are functionally relevant. Indeed, working at the edge of a critical point allows the system to rapidly compute a specific brain function by representing it in an attractor. This may be a fundamental reason why RSNs reflect cognitive functions and why RSNs are so interesting for basic and clinical neuroscience.
The present study allows us to refine the interpretation made in Deco and Jirsa (2012) (i.e., that ghost attractors influence the structure of correlation). In this previous work, the parameters of the spiking network were such that the destabilization of the spontaneous state and the appearance of active attractors coincided. Here we showed that the loss of stability of the spontaneous state is the relevant instability for studying the emergence of structured correlations. Increasing the noise amplitude in the DMF model, so that ghost attractors influence more the dynamics, leads to a slightly lower agreement between the model FC and the empirical FC. This suggests that small fluctuations around the spontaneous state explain better the resting-state dynamics of the brain, against large excursions that are shaped by the landscape of attraction. However, the differences in fitting depend on the model and noise details and, thus, are not conclusive. Hence, future experimental investigation of the magnitude and the time scales of the empirical fluctuations may offer a way to distinguish between the two scenarios.
Our study shows that the network's covariance can be calculated from (eigenvectors of) the Jacobian matrix and thus indicates that the covariance is essentially state-dependent: the interactions between cortical areas depend on the dynamical state of the global network at which the Jacobian matrix is evaluated. Indeed, propagation through the eigenvectors of the Jacobian matrix can be seen as the effective connectivity (i.e., how different cortical areas influence each other, even in the absence of direct connection between them) (Aertsen et al., 1989; Bullmore and Sporns, 2009; Friston, 2011; Robinson, 2012). The key characteristic of effective connectivity is that it is not fixed because synapses in the brain are activity- and context-dependent, resulting from synaptic plasticity, neuromodulation, or even the interplay between excitation and inhibition (Vogels and Abbott, 2009). The effective connectivity has been efficiently estimated using dynamic causal modeling (Friston, 2011) and Granger causality (Brovelli et al., 2004). In particular, in the framework of dynamic causal modeling, Bayesian inversion of observed data is used to estimate the effective connectivity in task condition, whereas anatomical connectivity is used for constraining prior distribution of effective connectivity weights (Stephan et al., 2009). Complementary to this, here we used the anatomical connectivity to constrain a dynamical system modeling the brain at rest, and we proposed that, in task condition, changes of activity in a set of cortical areas resulting from external inputs, attention, or learning, change the dynamical state of the brain network, thus changing the Jacobian matrix and consequently modifying the correlations between areas. Therefore, we expect that future investigation on how specific localized inputs change the network effective connectivity will deepen our understanding of the relation between brain activity during task and rest conditions.
In this study, we focused on the FC for characterizing the resting state. The FC represents a static map of interactions between brain areas. This purely spatial measure is likely to be insufficient to fully explain the underlying dynamics, as temporal patterns are not considered. Thus, different theoretical models cannot be distinguished if only spatial correlations are considered. Nevertheless, different models might show different spatiotemporal dynamics, with spatial patterns forming, dissolving, and reforming across time, and we therefore speculate that the time dimension will be crucial. The interplay of time and space is critical for characterizing the resting state. Indeed, recent work indicates that interactions within- and between-networks are nonstationary (Hutchison et al., 2012; de Pasquale et al., 2012), indicating that other factors in addition to structural connectivity (e.g., time delays) must control the dynamics of spontaneous ongoing activity.
We envision a new way of looking at resting state, based on the theory, which can serve as fertile test ground for basic neuroscience and especially for clinical applications. The successive reduction from spiking models to statistical moments' equations presented here provide a theoretical framework to directly test hypotheses on how a particular underlying connectivity and dynamical properties structure the fluctuations in neural networks. For example, this approach can be used for understanding how an anatomical lesion or a pharmacological treatment affects resting FC. Indeed, an advantage of spiking models is that they can incorporate additional dynamics of specific neurotransmitter receptors that, for example, implement neuromodulation in a biophysically plausible way. Moreover, recently, Eckhoff et al. (2011) modified the dynamic mean field procedure to incorporate norepinephrine modulation of synaptic conductances in a spiking model, thus providing a realistic implementation of neuromodulation in a low-dimensional model. Therefore, we expect this approach to be essential for future investigations and especially for clinical applications.
Footnotes
This work was supported by the Brain Network Recovery Group through the James S. McDonnell Foundation. G.D. was supported by the European Research Council Advanced Grant DYSTRUCTURE (295129), the Spanish Research Project SAF2010-16085, the CONSOLIDER-INGENIO 2010 Program CSD2007-00012, the Foundation La Marato (Catalonia), and the FP7-ICT BrainScales. P.H. was supported by the Leenaards Foundation. M.C. was supported by National Institutes of Health Grants R01HD061117 and R01MH096482. We thank E. Hugues, R. Hindriks, and M. Adhikari for critical comments and valuable discussions.
The authors declare no competing financial interests.
- Correspondence should be addressed to Dr. Adrián Ponce-Alvarez, Universitat Pompeu Fabra, C/ Roc Boronat, 138, 08018 Barcelona, Spain. adrian.ponce{at}upf.edu
