Abstract

Adaptive behaviour, cognition and emotion are the result of a bewildering variety of brain spatio-temporal activity patterns. An important problem in neuroscience is to understand the mechanism by which the human brain's 100 billion neurons and 100 trillion synapses manage to produce this large repertoire of cortical configurations in a flexible manner. In addition, it is recognized that temporal correlations across such configurations cannot be arbitrary, but they need to meet two conflicting demands: while diverse cortical areas should remain functionally segregated from each other, they must still perform as a collective, i.e. they are functionally integrated. Here, we investigate these large-scale dynamical properties by inspecting the character of the spatio-temporal correlations of brain resting-state activity. In physical systems, these correlations in space and time are captured by measuring the correlation coefficient between a signal recorded at two different points in space at two different times. We show that this two-point correlation function extracted from resting-state functional magnetic resonance imaging data exhibits self-similarity in space and time. In space, self-similarity is revealed by considering three successive spatial coarse-graining steps while in time it is revealed by the 1/f frequency behaviour of the power spectrum. The uncovered dynamical self-similarity implies that the brain is spontaneously at a continuously changing (in space and time) intermediate state between two extremes, one of excessive cortical integration and the other of complete segregation. This dynamical property may be seen as an important marker of brain well-being in both health and disease.

1. Introduction

It is increasingly evident that brain regions are continuously interacting even when the brain is ‘at rest’ and, more importantly, that the functional networks uncovered from resting data closely match those derived from a wide variety of different activation conditions [1,2]. Starting with the uncovering of coherent fluctuations of functional magnetic resonance imaging (fMRI) in time series of motor cortex [3], many other findings have validated the notion of correlated networks as a dynamical substrate of the resting brain. It has been established that these networks, which can be separated on the basis of their temporal features [4–6], are located at consistent locations across subjects and are equally detectable even during sleep [7] and anaesthesia [8].

These exciting findings provide a novel window to observe the brain at work and, at the same time, highlight our limited understanding of the functional organization of the brain at large scales [9], compared with the, often precise, knowledge we have of the small (neural circuit level) scale. In particular, little is known about how the cortex is able to solve the conflicting dynamical demands imposed by the functional segregation of local areas differing in their anatomy and physiology, on one hand, and, on the other, their global integration shown during perception and behaviour. This riddle is clearly pointed out by Tononi et al. [10]: ‘traditionally, localizationist and holist views of brain function have exclusively emphasized evidence for either functional segregation or for functional integration among components of the nervous system. Neither of these views alone adequately accounts for the multiple levels at which interactions occur during brain activity’ (p. 5037). In connection with this dilemma, it has been suggested that the brain's conflicting demands are a generic property of many collectives, regardless of being composed by neurons, genes, individuals, etc. [11,12]. This implies that the search for the mechanism behind this dilemma must be guided by the general properties of the system, rather than by the details of the neurobiology. In that regard, the study of large-scale collective properties has a long tradition in statistical physics, allowing the identification of different dynamical regimes through the study of ubiquitous correlation properties. In this work, we are guided along that direction, and the paper is dedicated to reporting the spatio-temporal correlation properties of fMRI resting-state data that are found to exhibit robust self-similarity signatures in space and time. This finding has profound significance, because it demonstrates a continuous range of correlations between the local cortical circuits up to the entire brain. Thus, this is the first direct empirical demonstration of the brain resting-state activity exhibiting simultaneously functional segregation and integration.

A variety of experiments have already provided indications of self-similarity in spatial and temporal scales of the brain dynamics. For example, at small scale, avalanches of neuronal activity in rat [13,14] and monkey [15] cortex are known to be scale-free, a finding that has been modelled by de Arcangelis et al. [16], Levina et al. [17] and Buice & Cowan [18] among others. At a larger scale, it was reported that brain fMRI networks are characterized by scale-invariant degree distributions, as described by Eguíluz et al. [19] and later replicated and extended by van den Heuvel et al. [20]. These networks are indistinguishable from those extracted from model systems at criticality [12,21]. Other findings include the observation of 1/f power spectra from simultaneously recorded magnetoencephalography and electroencephalography signals [22], from fMRI signals [23–25] as well as from cognitive responses [26], all indicative of long-range correlations in brain dynamics processes. It is these myriads of observations of power-law behaviour and long-range correlations that have led to the conjecture that the human brain as a whole behaves as a system at criticality [27–29] and that the framework of self-organized criticality [29–31] may be relevant to understand large-scale brain dynamics. In that sense, critical dynamics endows the system with a high susceptibility and a broad repertoire of responses, natural requirements for healthy brain function.

