Key role of coupling, delay, and noise in resting brain fluctuations

Brain's Intrinsic Properties.

The main aim of our investigation is to establish what particular intrinsic properties of brain networks play an essential role in the generation of the most typical aspects of brain dynamics at rest, namely slow oscillations and the emergence of anticorrelated subnetworks. In particular, we will consider 3 different intrinsic properties: (i) neuroanatomical connectivity structure, (ii) delays in the transmission of information between different brain nodes, and (iii) role of noisy fluctuations.

All simulations and analyses were performed by using a realistic connectivity matrix of the primate brain based on the CoCoMac neuroinformatics tool (14). Kötter and Wanke (16) proposed a coarse parcellation of cerebral cortex into 38 regions, which deliberately reflected broad and rather uncontroversial divisions so that a rough mapping to the human brain appeared feasible. For subsequent activation studies the regional map comprised in addition 2 subcortical thalamic regions, the pulvinar and anterior thalamic nucleus. Connectivity data from tracer studies collated in CoCoMac were transformed to the regional map by using the ORT procedure as described by Stephan et al. (17).

In addition, the center coordinates of the 38 cortical areas were calculated and their distances obtained from the geometry defined in the AAL cortical surface template of a human hemisphere (18). Assuming a uniform velocity of transmission v, we derived approximate delays. The velocity v is one of the free parameters that we consider in our parameter space study. The second parameter that we consider is the global coupling strength α between connected nodes (See Methods and SI for details).

The level of noisy fluctuation was also studied parametrically in the next section. We modeled random fluctuations using uncorrelated Gaussian noise that perturbed the population dynamics of each cortical node. Mathematically, this meant we simulated cortical activity by integrating stochastic differential equations based on a simple Wilson–Cowan model of population activity (see below). The origin of this noise could have different sources (see ref. 19). One realistic assumption might be spiking noise. Spiking fluctuations make a significant contribution, because this noise is a significant factor in a network with a finite (i.e., limited) number of neurons. It is important to note that these statistical fluctuations influence, on each trial, the dynamical characteristics of the outcome and not just its time course.

Collective Neurodynamics.

We consider in our simulations a very simple neurodynamic model for each node. We assume that each node's dynamics can be captured by a mean-field-like rate model expressing the coupling between excitatory and inhibitory neurons. In particular, we consider the Wilson–Cowan model, which is tuned such that each independent node, if disconnected, is silent (low-activity regime); but because their working point is very near to a Hopf bifurcation, when coupled, each node starts to oscillate. In particular, we choose a working point such that the oscillation of each node, which arises because of coupling, was in the γ-band-range of 40 Hz (see Methods and SI for details).

The reason for this choice is that we would like to keep the single node dynamics as simple as possible (oscillatory dynamics) and to concentrate our study on the emergence of a complex collective brain dynamics because of the intrinsic properties mentioned above. Furthermore, by considering simple 40-Hz fast oscillations at the single-node dynamics, we are able to investigate the link between fast local dynamics and slow global fluctuations (10, 11).

First, we study the appropriate working point for our network, i.e., we study the dependence of the collective dynamics as a function of the global coupling strength α and the delays through the velocity parameter v. In particular, because we are interested in cluster synchronization as a possible mechanism for generating the underlying anticorrelated subnetworks typical of the resting state, we first identified a division of the network in clusters using a modularity algorithm (22) (see Methods and SI for details). We found that the network can be subdivided in 2 communities (shown in Fig. 1). We note that these 2 communities are highly similar to the ones found in ref. 12.

To study the collective dynamics, we study the level of cluster synchronization in each community as a function of our free parameters. Two hundred forty seconds of the whole network dynamics were simulated, by employing an optimized Matlab routine (DDESD) based on Runge–Kutta's algorithm. Fig. 2 shows the level of synchronization in each of the 2 extracted communities [i.e., the figure plots the maximum of the Kuramoto indices (defined in Methods and SI) of both communities]. The figure shows that for a critical coupling α, there is a transition from a collective silent state (all nodes show low activation) to a synchronized global regime. Nevertheless, the synchronization is relatively low for most of the parameters combinations. However, there are 2 regions of parameter space that show elevated levels of synchronization that correspond to the increase of synchronization in either one of the community clusters: The left bump corresponds to one of the communities and the right bump to the other community. We fix our working point P (indicated Fig. 2 by a black asterisk) between the 2 synchronization bumps (α = 0.007 and v = 1.65 m/s). The reason is that we expect that in this region, we would find maximal cluster synchronization caused by fluctuations between the synchronization states of the 2 clusters.

