Despite the very large data reduction, we found that the information content of the few remaining points is very high. As a proof of principle, we first used the point process to calculate the spatial location of six well known RSN maps. These maps describe the major independent components of brain spontaneous activity, and as such they can be used as a relevant benchmark. We used the point process to obtain the RSN maps and compare them with maps computed from the full BOLD signal using a well established method (probabilistic independent component analysis – PICA; Beckmann et al., 2005). This is done by calculating in six RSNs the rate of points co-occurrence (up to 2 time units later in this case) between representative sites (“seeds”) and all other brain voxels and presented as maps in Figure 2A–C (see Materials and Methods for a detailed explanation of the computation). The seeds locations were selected according with previous work (see coordinates in Table 1 of Materials and Methods Section). The similarities between our conditional rate maps and the respective PICA maps (rightmost three columns and left column of Figure 2A respectively) is already obvious to the naked eye and confirmed by the correlation plotted in Figure 2C. The calculation shows that despite using less than 6% of the raw fMRI information, about 5 points (on average) are enough to obtain RSN maps that are highly correlated (95% confidence) with those obtained using PICA of the full BOLD signals. Similar good performance can be demonstrated in tracking physio-pathological changes of brain activation. This is presented in Figure 2D which shows the statistical differences in functional connectivity between a group of chronic back pain (CBP) patients and healthy controls already reported in Tagliazucchi et al., 2010a; comparison with seed correlation based in the DMN, increased correlation with bilateral insula in CBP). Finally, the data analysis from a finger tapping task (Tagliazucchi et al., 2010b) demonstrates also the merits of the current approach when compared with a seed correlation based in primary motor cortex contralateral to the tapping hand (Figure 2E).

The results in the previous section show that the point process can efficiently compress the information needed to reproduce the underlying brain activity in a way comparable with conventional methods such as seed correlation and independent component analysis. Importantly, while the former methods represent averages over the entire data sets, the point process, by construction, compresses, and preserves the temporal information. This potential advantage, unique of the current approach, may provide additional clues on brain dynamics. This is explored here by compiling the statistics and dynamics of clusters of points both in space and time. Clusters are groups of contiguous voxels with signal above the threshold at a given time, identified by a scanning algorithm in each fMRI volume (see Materials and Methods for details). Figure 3A shows examples of clusters (in this case non-consecutive in time) depicted with different colors. Typically (Figure 3B top) the number of clusters at any given time varies only an order of magnitude around the mean (~50). In contrast, the size of the largest active cluster fluctuates widely, spanning more than four orders of magnitude.

The analysis reveals four novel dynamical aspects of the cluster variability which hardly could have been uncovered with previous methods. (1) At any given time, the number of clusters and the total activity (i.e., the number of active voxels) follows a non-linear relation resembling that of percolation (Stauffer and Aharony, 1992). At a critical level of global activity (~2500 voxels, dashed horizontal line in Figure 3B, vertical in Figure 3C) the number of clusters reaches a maximum (~100–150), together with its variability. (2) The correlation between the number of active sites (an index of total activity) and the number of clusters reverses above a critical level of activity, a feature already described in other complex systems in which some increasing density competes with limited capacity (Stauffer and Aharony, 1992; Bak, 1996). (3) The rate at which the very large clusters (i.e., those above the dashed line in 3B) occurs (~ one every 30–50 s) corresponds to the low frequency range at which RSN are typically detected using PICA (Beckmann and Smith, 2004; Beckmann et al., 2005). (4) The distribution of cluster sizes (Figure 3D) reveals a scale-free distribution (whose cut-off depends on the activity level, see Figure 3F).