Here, we analyse directly the correlation function of the voxel signals in space and time. A voxel is assigned a set of three weights describing the density of grey matter, white matter and cerebrospinal fluid (CSF) within the voxel. Up to now, the standard analysis of fMRI functional correlations only included voxels with a density of grey matter above a certain threshold value. In the present analysis, we do not discard any voxels, nor do we apply a threshold on the correlation coefficients. Hence, the present analysis corresponds to the standard analysis tool applied in physics and materials science when structural properties are investigated.

To investigate whether the system is self-similar in space, we perform a coarse-graining analysis directly on the voxel data and then extract the associated correlation function. We will show that the coarse-grained correlation function exhibits a power-law decay of correlations at length scales between the microscopic scale of a voxel and the macroscopic scale of the size of the brain. This self-similarity is reminiscent of the multi-scale modular organization [32,33] observed in fMRI data [34], but is also a robust finding in systems near the critical point [29–31].

2. Results

2.1. Spatial self-similarity

Starting from the computation of the correlation matrix, the average correlation function for each separation r can be estimated. To study the degree of self-similarity present in the correlation function, we perform a coarse-graining procedure and check for self-similarity via equation (4.6) as detailed in §4. The result of this analysis for a single subject is shown in figure 1. The electronic supplementary material contains the results for six additional subjects. In all cases, the data were binned into integer distance values to decrease the statistical fluctuations (especially at long distances) without having to average across several subjects. Two different regimes are observed in the correlation function. At short distances, we see a collapse of the different curves so equation (4.6) holds, indicating that the system is self-similar. At longer distances, instead of the usual decay in the correlation function that is observed in physical models, there is a plateau, probably caused by the brain's bilateral symmetry, as seen previously across cortical bilaterally homologous regions by Salvador et al. [35,36], and then an increase in the correlations, probably caused by surface effects.

After one coarse-graining step, the length of a side of a block-voxel is twice the size of a voxel, but then we need to renormalize the Euclidean distance in the coarse-grained system by a factor 2 to recover distance in voxels. If the system under consideration were infinite, then 〈C(r)〉 would be defined for all r irrespective of the level of description and they would be identical (equation (4.6)). However, as the system we consider is finite, there exists a maximum distance rm in the original system. After one coarse-graining step, this maximal distance is equal to rm/2, and after n coarse-graining steps this maximum distance equals rm/2n, measured in units of voxels at the nth coarse-graining level. If equation (4.6) holds, then it means that the system behaves in the same manner at all scales, up to the maximum distance possible. Furthermore, it allows for a data collapse of the correlation functions by first rescaling the distance r in voxel units to Euclidean distances, r ↦ r2na, where a is the voxel spacing in Euclidean distance (say millimetres) in the original data and then renormalizing the correlation functions by the factor 2−αn, where α is a subject-dependent parameter. It is easy to check that if the correlation function depends on distance through a power law with exponent β, that is, 〈C(r)〉 ∝ r−β, then the exponent α in the renormalizing factor of the correlation function must be equal to β. The relation α = β is fulfilled for this regime.

Scaling is only expected to apply in an intermediate regime where the sizes of the coarse-grained voxels are big compared with the smallest length scale cut-off in the data, that is, the voxel scale set by the scanner, and the largest distance covered by the data, that is, the size of the brain. This is indeed observed in figure 1, where the data collapse is best for the three curves corresponding to the bigger block voxels.

To verify that the shape of the average two-point correlation function is not due to artefacts in the data, we randomized the positions of the voxels or the time series (figure 2). By randomizing the positions, we expect and find a constant profile for the correlations. In an infinite system or in a system where the correlation length is much smaller than the system's size, this constant would be equal to zero; however, in the system we study, strong long-range correlations are present, hence the constant is non-zero. By randomizing the time series of each voxel independently, we expect to destroy all correlation among voxels and thus expect and find close to zero average correlation for all distances. Furthermore, we also scanned a phantom to test our finding against artefacts coming from the scanner. The phantom data show an exponential decay and no long-range correlation, except a peak at the longest distance that is also observed in brain data, confirming that it is a surface effect.