Fig. 3A Left shows that this working point can reproduce the typical collective brain dynamics found at rest conditions. The black and red curves, respectively, plot the level of synchronization for each of the 2 communities as measured by the Kuramoto order parameter (see Methods and SI). The blue curve indicates the differences between the level of synchronization in the 2 clusters. Fig. 3A Right shows the power spectrum of the signal given by differences between the level of synchronization between both communities. The figure illustrates that at the chosen optimal working point, both slow 0.1-Hz oscillations of the synchronization signal and anticorrelation of the level of synchronization between the communities occur. However, each community does not show individually a 0.1-Hz modulation of their neural population activity, which underscores the relevance of neural synchronization as a mechanism for the emergence of the ultraslow fluctuations in the BOLD signal. In this figure, the level of noise is optimal (〈ν²〉 = 0.1), as we will see in the next section. In all these figures, we normalized the results to relative variations with respect to the mean (i.e., z = (z − 〈z〉)/〈z〉). Note that this normalization is done with respect to the mean of the particular time series of the community under consideration and not with respect to the global brain activity, which may cause artifactual anticorrelation. In other words, the anticorrelation patterns that we find are genuine and not a product of a normalization with respect to the global activity of the whole brain (23).

To study the relevance of the different intrinsic properties of the network we perform the same simulations but with a different level of noise (〈ν²〉 = 2) (Fig. 3B), eliminate the delays (Fig. 3C, note the different scaling of the y axis), and choose different working points (Fig. 3D shows just only 1 case given for α = 0.007 and v = 3.5 ms, but similar results were obtained for different working points). These results demonstrate that all these 3 factors are extremely relevant for obtaining the resting-state dynamics. In particular, the velocity parameters are probably constrained by the fact that the relevant emergent resting-state effects result from the equilibrated coordination between the fast local dynamics and the delays in the network. In our case, for the realistic 40-Hz range, we obtained the optimal synchronization level at the above selected working point P.