These four features remind of other complex systems undergoing an order-disorder phase transition (Bak, 1996; Jensen, 1998; Tsang and Tsang, 1999; Chialvo, 2010; Tagliazucchi and Chialvo, 2011) thus suggesting further exploration. Following standard techniques in statistical physics, two parameters were defined and computed from the same data plotted in Figure 3C. To represent the degree of order (i.e., the order parameter), the size of the largest cluster (normalized by the number of active sites) in the entire brain was computed and plotted as a function of the number of active points (i.e., the control parameter). This was done for all time steps and plotted in Figure 3E (small circles). We avoided the use of the branching ratio (Beggs and Plenz, 2003) as a control parameter because its estimation from the data is less than straightforward. It cannot be computed for each fMRI volume as required here and only converges to a stable quantity for relatively long time series. In addition, it requires the ad hoc definition of the number of bins and a suitable bin-width for its computation (Beggs and Plenz, 2003), therefore making its use cumbersome for the spatiotemporal resolution of the present study. On the other hand, the parameter used here (i.e., global level of activity) is computed in a straightforward manner, converges relatively fast, requires no fine tuning of parameters and has clear analogies to control parameters of well studied models of order-disorder transitions (the clearest example being percolation; Stauffer and Aharony, 1992).

Several key features are worth to mention, all highly suggestive of a phase transition: First, there is sharp increase in the average order parameter (empty circles), accompanied by an increase of its variability (empty squares). Second, the transition coincides with the peak in the function discussed in Figure 3C, which accounts for the number (not the size) of the clusters. Finally, a calculation of the relative frequency of the number of active sites was performed (i.e., residence times, filled circles) showing that the brain spends, on the average, more time near the transition than in the highly ordered or the highly disordered states. This is a remarkable support for earlier conjectures suggesting that the brain at large-scale works at criticality (Bak, 1996; Chialvo, 2010; Expert et al., 2011; Tagliazucchi and Chialvo, 2011).

The identification of a phase transition in the resting brain suggested additional work to characterize its properties, including a quantification of the dynamical properties of cluster spatial evolution. As shown in the example of Figure 4A an activated cluster can appear, grow to achieve a maximum size and then disappear (or translate or divide into sub-clusters). The present approach allows the study of two properties of the process. For each cluster we first measured a static space filling property, the average fractal dimension D. This is shown in Figure 4B which illustrates that D ~ 2.15 ± 0.02. While the fractal dimension D departs from this value for the highly ordered and disordered regimes (Figures 4D,E), the residence time distribution computed in Figure 3E indicates that most of the time the level of activity is around the critical value, thus on average D ~ 2. Second, we looked at the dynamics of the cluster propagation, which was found to happen in bursts. The statistics in Figure 4C shows that avalanches could last up to 30 s. with sizes up to 10³ and have no preferred scale, a behavior very similar to that of neuronal avalanches described previously in smaller scales (Beggs and Plenz, 2003; Petermann et al., 2009; Chialvo, 2010).

At first sight, the continuous nature of the fMRI BOLD signal, imposed by the nature of the neurovascular coupling itself (Friston et al., 1995, 1998), might have hindered the introduction of point process methods for its analysis. The situation is analogous to that of continuous rhythmic activity arising in scalp EEG due to predominant frequencies in the spiking activity and sub-threshold oscillations which underlie the generation of discrete action potentials (Traub et al., 1989; Steriade et al., 1993; Contreras and Steriade, 1995). We have shown that the application of HRF deconvolution gives a way to invert the process and find the train of impulse-like signals (of whatever origin) which closely resembles the point process. The fact that the majority of the points coincide with the peaks of the BOLD HRF-de-convoluted signal (i.e., Figure 1C) reinforces the view that upward going BOLD signals are non-linear events where the crossing times preserve the most relevant information. This is in line with recent findings of all-or-none “coherence potentials” macroscopically propagating in the monkey cortex, as observed in local field potential recordings (Thiagarajan et al., 2010). Therefore it seems reasonable the conjecture that, at this level of coarse graining, we are dealing with all-or-none intermittent avalanching events which involve short and long range cortical co-activations.

Another remark needs to be made concerning the HRF: while it is true that extensive work established the fundamental details of the brain’s BOLD HRF responding to a well defined (single or repetitive) stimulus, less is known about the BOLD response under the non-stationary conditions of resting state, in other words, the nature of the resting state HRF remains unknown. A theoretical formalism for the neuro-BOLD coupling at rest, as far as we know, has not been attempted but probably deserves to be considered in the future. Such studies should clarify up to what extent the HRF function obtained from stimuli spaced by relatively long intervals can predict the temporal evolution of the BOLD signal measured during resting state.