The average correlation function 〈C〉 versus the voxel distance after one step of coarse graining (64 × 64 × 16) for fMRI data (solid line), spatially randomized voxels (dashed-dotted line) and temporal-randomized voxels (dotted line). Also shown is the average correlation function for phantom data (dashed line). (a) Linear–linear, (b) linear–log and (c) log–log axis, respectively. The correlation function for the fMRI data decays as a power law, while the correlation function for the phantom data decays exponentially fast. The correlation function for the spatially randomized data is constant while it is zero for the temporally randomized data.

When analysing fMRI data, the focus is, unlike our approach, on grey matter voxels. Owing to the large size of the voxels, the discrimination between grey and white matter voxels is based on a high-resolution anatomical scan that gives the percentage of grey matter, white matter and CSF for each voxels. A voxel is said to belong to grey/white matter if its content in the latter exceeds a threshold that has to be fixed; the others are dismissed. This procedure has the drawback that information is lost in the process because the grey matter signal in dismissed voxels is lost.

Pure grey and white matter voxels are rare, but, from a high-resolution scan, we know the percentage of grey and white matter of each voxel. Hence, we can weight the correlations between two voxels according to their relative content of grey and white matter with the following normalization: 1 = wg + ww + wCSF. We have performed the same analysis as we did above. Using all voxels weighted with their contents of grey and white matter, respectively, we calculated the average correlation functions
2.1

The correlation of grey matter is shown in figure 3 and the correlation of white matter is shown in figure 4. The bump at the end of the distribution survives only in the grey matter, which is to be expected, since it is probably due to a surface effect. White matter correlations die off smoothly before the surface of the brain is reached, which was expected since white matter is supposed to be enclosed in grey matter. The fact that the exponents α ≠ β in both figures 3 and 4 is related to the deviation from power-law behaviour of the correlation functions of the grey and white matter. The scaling collapse is better for the correlation function when all voxels are included (figure 1) than for the grey matter and white matter correlation functions displayed in figures 3 and 4. It is remarkable that, only when the entire fMRI BOLD signal (all voxels without thresholding) is included in the correlation matrix, we obtain a good quality data collapse consistent with self-similarity. In general, we associate the essential dynamics and processing of the brain with grey matter, so our findings indicate that the separation of the fMRI BOLD signal using the weights wg and ww is more subtle than expected.

In the electronic supplementary material, we display the average correlation functions 〈C〉, 〈Cgg〉 and 〈Cww〉 for an additional six subjects. These correlation functions are consistent with those displayed in figures 1, 3 and 4. Table 1 summarizes the associated exponents for all seven subjects.

The associated exponents α and β for the average correlation functions 〈C〉, 〈Cgg〉 and 〈Cww〉, respectively, for seven young healthy subjects. The uncertainty of the last digits is given by the figures in the brackets. The corresponding graphs are displayed in figures 1, 3 and 4 and electronic supplementary material, figures S1–S7. The data can be requested on a CD by contacting H.J.J.

2.3. Temporal self-similarity

To complement the spatial correlation analysis, we present the average power spectra of the fMRI time series. Temporal auto-correlations and power spectra have been studied by use of fMRI for a while. A comprehensive review is given in Woolrich et al. [37] and the use of these as a diagnostic tool was discussed in Maxim et al. [24]. To obtain the spectra, we computed the spectrum of each voxel individually and normalized its integral to 1. This corresponds to analysing the Fourier transform of the auto-correlation function. We note that the term ‘power spectrum’ is often used to denote the Fourier transform of the auto-covariance function. However, it is better to operate with the normalized correlation function to avoid signals with more power to totally dominate the average power spectra. Although the frequency range we can explore is rather small, the power spectrum has a reasonable 1/f dependence, thus showing criticality in time as well as in space (figure 5). To check against potential artefacts of the measurement, we calculated the power spectra of phantom data and surrogate data, constructed by randomizing the brain fMRI time series. In both cases, we obtain a flat spectrum, as shown in figure 5, rejecting the possibility that the spectral features arise from scanning artefacts. The peak at 0.03 Hz superimposed on the 1/f slope corresponds to the low-frequency fluctuation dynamics already described for the so-called brain resting-state networks [6].

Average power spectrum for the fMRI time series (solid line), the phantom time series (dashed line) and the randomized time series (dotted line). Solid straight line is a guide to the eye for a 1/f spectrum.