We also calculated the BOLD-signal using the Balloon–Windkessel hemodynamic model of Friston et al. (24), which specifies the coupling of perfusion to the BOLD signal, with a dynamical model of the transduction of neural activity into perfusion changes. Fig. 4A Left plots the BOLD signal calculated from the model at the optimal working point P. The figure shows that the model can reproduce both the slow 0.1-Hz oscillations and the anticorrelation of the BOLD signals of both communities. The black and red curves plot, respectively, the BOLD signal for each of the 2 single communities. The blue curve represents the difference between the BOLD level in the 2 clusters. Fig. 4A Right shows the power spectrum of the BOLD signal given by the differences between the level of BOLD signal in the 2 communities. In all these figures, we normalized the results to relative variations with respect to the mean (i.e., z = (z − 〈z〉/〈z〉) (again, note that the normalization is with respect to the mean value of the particular time series, i.e., community under consideration, and not with respect to the whole brain activity). Fig. 4B contrasts the relationship between the BOLD signal (black curves) and the level of synchronization (blue curves) on both communities. The curves show that a peak in the Kuramoto synchronization parameter computed from fast voltage–time data reliably precede peaks in the simulated BOLD response. The relationship is offset by a 1- to 3-s hemodynamic delay. Let us note that Honey et al. (12) have also detected a relationship between fluctuations in synchrony and BOLD response, although in the absence of time delays, which is crucial for the mechanism presented here. In conclusion, the 0.1-Hz slow oscillations resulting from alternancy in the level of synchronicity of the 2 communities is the origin of the observed 0.1-Hz slow oscillations of the BOLD signal.

The fact that synchronization predicts BOLD activity is not trivial. This is because the drive to the hemodynamic responses reflects mean population activity and not its synchronization. Our results, therefore, mean that there is a coupling between the degree of synchronization and neural activity that is manifest in elevated BOLD signals. This coupling has been studied in the context of evoked responses (25) and in terms of endogenous fluctuations (26). These analyses of simulated spike trends and local field potentials show that in nearly every domain of parameter space, mean activity and synchronization are tightly coupled, allowing us to conclude that indices of brain activity that are based purely on synaptic activity (e.g., functional magnetic resonance imaging) may also be sensitive to changes in synchronous coupling. Thus, our simulations explain why BOLD might be particularly sensitive to slow fluctuations in fast synchronized dynamics.

Role of Fluctuations: Stochastic Resonance.

To study the role and relevance of noise on the collective dynamics of the brain networks, we simulated systematically the behavior of the brain network for different levels of noise. We fixed all parameters according to the optimal working point P (α = 0.007 and v = 1.65 m/s) defined in the previous section and performed the simulations for 240 s. Fig. 5A plots the dependence of the maximum of the power spectrum peak of the signal given by differences between the level of synchronization between both communities (measured at the neuronal level as specified above) versus the noise level (variance of the stochastic fluctuations). This gives us a measure of the level of fluctuation that has a maximum effect on the emergence of global oscillations. In all plots, points (diamonds) correspond to numerical simulation results, whereas the lines correspond to a nonlinear least-squared fitting by using an α-function. As the figure shows, there is a stochastic resonance effect, i.e., there is a specified level of noise for which the optimum is reached. Lower or higher levels of noise attenuate the global 0.1-Hz oscillations. Fig. 5B plots the dependence of the location (in frequency domain) of the maximum in the power spectrum of the signal given by differences between the level of synchronization between both communities versus the noise level. This measure specifies the position of the maximum of the global oscillation. The figure shows again a stochastic resonance effect for the same level of noise. Furthermore, at this optimal level of noise, the maximum of the spectrum is given by 0.1-Hz global slow oscillations. Finally, Fig. 5C plots the level of correlation between the level of synchronization between both communities versus the noise level. Stochastic resonance at the same level of noise reveals a maximum of the anticorrelation between the 2 subnetworks, consistent with the experimental data. It is important to remark that not only the essential role of fluctuations, as documented by the presence of a stochastic resonance effect, but also the fact that the optimum level of fluctuations is given simultaneously for the emergence of 0.1-Hz global slow oscillations and the emergence of anticorrelated spatiotemporal patterns for both communities.

Connectivity Data on the Macaque Brain.

To study the intrinsic properties of the brain at rest, we performed all of the simulations and analyses using a connectivity matrix for 1 macaque hemisphere based on data from the neuroinformatics database CoCoMac (http://cocomac.org).

Network Dynamics.

The collective dynamics of a network of identical neuronal populations is determined only by the neuroanatomical connectivity matrix, if the dynamics of the single nodes is simple (e.g., Kuramoto oscillators) and if the transmission of information between different cortical nodes is instantaneous. Honey et al. (12) have shown that a much richer and more complex behavior (like the one evidenced during resting state) could emerge if the single-node dynamics are more complex (in particular, they used a more elaborate neuronal population dynamics that shows chaotic behavior). Another method to get a more complex collective network dynamics is by using a simple dynamics for each node but assuming realistic delay in the signal transmission between nodes in the network.

In this article, we concentrate on this last alternative. We assume realistic delays and take a simple realistic dynamics given by a Wilson–Cowan oscillator (35). The reason for this is that we would like to focus on the role of delays and, at the same time, consider how the typical global slow oscillation at rest could emerge from a network built up with simple fast oscillators (in the γ-band of 40 Hz). The equations describing the dynamics are given in the SI.

Cluster Synchronization.

Cluster synchronization is studied by defining a Kuramoto order parameter for each community in a sliding time window of 500 ms shifted by steps of 50 ms. We first shift the excitatory (x) and inhibitory (y) component of all nodes in a value μx and μy, respectively (i.e., x = x − μ^x and y = y − μ^y), such that the oscillations are centered around the origin. In each time window starting at time t_i and ending at time t_f, a measure of the degree of synchronization in each community M is given by,

where 〈〉 denotes average over time along the corresponding time window, and î is the imaginary unit. In all cases, we plot the normalized Kuramoto parameter normalized over all time windows along, i.e.

where now the average is taken over all time windows.