The reduction of the fMRI BOLD signal to discrete events not only allows for the identification of well-described resting state networks as shown in Figure 2, but also reveals that large-scale brain activity organizes in avalanches of activity with power law size distributions. The point process approach allowed for the first time to identify explicitly the order and control parameters and to define the state of the resting brain as a fluctuation around a phase transition. The analysis shows not only that activity spreads as scale-free avalanches resembling those seen in smaller scales (Beggs and Plenz, 2003) but – and importantly – that the brain spends most of the time at a level of activity which corresponds to the critical point. These new findings add to the previous observations that the correlation function of fMRI BOLD signals exhibits fractal properties (Expert et al., 2011) and that the correlation length of the activity measured with fMRI diverges as predicted by the theory of phase transitions (Fraiman and Chialvo, 2010), supporting the hypothesis that brain dynamics operates at a critical point of a second order phase transition (Bak, 1996; Chialvo, 2010; Tagliazucchi and Chialvo, 2011).

In connection with previous experiments, one must emphasize that for a non-equilibrium system in a critical state avalanches are observed at a wide range of temporal and spatial scales. Observations at smaller scales (Beggs and Plenz, 2003) show a clear cut-off of the avalanche distribution at the size of the electrode array used for the recordings (signaling that the experimental technique is unable to sample larger events) as well as a distortion of the distribution caused by sub-sampling effects (Priesemann et al., 2009). Due to very good spatial resolution and whole brain coverage fMRI allows to overcome these issues in the macroscopic domain (≈1 mm). The observation of identically distributed avalanches at this level is direct evidence that the brain spatiotemporal dynamics is scale-free, as expected for a critical system. The present work also suggests the study of intermediate scales accessible by means of other experimental techniques to give further support for or against this hypothesis.

The observation that large-scale brain dynamics can be traced as discrete scale-free avalanches of activity raises the question of the physiologically relevance encoded in the timing of these large-scale events, already suggested by observations at smaller scales and computational models (Kinouchi and Copelli, 2006; Shew et al., 2009, 2011; de Arcangelis and Herrmann, 2010). For instance, although relatively rare, avalanches in the tail of the power law distribution emerge from a local origin and propagate as far as the length of the entire cortex, suggesting a role in the binding processes of far apart cortical regions. It would be interesting to investigate whether total or partial disruption of these large events, as well as alterations in the balance between activation and segregation into clusters are correlated with pathological conditions and with the level of awareness of the subject. Additionally, the non-linear relation between activated cortical tissue and number of clusters exhibits an optimal point, in which the level of brain activity is segregated into the maximum number of spatially isolated activations. We can hypothesize that this result is relevant to the solution of the integration/segregation dilemma long advocated by Tononi et al., 1994; Sporns, 2010) as the fundamental conundrum that the healthy cortex needs to be executing at any given time. If our hypothesis is true, we can predict, together with integration/segregation theories of consciousness (Tononi et al., 1994; Tononi, 2004), that a displacement of the optimal point should be observed for brain states of diminished conscious content such as deep sleep, anesthesia, or coma (Lee et al., 2009).

A special place in the discussion should be dedicated to analyze the similarities and differences between the definition of the order parameter and phase transition used here and the concept of synchrony widely used in previous studies. To place this point in context, it is appropriate to recall the earlier studies a decade ago, by Varela and colleagues (Rodriguez et al., 1999) investigating the brain electrical activity of subjects viewing ambiguous visual stimuli (perceived either as faces or as meaningless shapes). They were able to show for the first time that “only face perception induces a long-distance pattern of synchronization, corresponding to the moment of perception itself and to the ensuing motor response. A period of strong desynchronization marks the transition between the moment of perception and the motor response” (Rodriguez et al., 1999). These results lead to the authors to suggest that “this desynchronization reflects a process of active uncoupling of the underlying neural ensembles that is necessary to proceed from one cognitive state to another” (Rodriguez et al., 1999). A number of papers followed Varela et al. idea (Rodriguez et al., 1999) that synchrony is physiologically relevant. At the fMRI level the timing and length of these epochs of synchrony were used recently to infer the presence of criticality (Kitzbichler et al., 2009). The present results indicate that while order (as defined here) implies always synchrony, the reverse is not always true, since space is not manifest in the definition of synchrony and then one can have a very synchronic but (spatially) disordered pattern of brain activity. Since our results indicate that the brain at rest spends most of its time in a mix of order and disorder, it would be very interesting in future studies to relate Varela’s synchrony-asynchrony concept with the current ideas of proximity to a order-disorder transition.