3. Summary and discussion

These results are the first direct demonstration of spatial and temporal self-similarity of the brain resting-state dynamics. Although several reports have already hinted about the possibility that at large scale the brain may be operating in such a state, previous measures have been indirect and/or model dependent. By contrast, in the present work, scale invariance in space and time is determined from the entire dataset without discarding any voxels and using the same renormalization techniques championed in the study of critical phenomena in physical systems and models. In addition, to retain all the information and focus on the global state of the brain, we do not impose a threshold on the correlation coefficients. Thresholding would imply throwing away information to highlight patterns of correlated activity.

However, it is interesting to relate our finding of scale invariance with the spatial and temporal localized functional activity extracted from fMRI data by independent component analysis (ICA) methods. By construction, ICA emphasize the identification of strongly independent components and ‘filter’ out less important contributions to the total BOLD signal. Recent studies have investigated the growing complexity of the spatial structures extracted by ICA as the total numbers of components are allowed to increase [2,38]. The additional structure obtained this way is indeed consistent with the self-similar character we have identified in the correlation function of the BOLD signal.

Consistent with the present results is a recent report of altered BOLD correlations in certain chronic brain diseases. Criticality implies magnetization zero; in other words, fluctuations are such that the system is broken in opposite domains of approximately equal size. This is precisely a striking feature demonstrable in the BOLD spatial correlations by computing the ratio between the brain area covered by positively correlated regions and those by negatively correlated sites. This ratio, analogous to the magnetization parameter above, was found to be about unity for healthy subjects but orders of magnitude larger in patients with chronic low back pain [39]. These findings, together with the present results, suggest that scale invariance could be an important objective marker of brain well-being in both health and disease.

We would expect similar findings for subjects performing a particular task paralleling extensive work on resting-state networks, which have shown that the fMRI spatio-temporal patterns seen during rest and execution of a task do not differ [2] . This is due to the fact that external inputs only produce minute fluctuations in brain activity, thus the large-scale dynamics remain essentially unchanged.

From a dynamical systems perspective, the uncovered self-similarity implies that the brain dynamics is permanently at an intermediate state between two extremes, one that is strongly correlated across large distances, producing transient highly integrated cortex states, and the other in which only nearby clusters are acting in sync. This scenario, of long-range correlations in space and time, is only conceivable in dynamical systems at criticality and could be the manner in which the cortex can manage to produce an arbitrarily large repertoire of interaction patterns among arbitrarily distant cortical sites.

4. Material and methods

4.1. Image acquisition

A 3T Philips MRI scanner was used to acquire T2*-weighted echo-planer images (EPI) in 128 × 128 × 31 voxels of dimension 2.5 × 2.5 × 5.0 mm3 using a repetition time of 2000 ms, echo time of 60 ms and a flip angle 90°. An eight-channel array coil and SENSE factor 2 were used as well as second-order shims. Finally, a three-dimensional mask was used to identify the content of white matter, grey matter and CSF.

4.2. Image pre-processing

The brain data consist of fMRI datasets of seven young healthy adults in the resting state. The data can be requested on a CD by contacting H.J.J. The subjects were instructed to lie still in the scanner with their eyes closed, avoiding falling asleep. A total of 305 functional volumes of each subject were acquired from each session, the first five of which were discarded in order to remove the effect of T1 equilibration. Image pre-processing was performed using the University of Oxford's FMRIB software library (FSL) involving realigning to account for movement using FMRIB's Linear Image Registration Tool (MCFLIRT) [40], high-pass temporal filtering using FEAT to remove low-frequency artefacts and a high-pass filter cut-off preset to 100 s.

4.3. Coarse graining of the signal

Coarse graining is a well-established technique to describe the system's changes at different levels of spatial or temporal observation. This allows one to investigate to what extent the spatial or temporal structures of a given phenomenon exhibit scale invariance (e.g. [41]). A self-similar or scale-invariance object, for instance a cauliflower, looks like itself at all scales. Each little floret of the cauliflower looks like a miniature version of the entire cauliflower. Hence, a short length-scale study of the cauliflower in which one probes the substructure of an individual floret will provide the same information as the study at a large length scale in which one observes the cauliflower as a composition of florets.