Attempts to construct biologically realistic equilibrium models of brain networks require as a main ingredient the introduction of (sometimes finely tuned) noise (Deco et al., 2009; Rolls and Deco, 2010). In this type of models, without the external noise the dynamics are stuck in a stable equilibrium state, thus noise must be introduced ad hoc to allow sufficient variability in the dynamical behavior of the system. One should be very careful, however, not to over-emphasize the biological relevance of a construct needed to overcome the shortcomings of a restricted class of models. Statistical physics results tell us that dynamical fluctuations around stationary states are small except near critical points (Prigogine, 1962).On the contrary, a non-equilibrium system undergoing criticality does not need the introduction of noise: variability is self-generated by the collective dynamics which spontaneously fluctuate near the critical point (for further discussion, see (Tagliazucchi and Chialvo, 2011)). Coincidently, the present results show that the spatiotemporal organization of the resting brain dynamics achieves maximum variability (i.e., Figures 3C,E) at a particular level of activation, and the analysis of the order and control parameters reveals that the origin of such variability can, in fact, be traced to a phase transition. Furthermore, the level of activity spends the largest amount of time around such transition. Then, these results point out that a different class of models is needed: one that emphasizes non-equilibrium self-generated variability over ad hoc introduced noise of uncertain origin.

The present results gathered in a large-scale domain can be also analyzed at the light of earlier observations of transient states at faster time scales, in which the scalp EEG is reduced to a certain number of stereotypical topographical maps (i.e., EEG microstates) (Koenig et al., 2002) and non-stationarities which define discrete segments of electrical activity are observed (Kaplan et al., 2005). Both descriptions of electrical scalp activity have also been shown to exhibit properties consistent with critical dynamics (Allegrini et al., 2010; Van De Ville et al., 2010). Further multimodal imaging studies could link these observations together in the context of discrete avalanches of neural activity propagating through the cortex and determine their functional relevance for health and disease (Greicius et al., 2004). Also, future work on the analysis of large-scale spontaneous fMRI signals as a train of activations should take advantage of the fact that the temporal information is not completely discarded (as in a straightforward correlation analysis) but kept in the timing of the events. The point process extracted from the BOLD signal can thus provide valuable information on the transient co-activations (or co-participation in an avalanche) of different brain regions. This measure can then be of value if correlated with the aforementioned EEG measures of instantaneous synchronization, as well with spontaneous index of perception or task performance.

The most obvious limitations of the point process approach stem from the spatiotemporal resolution of the fMRI recordings (i.e., TR and voxels dimensions) as well as the time constant of the BOLD HRF. Because of these limitations it is in principle impossible to distinguish two points in the process which are spaced by less than a characteristic time, as well as to detect very small clusters or avalanches whose size is smaller than the voxel dimensions. Therefore, it is impossible to guarantee that all points and clusters are included in the statistical analysis. However, since those which may be left out have (by definition) the smallest contributions, results are unlikely to be affected by this limitation. Another possible drawback of the method is that, while it yields more information than other methods such as linear correlation, there is a free parameter (threshold) to select. Nevertheless, in the Materials and Methods section we show that results are robust against different threshold choices.

Overall, the results show that the location and timing of the largest BOLD fluctuations define a spatial point process containing substantial information of the underlying brain dynamics. Despite the very large data reduction (>94%), the approach was validated by the favorable comparison of the conditional rate maps of avalanching activity with those constructed with the full fMRI BOLD signals using PICA as well by comparison with two distinct pathophysiological conditions (resting state in CBP patients and a finger tapping task). In addition to uncover new dynamical properties for the activated clusters, the method exposed scale-invariant features conjectured in the past (Chialvo, 2010) which are identical to those seen at smaller scales (Beggs and Plenz, 2003; Petermann et al., 2009; Chialvo, 2010). For the first time, the order and control parameters have been derived from human fMRI data allowing the identification of a phase transition and the demonstration that the resting brain spends most of the time near criticality. Beyond its potential value for fMRI signal processing, the ability of the present approach to capture relevant spatiotemporal brain dynamics underscoring non-linear aspects of the BOLD signal deserves further exploration.