In general, the approach to probe the behaviour at different length scales consists in replacing the signal at a given point by a spatial average over a region of a given spatial extension about this point (figure 6). In the case of an fMRI signal, the procedure consists in aggregating 2 × 2 × 2 = 8 adjacent voxels to create a block-voxel ℬ′ whose BOLD intensity at time t, vℬ′(t), is the average intensity of the BOLD intensities at time t, vℬ(t), in its constituting voxels B ∈ ℬ′, that is,
4.1
see figure 6 for an illustration.

(a) Example of coarse graining in two dimensions where there are four boxes ℬ within a block-box ℬ′. (b) The four dashed-coloured signals from the four original boxes ℬ are averaged to produce the solid-black coarse-grained signal of ℬ′.

Computing equation (4.1) for all times t = 1,2, … , 300 gives the time series for the BOLD signal in the block-voxels {vℬ′(t)}. We now repeat the coarse-graining procedure on the time series for the BOLD signal in the block-voxels {vℬ′(t)} to produce a time series for the BOLD signal in the block-voxels {vℬ″(t)} and so on. For the spatial resolution of our fMRI data, the coarse-graining procedure applied three times yields voxels ℬ(128 × 128 × 31)↦ℬ′(64 × 64 × 16)↦ℬ″(32 × 32 × 8)↦ℬ″′(16 × 16 × 4) with associated signals vℬ(t)↦vℬ′(t)↦vℬ″(t)↦vℬ″′(t). In the first coarse-graining step, the ultimate plane of voxels at the 31st slice is coarse grained only along the x–y plane. Averages have been calculated over active voxels only to prevent the inclusion of non-active voxels at the boundaries.

Next in our analysis, we consider the time series of each (block) voxel as a realization of a random variable and compute the time correlation matrix at the four levels of coarse graining. This matrix is obtained by computing the correlation between all pairs of (block) voxels ℬi at position ri and ℬj at position rj,
4.2
where is the temporal average value of (block) voxel's ℬi BOLD signal, and is the standard deviation of the BOLD signal of (block) voxel ℬi and T = 300.

4.4. Scaling

From the correlation matrix, the correlation function at equal time Cℬ (vℬi, vℬj) is calculated as a function of the position in space of the voxels. We then compute the correlation function 〈Cℬ(r)〉 averaged over all voxels with r = |ri − rj| measured in units of (block) voxels. This is done for each level of coarse graining to obtain 〈Cℬ (r)〉, 〈Cℬ′(r)〉, 〈Cℬ″(r)〉 and 〈Cℬ″′(r)〉 and the indices ℬ, ℬ′, ℬ″ and ℬ″′ indicate the level of coarse graining.

Let us discuss what to expect concerning the relation of the correlation function calculated at one coarse-graining level, say 〈Cℬ(r)〉, compared with the one at the next level, 〈Cℬ′(r)〉. To clarify the effect of coarse graining, we relate the correlation function at this level of coarse graining to the covariance function at the previous level. For a homogeneous system, the covariance function is given by
4.3
where 〈·〉 denotes the spatial and temporal average. We consider, as in equation (4.1), the coarse graining over a block-voxel ℬ′ centred at r,
4.4
where the block-voxel ℬ′ contains |ℬ| voxels. We then have
4.5

For a scale-invariant (self-similar) system of infinite size, 〈Cℬ′(r)〉 will remain invariant if distance r is measured in units of coarse-graining boxes, that is,
4.6

Note that r is the distance expressed in the voxel unit corresponding to the given level of coarse graining.

The simplest case of scale invariance corresponds to a power-law behaviour of 〈C(r)〉. Figure 7 illustrates the different behaviour under rescaling of two three-dimensional artificial systems, where we imposed the correlation to have either a power law or an exponential behaviour,
4.7a4.7b

The average correlation function 〈C〉 versus the block distance for three coarse-graining levels of an hypothetical object with a power-law form of the correlation function. The graphs collapse in voxel distance. Inset: results from a hypothetical object having an exponential form of the correlation function where the curves do not collapse in block distance.

Note that only those correlations with a power-law form collapse in a block unit.

Acknowledgements

We thank Prof. Joseph Hajnal for providing the phantom data. It is a great pleasure to acknowledge very helpful discussions with Dr Christian F. Beckmann. The project was supported by the EPSRC under grant EP/E049451/1. D.R.C. is supported by NINDS grant NS58661(USA).