Data was obtained, after informed consent, from ten right-handed healthy volunteers (9 female, 1 male; mean age = 49, SD = 12), during 10 min, requested to keep their eyes closed and to avoid falling asleep. The study was approved by the Clinical Research Ethics Committee of the University of the Balearic Islands (Palma de Mallorca, Spain). fMRI data acquisition was performed with a GE Medical Systems Signa HDx 3 T scanner using echo-planar sequences, 240 volumes were acquired with a TR of 2500 ms, TE of 35 ms, and 90° flip angle. Thirty-six slices of 64 × 64 dimensions were obtained with a field of view of 200 mm and slice thickness of 3 mm. Structural images consisted of a T1-weighted scans of 176 × 512 × 512 voxels, with a TR of 7176 ms, TE of 3150 ms, flip angle of 12°, FOV 240 mm and slice thickness of 1 mm. Preprocessing of BOLD signal was performed using FMRIB Expert Analysis Tool (Jezzard et al., 2001)², including motion correction using MCFLIRT, slice-timing correction using Fourier-space time series phase-shifting, non-brain removal using BET and spatial smoothing using a Gaussian kernel of full-width-half-maximum 5 mm. Brain images were normalized to standard space with FLIRT using the MNI 152 template and resampled to 4 mm × 4 mm × 4 mm resolution. Resting functional data was filtered with a zero lag finite impulse response band pass filter (0.01–0.1 Hz; Cordes et al., 2000, 2001). fMRI data used for the Figures 2D,E, as well as the preprocessing steps, were the same than in (Tagliazucchi et al., 2010a) and (Tagliazucchi et al., 2010b). Melodic was used for the PICA calculation of RSN (Beckmann and Smith, 2004) in Figure 2 as well as for denoising motion artifacts.

The point process is defined by the sequences of time points at which the BOLD signal crosses a given threshold from below. Formally, the problem is defined in an autonomous system as

where the dot denotes time derivative and x→∈Rd. Let y=h(x→) be a scalar observable function (such as the BOLD signal, for instance) and consider the plot of y versus t. The times at which y(t) upward (or downward) crosses some predetermined threshold y = y_c determine a sequences of time points {t_k}_k=1,N which defines the so-called point process (see Figure 5).

The Poincaré section of any given dynamical system reduces a d-dimensional continuous time description into an associated (d-1)-dimensional discrete map by finding the intersections of trajectories in phase space with a surface S transverse to the flow. If x→k∈S denotes the k^th intersection of the trajectories with S, a Poincaré map is defined as

A sequence of crossing times may be taken as a Poincaré section as representative as the typical phase space coordinates, when certain conditions are satisfied: Let γ→(x→0,t) a solution of Eq.(1) in an open interval I = |t₀, T| and let |γ→˙(x→0,t)|≠0 for all t ∈ I. In terms of the underlying dynamical system, this means that the dynamics is not in a fixed-point or equilibrium of the system (as one can assume for the BOLD signal and neural dynamics in general). Under this condition, the arc-length, defined as s(x→0,t)=∫t0t|γ→˙(x→0,t′)|dt′, is a suitable observable in the sense of embedding theory (Hegger and Kantz, 1997) and it is possible to reconstruct the attractor of the system by measuring line segments s̃(ti+1,ti)=s(x→0,ti+1)-s(x→0,ti).

The derivative of s with respect to t, dsdt=|γ→˙(x→0,t)|=|f→(x→)|, allows us to rewrite equation 1 as a set of differential equations in s:dx→∕ds=f→∕|f→| and this gives the possibility to use the time t as an usual variable, which is no longer the independent parameter but an usual coordinate as the variable x→. Thus, the embedding theorems also apply to time or properly defined time sequences and it is possible to reconstruct the attractor (i.e., the full properties of the underlying dynamical system) from this sequence of times (i.e., the point process), as discussed in Hegger and Kantz (1997).

The fMRI BOLD signal was de-convoluted using the deconvlucy.m function from Matlab³. For the Hemodynamic Response Function (HRF) standard parameters were those provided in the SPM8 package⁴. The deconvolution function follows the Lucy-Richardson algorithm (Richardson, 1972; Lucy, 1974) which converges to the maximum likelihood estimate of the de-convoluted process assuming a Poissonian source of noise. The results of the deconvolution, as in the example presented in Figure 1C, are impulse-like signals, whose underlying neural mechanisms are beyond the scope of the present work. Note that the deconvolution of all voxel’s time series is in principle possible, however its numerical implementation is several orders of magnitude less efficient than the simple thresholding used here.

To construct the conditional rates reported in Figure 2 the point process is defined at a seed location and at the targets throughout the entire brain. Figure 6 illustrates the basic procedure. The BOLD signal is extracted from a seed region (top trace) and the points (arrows and vertical dashed lines) are defined by the crossings at 1 SD (horizontal dashed lines). Every time the signal at a target region crosses the threshold (asterisks) up to 2 time steps later than in the seed, the rate at the target is increased in one unit. This rate is normalized by the number of points in the seed. The top panel shows the location of the seed and of the two example targets, as well as the resulting average conditional rates maps (left) and DMN obtained from PICA (right). Medium panels show the BOLD signal at the seed and at the two example target regions. A similar procedure was used in (Tagliazucchi et al., 2010b) where the resting BOLD event triggered averages (rBeta) were calculated at similar seed and target regions. Table 1 contains the seed coordinates used to reproduce the RSNs.

Spatial clusters of activated voxels were identified using an algorithm implemented in MATLAB, based on the detection of connected components in a co-activated first neighbors graph. Clusters’ fractal dimension was calculated using a standard box-counting algorithm. Avalanches were defined (similar as in sandpile models, and others (Bak, 1996; Jensen, 1998)) as starting with the isolated activation (i.e., not by any of its neighbors) of a previously inactive voxel (or group of voxels), continuing while at least one contiguous voxel is active in the next time step and otherwise ends. The avalanche tracking algorithm implemented in this work uses as a criteria for avalanche membership the non-empty intersection with a previously identified cluster of the avalanche at a previous time. This is able to resolve shrinking and expanding of clusters, translation, and division, whenever there is spatial overlap at subsequent times.

To detect contiguous clusters of activated voxels (defined as those crossing the threshold), for each time step, the problem was reduced to the detection of connected components in a suitably defined graph or network. More precisely, for each volume, a graph was constructed having each voxel as a node, and two nodes connected with a link if they were both activated (BOLD signal above 1 SD) and also first neighbors in the spatial sense. The connected components of this graph correspond to clusters of contiguous activated voxels isolated from other similarly defined clusters.

In simple terms an avalanche starts with the activation of a previously inactive voxel, follows while in the next time step at least one contiguous voxels is active and otherwise ends. The avalanche detection algorithm is based on the connected cluster decomposition. Clusters are followed during different volumes, belonging to the same avalanche if they have spatial intersection during consecutive times. Formally, the algorithm is as follows: Let Cit be the i-th cluster at time t. We consider a cluster i₀ starting an avalanche at time t₀ if for all j, Cjt0-1∩Ci0t0=∅(i.e., no clusters were present in that region of the brain at the previous time step). An id is assigned to this avalanche and the same id is assigned to all clusters intersecting this cluster at the following time, this is all clusters i such that Ci0t0∩Cit0+1≠∅. The same procedure is applied recursively to all clusters satisfying the former condition until no more intersections are found. When this happens, all clusters labeled with this id constitute the avalanche.

In this work, the only free parameter used in the definition of the point process is the threshold. In this sense, it is important to know how the main spatiotemporal statistical properties of the point process dynamics, namely cluster size distributions, avalanche size and duration distributions depend on threshold values. Figure 7 shows that these results are robust against changes in threshold over a wide range of choices.