Spike Avalanches Exhibit Universal Dynamics across the Sleep-Wake Cycle

Abstract

Background

Scale-invariant neuronal avalanches have been observed in cell cultures and slices as well as anesthetized and awake brains, suggesting that the brain operates near criticality, i.e. within a narrow margin between avalanche propagation and extinction. In theory, criticality provides many desirable features for the behaving brain, optimizing computational capabilities, information transmission, sensitivity to sensory stimuli and size of memory repertoires. However, a thorough characterization of neuronal avalanches in freely-behaving (FB) animals is still missing, thus raising doubts about their relevance for brain function.

Methodology/Principal Findings

To address this issue, we employed chronically implanted multielectrode arrays (MEA) to record avalanches of action potentials (spikes) from the cerebral cortex and hippocampus of 14 rats, as they spontaneously traversed the wake-sleep cycle, explored novel objects or were subjected to anesthesia (AN). We then modeled spike avalanches to evaluate the impact of sparse MEA sampling on their statistics. We found that the size distribution of spike avalanches are well fit by lognormal distributions in FB animals, and by truncated power laws in the AN group. FB data surrogation markedly decreases the tail of the distribution, i.e. spike shuffling destroys the largest avalanches. The FB data are also characterized by multiple key features compatible with criticality in the temporal domain, such as 1/f spectra and long-term correlations as measured by detrended fluctuation analysis. These signatures are very stable across waking, slow-wave sleep and rapid-eye-movement sleep, but collapse during anesthesia. Likewise, waiting time distributions obey a single scaling function during all natural behavioral states, but not during anesthesia. Results are equivalent for neuronal ensembles recorded from visual and tactile areas of the cerebral cortex, as well as the hippocampus.

Conclusions/Significance

Altogether, the data provide a comprehensive link between behavior and brain criticality, revealing a unique scale-invariant regime of spike avalanches across all major behaviors.

About the Authors:

Tiago L. Ribeiro

Affiliation: Department of Physics, Federal University of Pernambuco (UFPE), Recife, Pernambuco, Brazil

Mauro Copelli

Affiliations Department of Physics, Federal University of Pernambuco (UFPE), Recife, Pernambuco, Brazil, Neuroscience Graduate Program, Federal University of Rio Grande do Norte (UFRN), Natal, Rio Grande do Norte, Brazil

Fábio Caixeta

Affiliations Neuroscience Graduate Program, Federal University of Rio Grande do Norte (UFRN), Natal, Rio Grande do Norte, Brazil, Edmond and Lily Safra International Institute of Neuroscience of Natal (ELS-IINN), Natal, Rio Grande do Norte, Brazil

Hindiael Belchior

Affiliations Edmond and Lily Safra International Institute of Neuroscience of Natal (ELS-IINN), Natal, Rio Grande do Norte, Brazil, Department of Physiology, Federal University of Rio Grande do Norte (UFRN), Natal, Rio Grande do Norte, Brazil

Dante R. Chialvo

Affiliation: Department of Physiology, Northwestern University, Chicago, Illinois, United States of America

Miguel A. L. Nicolelis

Affiliations Neuroscience Graduate Program, Federal University of Rio Grande do Norte (UFRN), Natal, Rio Grande do Norte, Brazil, Edmond and Lily Safra International Institute of Neuroscience of Natal (ELS-IINN), Natal, Rio Grande do Norte, Brazil, Department of Neurobiology, Center for Neuroengineering, Duke University Medical Center, Durham, North Carolina, United States of America, Department of Biomedical Engineering, Duke University, Durham, North Carolina, United States of America, Department of Psychological and Brain Sciences, Duke University, Durham, North Carolina, United States of America

Sidarta Ribeiro

* E-mail: [email protected]

Affiliations Neuroscience Graduate Program, Federal University of Rio Grande do Norte (UFRN), Natal, Rio Grande do Norte, Brazil, Edmond and Lily Safra International Institute of Neuroscience of Natal (ELS-IINN), Natal, Rio Grande do Norte, Brazil, Department of Physiology, Federal University of Rio Grande do Norte (UFRN), Natal, Rio Grande do Norte, Brazil

Introduction

Several recent studies have revealed that neuronal populations exhibit a type of activity termed neuronal avalanches, characterized by the occurrence of bursts of activity that, despite their wide variation in sizes and durations, still follow precise statistical properties. The main signature of these avalanches is their size distribution, which decays as a power law , with exponents α around 1.5. Two features of this type of distribution are particularly noteworthy. First, they are scale-invariant: if we know how likely it is to observe a burst of size s and ask how likely it would be to observe a burst of size k times this size, the answer is that the relative likelihood is , which is independent of s (i.e. changing the scale at which sizes are measured does not change the relative abundance of burst sizes). Second, such power law distributions are heavy-tailed, which implies that it does not make sense to speak of a typical (or characteristic) burst size (note that for α<2 the variance diverges). In other words, fluctuations rule the underlying dynamics.

Scale-invariant neuronal avalanches have been first observed in cell cultures and slices [1], but recent studies of anesthetized rats [2] and awake restrained monkeys [3] indicate that they also occur in intact brains. These results are important because scale-invariance in neuronal dynamics may provide a long-sought connection between brain functioning and self-organized critical (SOC) systems [4], [5]. These are systems that can self-tune to a balanced (critical) state, precisely at the transition between a (subcritical) regime of inactivity and one of (supercritical) runaway activity.

The hypothesis that tuning a biological system to a critical state would render it somehow optimal has a long history [6]. The underlying idea is that a system tuned to criticality presents a richer dynamical repertoire, being therefore able to react (i.e. process information) to a wider range of challenges (environmental or other). The experimental evidence in this direction ranges from gene expression patterns in response to stimulation of single macrophages [7] to collective ant foraging [8].

In particular, criticality was also suggested to play an important role in the brain [4], [5]. From the theoretical side, numerous results show that scale-invariant dynamics provide functionally desirable features for the behaving brain, such as optimal computational capabilities [9], information transmission [1], size of memory repertoires [10] and sensitivity to stimuli [11], [12]. Experiments, on the other hand, have both confirmed theoretical predictions [13] as well as provided evidence of scale-free dynamics that still need to be better explored from the modeling point of view. These include results obtained at the whole-brain scale, where functional networks compatible with a critical brain were observed via functional magnetic resonance imaging (fMRI) [14], [15], magnetoencephalographic (MEG) [16] and electroencephalographic (EEG) [17] data.

At a smaller scale, measurements of neuronal avalanches were mostly restricted to local field potentials (LFPs) recorded in vitro or anesthetized in vivo conditions [1], [2], [10], [13]. In these cases, consecutive avalanches are usually well separated in time, their duration typically lasting much less than the interval among them (this separation of time scales being a hallmark of SOC models). In this sense, avalanches have been interpreted as elementary collective excitations that occur at base level as ongoing activity [1], but constitute nevertheless stable templates of spatio-temporal activity with a repertoire potentially recruitable by behavior [10].

While criticality is well established for LFP data [1]–[3], [10], the findings regarding spike activity remain unsettled. Spike avalanches were studied in two very different scenarios. Recordings from dissociated neuronal cultures [18], [19] yielded distributions similar to those previously reported for LFPs. Recordings in vivo from the cat parietal cortex, however, led to size distributions that seemed incompatible with power laws [20]. In fact, this study suggested that the statistical properties observed in LFPs might be attributed to a nontrivial filtering caused by the complexity of the extracellular medium. Given this state of affairs, it is therefore crucial to understand spike avalanche dynamics in non-reduced preparations, across the full range of natural behaviors.

To address this issue, extracellular spiking activity was recorded with multielectrode arrays (MEA) from multiple brain regions of adult rats freely cycling across the major behavioral states: waking (WK), slow-wave sleep (SWS) and rapid-eye-movement sleep (REM). Behavioral sorting of these major states was automatically implemented by a computer program for spectral analysis of LFPs simultaneously recorded from the same electrodes [21], then confirmed by inspection of video recordings. Data were acquired at different stages: before, during and after exposure to novel objects, respectively referred to as PRE, EXP and POST periods (Fig. S1 and Methods). For comparison with this freely-behaving group (FB, n = 7), animals were recorded under deep anesthesia (AN, n = 7) (see Methods for details). MEAs were targeted to the primary visual (V1) and somatosensory (S1) areas of the cerebral cortex, which receive direct inputs from thalamic relays connected to the eyes and facial whiskers, respectively [22]. Recordings were also performed from the hippocampus (HP), a subcortical structure related to sensory integration, exploratory behavior and memory formation [23].

Results

As previously defined (Ref. [1]; see Methods), spike avalanches were extracted from the spike time series (Fig. 1a) and temporally divided in rate-normalized bins (Fig. 1b). The sizes of spike avalanches varied widely over time, spanning more than two orders of magnitude (Fig. 1c). To characterize this variation, we calculated the probability P(s) = Prob[size = s]. We obtained P(s) separately for the different brain regions (V1, S1, HP), behavioral states (WK, SWS, REM) and stages of the experiment (PRE, EXP and POST). For each rat, bin widths were separately calculated for each of these 27 combinations, ranging from 2 to 50 ms (Table S1) and thus reflecting the diversity of the number of sampled neurons (Table S2) and firing rates. With rate-normalized bins, these widely different situations could be cross-compared on fair grounds.

[Figure omitted. See PDF.]

Figure 1. Obtaining spike avalanches from raw data.

(a) Raster plot of neuronal spikes and LFPs traces recorded from a freely-behaving rat undergoing the three major behavioral states (first three panels, 2 s windows) or anesthesia (last panel, 3 s window). Note the clearly rhythmic spiking activity coupled with LFP oscillations during anesthesia. (b) To understand how spike avalanches were defined, consider a 40-ms excerpt sliced in 1.3-ms time bins. Adding up all spikes within each bin, one obtains a sequence of avalanches of sizes 2, 1, 2, 1, 4, 2, 4, 2, 1, and 2. To account for firing rates variations across behavioral states, experimental stages and brain structures, and to control for neuronal ensemble size, bin width corresponded to the average inter-event interval (IEI) in each dataset. (c) Time series of spike avalanche sizes in S1 cortex. Horizontal arrow shows waiting time between consecutive avalanches of minimum size sc.

https://doi.org/10.1371/journal.pone.0014129.g001

Pooling avalanches from all FB rats results in very similar size distributions, either across the sleep-wake cycle (for a given stage of the experiment) or across the stages of the experiment (for a given behavioral state), regardless of the brain area (Fig. 2). Differently from what has been observed previously [1], [2], [3], [13], [19], however, the FB size distributions were not compatible with (and decay faster than) a power law. They resemble spike avalanche size distributions recorded from the cat cortex [20]. However, contrary to what was proposed for those distributions, the FB distributions did not fall off exponentially. Rather, they were well fit by a lognormal distribution: (Fig. 2). When the data were surrogated by shuffling inter-spike intervals for each neuron, the distribution tails were substantially shortened for all states, experimental conditions and brain areas compared, reflecting the destruction of the largest avalanches (Fig. 3a; see also Fig. S3).

[Figure omitted. See PDF.]

Figure 2. Avalanche size distributions were stable across behavioral states, experimental stages and brain areas.

Avalanche size distributions. Each row represents a brain region, while columns distinguish stages of the experiment. For each combination, the three behavioral states are shown in a double logarithmic plot for data pooled from representative rats (single animal WK distributions in gray). Lines represent lognormal fits.

https://doi.org/10.1371/journal.pone.0014129.g002

[Figure omitted. See PDF.]

Figure 3. Size distributions from different conditions are not significantly different and deviate substantially from the surrogated data size distributions.

(a) Size distributions for original (colors) and surrogated (black) data, for three different conditions. Lines represent lognormal fits. Although surrogated spike trains have precisely the same firing rates as original data, larger avalanches consistently occur less frequently. (b) Comparison between cumulative size distributions for different cases. In gray, the QQ-plot for the same curves. P-value calculated by a KS test; note that the distributions are very similar despite failing the statistical test.

https://doi.org/10.1371/journal.pone.0014129.g003

These results were supported by various statistical analyses. The fits of the pooled size distributions were subjected to the Kolmogorov-Smirnov (KS) test at a p = 0.05 significance level: 23 out of 27 distributions from Fig. 2 were compatible with the fitted lognormals, and none was compatible with a power law or an exponential. The KS test was also employed to compare pairs of FB size distributions in two scenarios: 1) from different stages of the experiment (but the same brain region and behavioral state) and 2) from different behavioral states (but the same brain region and stage of the experiment). The fractions of equivalent comparisons were 36% and 22%, respectively. Note, however, that the KS test is extremely stringent, as it compares two distributions only on the basis of the maximum difference between them. We therefore employed a graphical method to better illustrate the similarity of the distributions. The Q-Q plots (Fig. 3b) display an excellent agreement of different distributions, even when the comparison fails the KS test (see also Fig. S5). Further statistical analysis of this issue can be found in the Supporting Information Text S1 and in Figs. S5 and S6.

Previous work has shown that avalanche size distributions can change considerably depending on whether a critical system is fully or partially sampled [24], [25]. Evidently, any avalanche size distribution obtained for neuronal ensembles recorded with MEAs corresponds to a severe undersampling of the total population of neurons in a given brain. To further investigate this issue, we built a probabilistic excitable cellular automaton model [11] tuned near the critical state (Fig. 4a), and deliberately undersampled it with a spatial structure equivalent to that of our MEAs (Fig. 4a, inset; see also Methods). Despite the fact that the model was critical by construction, simulated spike avalanches exhibited lognormal-like size distributions when undersampled (Fig. 4a, red triangles), in excellent agreement with the in vivo data (Fig. 4a, blue triangles). In contrast, the distribution of avalanche sizes using all neurons in the model lattice obeyed a power law (Fig. 4a, black circles). Minor deviations for the smallest samples were also observed in the model. The inset in Fig. S2 shows that this “saturation” effect, which increases as the system is progressively undersampled, was also well reproduced by the model (compare the 1–2% undersampling shown in the inset of Fig. S2 with Fig. 4a, that represents a 4% undersampling).

[Figure omitted. See PDF.]

Figure 4. Size distributions from undersampled critical systems interpolate between lognormals and power laws.

(a) Size distributions for model (red triangles: undersampling; circles: full sampling) and FB data (blue triangles). Lines are lognormal and power law fits. Inset: model lattice (black dots) and sampled sub-lattice that mimics the configuration of the neurons recorded by the MEA (red triangles). (b) Size distributions from AN animals are well fit by power laws. Inset: size distributions for different levels of undersampling using the model modified to simulate anesthesia. (c) Size distributions from three AN rats. From bottom to top, curves go from deeply anesthetized to fully recovered (each curve corresponds to 30–60 minutes of data). Red lines represent the best fit for the bottom (power law) and top (lognormal) distributions.

https://doi.org/10.1371/journal.pone.0014129.g004

This agreement between simulations and FB data is largely insensitive to changes in model parameters. For instance, the value of the simulated stimulus rate, h, can be changed by many orders of magnitude without altering the results. The only constraint is an upper limit, above which the network firing rate will be so high that the calculated time bin will be less than one time step, thus collapsing all avalanches onto a single one (for the 1024 model neurons of our simulations, we found that this upper limit is of the order of 10−3 ms−1). On the other hand, there is no lower limit for h. In fact, simulations in which avalanches are created by exciting a single neuron and waiting until network activity dies out (which corresponds to the limit h = 0) leads to the same lognormal distributions observed when the system is undersampled (see Methods).

Avalanche size distributions obtained from the FB group were equivalent for neuronal ensembles of very different sizes and, for a given size, insensitive to changes in the choice of the specific sampled neurons. Fig. S2 shows that reducing the number of neurons sampled in our original datasets yielded similar lognormal distributions. This was expected because the distributions observed in Fig. 2 are all very similar, despite the wide range of number of neurons sampled across animals (Table S2). For the smallest samples, slight non-monotonic deviations were observed. They were expected because the decrease in the number of neurons leads to larger time bins. This by its turn leads to fewer avalanches per time unit, yielding poorer statistics and increased variability.

Evidently, the AN group is subjected to the same undersampling constraints imposed by the MEA arrays on the recordings from FB animals. Measuring spike avalanches from animals deeply anesthetized with ketamine-xylazine, we obtained size distributions very similar to power laws (Fig. 4b), with exponents α comparable to those observed for LFPs in brain slices [1] and for spikes in dissociated neuronal cultures [19] as well as anesthetized cats [26] (AN exponents vary from 0.9 to 1.8, while in vitro exponents fluctuate around 1.5).

In order to deepen the understanding of the differences in avalanche size distributions between FB and AN data, we modified the cellular automaton model by applying weak, sparse and periodic stochastic drive. In contrast to the initial model, clearly the modified model was less affected by undersampling, yielding size distributions similar to a power law even when the system was not fully sampled (Fig. 4b, inset). In agreement with this scenario, size distributions for V1 and S1 neurons gradually returned to their lognormal-like shape as the animals recovered from anesthesia (Fig. 4c). The main difference between the distributions for anesthetized and recovered conditions was the cutoff in the bottom curve, which is expected because there is a clear separation of time scales in the AN data (Fig. 5b; compare with FB data in Fig. 5a), so each neuron typically spikes at most once per avalanche. Moreover, comparing lognormals and power laws via the normalized squared sum of the residuals (Nred2, see Methods), we observed that in all three cases of panel 4c the “recovered” distributions were better fit by lognormals (Nred2 = 0.08 vs 0.71 for AN2, 0.24 vs 0.33 for AN4 and 0.18 vs 0.52 for AN5), while the “anesthetized” distributions were better fit by power laws (Nred2 = 0.16 vs 0.34 for AN2, 0.19 vs 0.32 for AN4 and 0.33 vs 0.46 for AN5). The effect is very consistent across animals and demonstrates that the size distributions of spike avalanches in the FB and AN conditions have fundamental statistical differences.

[Figure omitted. See PDF.]

Figure 5. Avalanche duration and inter-avalanche interval distributions.

(a) Distributions for three different rats during SWS sleep. Inset: for the same rat, the same distributions for WK and REM. (b) Distributions for an animal from the AN group, during anesthesia (orange) and after recovery (brown). Note the separation of time scales between avalanche durations and inter-avalanche intervals during anesthesia.

https://doi.org/10.1371/journal.pone.0014129.g005

To further investigate signatures of criticality in FB rats, we analyzed the power spectra of the avalanche size time series (Fig. 1c), which are consistent with 1/f -like behavior (Fig. 6a, compare with distribution from surrogated data). A detrended fluctuation analysis (DFA) showed that the root mean squared deviation from the detrended time series increases as a power of the window width (Fig. 6b). Notice that an exponent close to one corresponds to a 1/f spectrum [27]. In contrast, AN data yielded a power spectrum with a Poisson-like decay (Fig. 6a, bottom) and no clear power law regime in the DFA analysis (Fig. 6b, bottom).

[Figure omitted. See PDF.]

Figure 6. Statistical fingerprints of criticality in spike avalanches recorded from freely-behaving animals.

(a) Power spectrum of the avalanche size time series for two FB rats and one AN rat. Though conservative, the shuffling procedure destroys the long-range correlations characterized by the 1/f spectrum seen for FB data. (b) Root-mean squared fluctuation F of the detrended avalanche size time series versus window width n, for two FB rats and one AN rat. In all cases, α denotes the exponent of a fitted power law. DFA exponents close to one are compatible with 1/f power spectrum. Note the poor quality of the power law fit for DFA AN data.

https://doi.org/10.1371/journal.pone.0014129.g006

Finally, we studied the probability density D(t;sc) of waiting times t between consecutive avalanches of size larger than or equal to sc (Fig. 1c). Clearly, larger values of sc increase the probability of longer waiting times (Fig. 7a). However if for each sc we plot as a function of , where is the mean interval between avalanches of size at least sc, the rescaled curves collapse reasonably onto a single function (right plot of Fig. 7a), thus . In other words, the recurrence of avalanches of different minimum size sc is governed by a single function. The collapse occurs in different brain regions, stages of the experiment and behavioral states (Fig. 7b, top), but not for anesthesia (Fig. 7b, bottom). Furthermore, since the rescaled axes are dimensionless, it is possible to directly compare results from different rats. When pooled, results for all seven FB animals during WK collapse over six orders of magnitude (Fig. 8a, right plot), and the same holds for the other natural behavioral states assessed (Fig. 9a). In contrast, a similar scaling regime does not apply to the AN data (Fig. 8b).

[Figure omitted. See PDF.]

Figure 7. Waiting time distributions for different minimum avalanche sizes collapse onto a single scaling function for each FB animal (but not for AN).

(a) Probability density of avalanche recurrence times (without rescaling in the left panel; rescaled in the right panel) for one FB rat. (b) The same collapse for different animals (FB top, AN bottom). Note that the collapse under this kind of scaling occurs for all major natural behaviors, stages of the experiment and brain areas, but not during anesthesia.

https://doi.org/10.1371/journal.pone.0014129.g007

[Figure omitted. See PDF.]

Figure 8. Data from all FB animals have a similar scaling function, which breaks down during anesthesia.

(a) Regular and rescaled waiting time distributions for all FB rats. The scaling function is well fit by a double power law (see also Fig. 9). (b) The same distributions for AN data show no sign of collapse under the same rescaling procedure. Note the presence of a characteristic waiting time for a range of minimum avalanche sizes.

https://doi.org/10.1371/journal.pone.0014129.g008

[Figure omitted. See PDF.]

Figure 9. The scaling function is very similar across all major natural behavioral states and brain regions.

(a) Rescaled waiting time distributions obtained from all FB rats for each behavioral state and brain region (all stages of the experiment included). Colors (black) represent original (surrogated) data. The lines represent the best fit in each case. (b) Comparison between double power law (DPL) and exponentially decaying gamma (EdG) fits. The quality of the fit increases as the sum of square residuals Nred2 decreases, showing that the DPL yields the best fit for all cases. (c) Scatter plot of the exponents of the DPL fit for all distributions in panel (a). Note that the dispersion is significantly larger for surrogated data.

https://doi.org/10.1371/journal.pone.0014129.g009

We fitted double power laws (DPL) and exponentially-decaying gamma functions (EdG) to both original and surrogated FB data. As shown in Fig. 9b, the DPL yielded the least Nred2 in all cases, even for surrogated data (solid lines in Figs. 8a and 9a). Note that, for different behavioral states, the difference Δ in Nred2 values between original and surrogated data decreases with decreasing sampling time: ΔWK (T = 62,500 s)>ΔSWS (T = 39,690 s)>ΔREM (T = 5,530 s). This can be explained by the surrogating method employed, which shuffles inter-spike intervals only within each window of a given behavioral state (see Methods). The exponents for original FB data were clustered around their mean values, and , in contrast with exponents for surrogated data, whose values had consistently larger variation: and (Fig. 9c).

Discussion

Avalanche size distributions in freely-behaving animals

The size distributions obtained from the FB group are remarkably similar across sleep-wake states, experimental stages and brain areas (Fig. 2). This is surprising, given that brain dynamics changes considerably in these different conditions. The behavioral states are not only characterized by different LFP spectral features (Fig. 1a) [21], but also the exposure to novel objects leads to very significant changes in firing rates [23]. The results are not a simple effect of firing rate normalization owing to our binning procedure: the size distributions are heavy-tailed, in the sense that large avalanches occur more frequently than would be observed for spike trains with identical mean firing rates but uncorrelated (see surrogated data in Figs. 3a, S3 and S4).

What could be the origin of these non-power law but heavy-tailed distributions? Insight into this issue came from sandpile and forest fire models of self-organized criticality, known to exhibit power law size distributions. Previous works showed that power laws in neuronal avalanches fail to emerge when the system is sparsely sampled [24], [25]. This occurs because in these systems the observables of interest (e.g. size of avalanches or forest fires) are derived from the spatio-temporal activity of a much more complex underlying dynamics. The MEA setup inevitably misses most of the spike activity in the implanted brain region. Moreover, since MEAs allow for spike sampling from multiple neurons in the vicinity of each electrode, a recording bias towards relatively small avalanches ensues. More specifically, an avalanche created somewhere in the targeted brain region may: 1) not be detected at all by the electrodes, 2) be detected by some of the electrodes as one smaller avalanche, or 3) be detected by different groups of electrodes as more than one smaller avalanche. The combination of these factors modifies the density distribution, morphing a putative power law into a lognormal.

In accordance with these results, FB distributions were lognormals and were insensitive to changes in the number or position of the neurons sampled by the MEA. A likely explanation for this invariance is that the full MEA already represents a substantial undersampling of the targeted brain region. Evidently, undersampling is present in any current large scale neuronal recording, therefore raising the question of how power law avalanche size distributions previously found in the literature could emerge under this constraint. In this regard, FB neural activity contrasts sharply with that of reduced preparations in the degree of redundancy across electrodes. Cultures and slices exhibit high synchronization, manifested in intermittent spike bursts [28], [29]. For this reason, spike avalanches in dissociated neuronal cultures display robust power laws even when sampled with sparse MEAs [19]. In slices and anesthetized intact animals [30], large and slow (<2 Hz) LFP oscillations constrain neurons to up and down states, which correspond respectively to the depolarizing and hyperpolarizing phases of the oscillation. Since spikes tend to occur during up states, anesthesia leads to discrete bursts of spiking activity across vast neuronal ensembles, in synchrony with LFP oscillations (Fig. 1a). This generates short, non-overlapping avalanches, in a temporal pattern similar to those observed in vitro [1], [19]. We therefore hypothesized that spike avalanches recorded from anesthetized animals exhibit size distributions more similar to power laws, in spite of the use of sparse MEAs.

Avalanche size distributions in anesthetized animals

The results showed that network dynamics during anesthesia effectively overcome the undersampling effects seen in FB. We obtained power law size distributions from the AN data, in accordance with a recent study of spike and LFP avalanches in the visual cortex of cats also under ketamine-xylazine anesthesia [26]. This can be explained by three immediate consequences of spike burst synchronization during anesthesia: First, the lower firing rates typical of AN lead to time bins which are large enough to ensure that a large avalanche will not be artificially split in smaller avalanches due to brief periods of silence; second, it implies a separation of time scales between the avalanche durations and the intervals between them (a ubiquitous feature in SOC systems), making it less likely that different avalanches will artificially merge (compare Figs. 5a and b); third, large-scale synchronization leads to redundancy in the MEA, increasing the probability that two or more neurons far from each other will fire within the same avalanche (thus attenuating the undersampling effect).

Anesthesia was simulated in the cellular automaton model with the use of weak, sparse and periodic stochastic inputs in order to mimic higher spike correlations and enhanced rhythmic activity that are characteristic of the anesthetized state [31], [32]. The particular anesthetic drugs used in the experiment comprise a two-fold action: ketamine decreases the net levels of excitation by antagonizing glutamatergic N-methyl-d-aspartate (NMDA) channels [33], [34], whereas xylazine decreases noradrenergic modulation by activating a2 adrenergic receptors [35]. Since the pattern of noradrenergic modulation in the telencephalon is globally widespread but locally scattered [36], xylazine can be presumed to affect neuronal activity in a sparse manner. The interaction of gabaergic and non-NMDA glutamatergic systems spared by ketamine and xylazine generates slow membrane potential oscillations that drive cortical neurons periodically [30]. Likewise, reduced preparations (culture, slices) deprived from neuromodulatory inputs develop slow waves of activity as a result of the interaction between glutamatergic and gabaergic circuits [28]. The relationship of these waves to in vitro avalanches is still unclear.

The deviations from the α = 1.5 power law exponent under our AN conditions do not show under anesthesia with urethane [2], a broad-action anesthetic that potentiates gabaergic, glycinergic and nicotinic cholinergic receptors, while inhibiting NMDA and non-NMDA glutamatergic receptors [37]. Neural processing can be quite different under ketamine-xylazine and urethane, as reviewed in [38]. Further investigation is required to determine which anesthetic best mimics the dynamics of up and down states that characterizes in vitro preparations [1].

Temporal signatures of criticality in freely-behaving animals

All existing evidence indicates that a size distribution undersampled from a system that follows a power law will have a different shape, most importantly a reduction of the weight of its tail [24], [25]. Therefore it is likely that power laws, which are the most commonly sought signature of critical behavior, may not be directly detected in the case of spike avalanche size distributions in FB animals, owing to the inevitable undersampling of the MEA method. This does not imply, however, that other signatures of criticality cannot be found in spike data recorded from FB animals. We searched for alternative statistical fingerprints that could test the hypothesis that the freely-behaving brain operates near a critical regime. We found that FB distributions display 1/f -like behavior (Fig. 6a), which indicates that the system has long-term correlations, in agreement with other SOC systems [39]. Consistently, a detrended fluctuation analysis (DFA) revealed signatures of long-term correlations for natural behavioral states, but not for anesthesia (Fig. 6b).

Finally, we verified that the recurrence of avalanches of minimum size sc is governed by a scaling function. In other words, a single function describes recurrence times from a few milliseconds to hundreds of seconds, for any size sc. The kind of scaling we obtained (dependence of the scaling function only on ) is akin to what is observed in self-organized critical systems, such as solar flares [40], fractures [41] and forest fires [42]. Moreover, we found that the scaling function is well fit by a double power law, with remarkably similar exponents across all brain regions and behavioral states (Figs. 9a and 9c), suggesting the existence of universal mechanisms underlying the dynamics of spike avalanches in the brain. The exponent values are very similar to those observed for earthquakes, where double power laws have also been observed [43], [44]. This particular scaling function for the waiting time distribution can be interpreted following Ref. [43]. For “shorter” waiting times, the distribution is dominated by the exponent α∼0.7, implying that consecutive avalanches are correlated (like in the Omori law for earthquakes [45], which has recently been observed for neuronal avalanches [46]). For “longer” waiting times, the distribution is dominated by the much larger exponent β∼3, a regime in which consecutive avalanches would be independent. The meaning of “shorter” and “longer” waiting times, however, is not absolute, but size-dependent. This absence of a characteristic time scale is suggestive of a critical system.

Criticality during natural behavior or anesthesia?

With regard to size distributions during anesthesia, the existence of power laws supports spike criticality despite undersampling. This likely means that in the anesthetized brain the local connectivity (at the MEA scale) is preserved, i.e. the main avalanche pathways remain active. On the other hand, anesthesia disrupts the temporal dynamics of the system, transforming the critical temporal dynamics seen in freely-behaving animals (as indicated by Fourier, DFA and waiting time collapse analyses) into a non-critical temporal process characterized by rhythmic activity, a typical inter-avalanche interval, and impoverished temporal repertoire.

The results in freely-behaving animals indicate that a single mechanism produces small and large spike avalanches, as well as short and long inter-avalanche intervals, during WK, SWS and REM. This finding is far from trivial, because behavioral state variations are associated with marked changes in membrane resting potentials, neuronal firing rates, and LFP oscillations [21], [30], [47]. Our results provide strong evidence, at the level of neuronal ensembles, that the behaving brain operates near a temporally complex regime that is maintained across all major behaviors but collapses during ketamine-xylazine anesthesia.

Could the timescale separation in the AN data be sufficient to explain the differences between AN and FB states? According to the results obtained in the model for the FB data, the answer is no. When very low rates of external stimulation were applied, and therefore when avalanches were well separated, we still observed the undersampling effect. In fact, note that in Ref. [24] the models used implied an infinite separation between avalanches, but the undersampling effect was still present (the same holds for our model). The fundamental distinction between anesthesia and any natural behavioral state has been recently underscored by the discovery that, while some comatose patients are capable of learning, subjects anesthetized with propofol, a gabaergic agonist and sodium channel blocker [48]–[51], are not [52].

Clearly, our simple model is unable to account for the scaling function observed in the waiting time distributions of the FB group. We are unaware of any model with neurobiological plausibility that simultaneously reproduces scale-free size distributions and critical time-domain measures (e.g. 1/f spectra and DFA). As shown in previous work on a critical sandpile model, a nontrivial drive can modify substantially the resulting statistics of waiting time intervals [53]. Given the complex input to which any brain region is subjected, the modeling of such a system remains a major challenge.

We also observed that spikes recorded before and after the exploration of novel objects showed similar avalanche statistics. Exposure to novel objects is a procedure known to increase firing rates, induce plasticity factors and promote dendritic sprouting in the cerebral cortex and hippocampus, leading to memory formation and learning of object identity [23], [54]–[56]. Our results argue directly against the notion that the encoding of new memories is produced by gross changes in avalanche regime. Rather, the data support the view that behaving brains are optimized for the encoding of memory patterns across all natural states, coping with major changes in neuronal activity without major departures from a single distribution of avalanche waiting times. Indeed, the results are compatible with the hypothesis that individual memories are encoded by specific spike avalanches, i.e. by stereotyped firing sequences within a given neuronal ensemble [10], [57]–[59]. Ongoing investigation of experience-dependent changes in avalanche repertoire shall clarify this issue.

Why is there a single regime of spike avalanches across all major behavioral states? A candidate common mechanism capable of unifying the dynamics of spike avalanches during natural behavior is the diffuse neuromodulatory drive from deep-brain centers crucially involved with attention, movement, motivation, sleep and arousal [60]–[64]. We propose that the severing of deep-brain neuromodulatory inputs by chemical (anesthetics) or physical methods (cell cultures, slices) abolishes long-range telencephalic coordination at high frequencies [65], [66], preventing the overlap of neuronal avalanches and disrupting the dynamic recruitment of distributed neuronal ensembles that characterizes behavior. Computer models show that the cooperative performance of neurons electrically connected by gap junctions is favored by critically-tuned coupling [11]. In this regard, the existence of gap junctions within deep-brain neuromodulatory centers with diffuse projections [67], [68] may provide a very apt mechanism to generate a single critical spiking regime throughout the telencephalon. Further experimentation is required to elucidate this hypothesis.

Materials and Methods

Ethics Statement

All animal work including housing, surgical and recording procedures were in strict accordance with the National Institutes of Health guidelines, and the Duke University Institutional Animal Care and Use Committee, and was approved by the Edmond and Lily Safra International Institute of Neuroscience of Natal Committee for Ethics in Animal Experimentation (permit # 04/2009).

MEA implants

A total of 14 adult male Long-Evans rats (300–350g) were used for electrophysiological recordings. Multielectrode Arrays (MEA; 35 µm tungsten wires, 16–32 wires per array, 250 or 500 µm spacing, 1 MΩ at 1 kHz) were surgically positioned within HP, S1 and V1 on the left hemisphere, according to the following stereotaxic coordinates in mm from Bregma with respect to the antero-posterior (AP), medio-lateral (ML) and dorso-ventral (DV) axes [69]: HP (AP: −2.80; ML: +1.5; DV: −2.80); S1 (AP: −3.00; ML: +5.5; DV: −1.40); V1 (AP: −7.30; ML: +4.00; DV: −1.30). DV measurements were taken with respect to the pial surface. Positioning was verified during or after surgery by spontaneous and evoked activity profiles, and confirmed by post-mortem histological analysis [70].

Neuronal recordings

One to five weeks after a 10-day recovery period, animals were recorded across the spontaneous wake-sleep cycle before and after object exposure (n = 7), or during anesthesia (n = 7). From each electrode spike times from up to 4 nearby neurons were sampled at 40KHz, whereas LFP were sampled at 500 Hz. Multiple action potentials (spikes) and local field potentials (LFP) were simultaneously recorded using a 96-channel Multi-Neuron Acquisition Processor (MAP, Plexon Inc, Dallas, TX), as previously described [23], [70]. Briefly, single-unit recordings were performed using a software package for real-time supervised spike sorting (see Fig. S1) (SortClient 2002, Plexon Inc, Dallas, TX). Spike sorting was based on waveform shape differences, peak-to-peak spike amplitudes plotted in principal component space, characteristic inter-spike-interval distributions, and a maximum 1% of spike collisions assuming a refractory period of 1 ms. Candidate spikes with signal-to-noise ratio lower than 2.5 were discarded. A waveform-tracking technique with periodic template adjustment was employed for the continuous recording of individual units over time. In order to ensure the stability of individual neurons throughout the experiment, waveform shape and single neuron clustering in principal component space were evaluated using graphical routines (WaveTracker software, Plexon, Dallas, TX). Ellipsoids were calculated by the cluster mean and 3 standard deviations corresponding to two-dimensional projections of the first and second principal components over consecutive 30 min data recordings. Strict superimposition of waveform ellipsoids indicated units that remained stable throughout the recording session and were therefore used for analyses, while units with nonstationary waveforms were discarded. Spike and LFP recordings were continuously performed before, during and after a 20 minutes experimental session in which animals engaged in the free exploration of four novel objects, as previously described [23], [70]. Neuronal ensembles of 45 to 126 neurons per rat were recorded for 4–6 hours. Visible lights were kept off throughout the experiment. For anesthetized recordings, animals received a single intramuscular administration of ketamine chlorhydrate (100 mg/kg) and xylazine (8 mg/kg), plus a subcutaneous injection of atropine sulfate (0.04 mg/kg) to prevent breathing problems. Anesthetized animals were placed inside a dark chamber and recorded for 4–6 hours, until they recovered waking behavior.

LFP-based classification of the major behavioral states

LFPs simultaneously recorded from S1, V1 and HP were used for the semi-automatic spectral classification of the three major behavioral states, WK, SWS and REM, as detailed in Ref. [21]. This method takes advantage of state-specific LFP power variations within different frequency bands as the sleep-wake cycle progresses, and has been successfully employed in the high throughput sorting of wake-sleep states in rodents [23], [70]–[72]. Briefly, two LFP amplitude ratios within specific spectral bands (0.5–20/0.5–55 Hz for ratio 1 and 0.5–4.5/0.5–9 Hz for ratio 2) were plotted in 2D principal component space to separate and sort data clusters corresponding to each of the three major behavioral states. In comparison with visual coding, this semi-automated method has >90% of accuracy, sensitivity and specificity [21].

Spike avalanche measurement

For each rat, let tij be the time of occurrence of the i-th spike of the j-th neuron. In order to define a neuronal avalanche, the spike time series {tij} (Fig. 1a) were divided in bins of duration Δt, as exemplified in Fig. 1b. The beginning of a neuronal avalanche is formally defined by the occurrence of a time bin without any spikes (in any neuron) followed by a time bin with at least one spike (in at least one neuron). The end of the avalanche is reached when another empty time bin occurs. The duration of the avalanche corresponds to the number of non-empty bins (times Δt), while its size is defined as the total number of spikes surrounded by empty bins [1].

Clearly, choosing larger (or smaller) values of Δt favors larger (or smaller) avalanche sizes. Results will therefore depend on the particular choice of time bin; for instance, all the avalanches collapse into a few large ones when Δt is large enough, whereas for very small Δt avalanches are split into smaller ones with a few spikes each. To rule out a systematic bias owing to the choice of time bin, we employed the same heuristic prescription as that of Ref. [1], namely to create a pooled time series with spikes from all neurons, and to use as time bin the average inter-event interval (IEI, see Fig. 1b), i.e. the time between consecutive spikes (whether or not from the same neuron). These rate-normalized time bins were therefore independently determined by the data, being specifically calculated for different rats, brain areas and behavioral states.

Surrogated data

Surrogated data were obtained by shuffling inter-spike intervals of each neuron independently, within single episodes of WK, SWS or REM states. This is a conservative procedure because these episodes are typically short (from seconds to minutes), and therefore the shuffling is limited. Furthermore, the inter-spike interval distribution of each neuron remains unchanged. Still, since the neurons are independently shuffled, across-neuron correlations are severely attenuated. Also note that the shuffling method employed does not change the average firing rate of each neuron. Therefore, the rate-normalized time bin for each scenario will be the same for original and surrogated data.

Cellular automaton model

We simulated a two-dimensional model where each site i (i = 1, …, L2, L = 32) is an excitable cellular automaton which cyclically goes through its four states (representing quiescence, excitation and two states for refractoriness). Quiescent neurons fire by excitation from firing nearest neighbors (with probability p per neighbor) or by external stimulus (with probability ph).

Above p = pc∼0.38 self-sustained activity becomes stable [11], [73]. We tune p to pc and employ a Poisson process to mimic stimuli independently arriving at electrode sites from the environment and from other brain regions. Choosing a small stimulus intensity ms−1, avalanches are continuously created, eventually colliding and/or overlapping.

For simulation of the AN group, the model was modified by periodically modulating the Poisson rate (with ms−1 and f = 2 Hz), which impinged only on 10% of the sites. The periodic modulation mimics the synchronization of spike bursts with LFPs, whereas stimulating only a fraction of the network mimics the reduction of synaptic input owing to the effect of ketamine and xylazine.

The spatial arrangement of the simulated neurons employed to measure the avalanches were the same as those of the experiment (inset of Fig. 3a). The ratio between electrode spacing and electrode measurement range matched the one estimated for the experiment. To prevent any bias deriving from the location of the simulated MEA in the network [24], periodic boundary conditions were used in simulations.

Power spectrum and DFA analysis

The power spectrum of the avalanche size time series was calculated with the Fast Fourier Transform [74]. Only continuous series longer than 1000 s were used to ensure a large enough number of avalanches, which restricted the analysis to WK states.

DFA analysis was performed following the standard procedures described in Ref. [27] and employing the software freely available at www.physionet.org/physiotools/dfa.

Supporting Information

[Figure omitted. See PDF.]

Figure S1.

Spike sorting and experimental design (a) The top panel shows the waveforms of two single units recorded from one electrode. The bottom panel shows that the two units can be separated as distinct clusters in a PC space. (b) The top panel shows the waveforms of multiple single units recorded from 16 channels. The bottom panel shows a rastergram of the sorted units. (c) Waveform stability was tracked throughout the experiment. Spike data (voltage-time ellipsoids, left panels) were sampled regularly from eight epochs of the total recording time (waveforms, right panels). The top left panels show good superposition of the ellipsoids, which indicates stability of a unit included in the study. The bottom left panels show discontinuity of the ellipsoids over time, indicating instability of a unit discarded from the study. (d) The FB animals were recorded across their spontaneous wakesleep cycle, comprising WK, SWS and REM. Recordings were performed before, during and after exposure to novel objects. This exposure consisted of a 20 minutes session in which four novel objects were placed inside the recording box (middle panel). Recordings made before (PRE, left panel) and after (POST, right panel) the exploration session lasted for up to 3h. Figure adapted from Ref. [23].

https://doi.org/10.1371/journal.pone.0014129.s001

(0.31 MB PDF)

Figure S2.

Size distributions obtained from a decreasing number of sampled neurons in the MEA are not qualitatively different. The plots show size distributions obtained from decreasing subsamples of the total number of neurons recorded in the experiment (FB1, WK, S1, PRE). Percentages indicate the sampled fraction of the recorded neurons. Inset: Size distributions obtained from the model.

https://doi.org/10.1371/journal.pone.0014129.s002

(0.10 MB PDF)

Figure S3.

Original vs surrogated FB avalanche size distributions. Comparison of the original (red) and surrogated (black) WK size distributions for different brain regions and stages of the experiment (in log-linear plots). Distributions were obtained by pooling avalanches from all FB rats.

https://doi.org/10.1371/journal.pone.0014129.s003

(0.13 MB PDF)

Figure S4.

Original vs surrogated AN size distributions. Comparison of the original (red) and surrogated (black) size distributions for some of the anesthetized animals (in log-linear plots).

https://doi.org/10.1371/journal.pone.0014129.s004

(0.10 MB PDF)

Figure S5.

Comparison of FB cumulative avalanche size distributions for different states and stages of the experiment. Cumulative distributions are shown together with the p values calculated from the KS tests. Note that the distributions are very similar in all cases, but only the ones in the left column pass the KS test.

https://doi.org/10.1371/journal.pone.0014129.s005

(0.11 MB PDF)

Figure S6.

Size distribution from AN data and the truncated power-law. (a) Avalanche size distribution for one AN rat and the truncated power-law fit (see text). (b) QQ-plot for the same AN data and fit. The solid line represents the linear fit which resulted in the slope a = 1.

https://doi.org/10.1371/journal.pone.0014129.s006

(0.11 MB PDF)

Table S1.

Time bin Δt, in milliseconds, calculated in each case. Since Δt is the mean interevent interval for a given condition, smaller sets of neurons typically lead to larger time bins (see also Table S2).

https://doi.org/10.1371/journal.pone.0014129.s007

(0.06 MB PDF)

Table S2.

Number of neurons sampled by the MEA per brain region for FB and AN rats.

https://doi.org/10.1371/journal.pone.0014129.s008

(0.04 MB PDF)

Text S1.

Statistical analysis of avalanche size distributions.

https://doi.org/10.1371/journal.pone.0014129.s009

(0.08 MB PDF)

Acknowledgments

We acknowledge fruitful discussions with O. Kinouchi, S. Coutinho and F. Cysneiros. We thank N. Vasconcelos for early help with data processing, J. Meloy, G. Lehew and G. Filho for manufacturing electrode arrays and stimulation devices, and A. Ragoni, M. Pacheco, L. Oliveira and S. Halkiotis for laboratory management.

Author Contributions

Conceived and designed the experiments: SR. Performed the experiments: FC HB SR. Analyzed the data: TLR MC. Contributed reagents/materials/analysis tools: TLR MC FC HB DRC MALN SR. Wrote the paper: TLR MC DRC MALN SR.

Citation: Ribeiro TL, Copelli M, Caixeta F, Belchior H, Chialvo DR, Nicolelis MAL, et al. (2010) Spike Avalanches Exhibit Universal Dynamics across the Sleep-Wake Cycle. PLoS ONE5(11): e14129. https://doi.org/10.1371/journal.pone.0014129

References

1. Beggs JM, Plenz D (2003) Neuronal avalanches in neocortical circuits. J Neurosci 23: 11167–11177.JM BeggsD. Plenz2003Neuronal avalanches in neocortical circuits.J Neurosci231116711177

2. Gireesh ED, Plenz D (2008) Neuronal avalanches organize as nested theta and beta/gamma-oscillations during development of cortical layer 2/3. Proc Natl Acad Sci U S A 105: 7576–7581.ED GireeshD. Plenz2008Neuronal avalanches organize as nested theta and beta/gamma-oscillations during development of cortical layer 2/3.Proc Natl Acad Sci U S A10575767581

3. Petermann T, Thiagarajan TC, Lebedev MA, Nicolelis MAL, Chialvo DR, et al. (2009) Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proc Natl Acad Sci U S A 106: 15921–15926.T. PetermannTC ThiagarajanMA LebedevMAL NicolelisDR Chialvo2009Spontaneous cortical activity in awake monkeys composed of neuronal avalanches.Proc Natl Acad Sci U S A1061592115926

4. Bak P (1997) How nature works: The science of self-organized criticality. Oxford: Oxford University Press. P. Bak1997How nature works: The science of self-organized criticalityOxfordOxford University Press

5. Chialvo DR (2004) Critical brain networks. Physica A 340: 756–765.DR Chialvo2004Critical brain networks.Physica A340756765

6. Langton CG (1990) Computation at the edge of chaos: phase transitions and emergent computation. Physica D 42: 12–37.CG Langton1990Computation at the edge of chaos: phase transitions and emergent computation.Physica D421237

7. Nykter M, Price ND, Aldana M, Ramsey SA, Kauffman SA, et al. (2008) Gene expression dynamics in the macrophage exhibit criticality. Proc Natl Acad Sci U S A 105: 1897–1900.M. NykterND PriceM. AldanaSA RamseySA Kauffman2008Gene expression dynamics in the macrophage exhibit criticality.Proc Natl Acad Sci U S A10518971900

8. Beekman M, Sumpter DJT, Ratnieks FLW (2001) Phase transition between disordered and ordered foraging in Pharaoh's ants. Proc Natl Acad Sci U S A 98: 9703–9706.M. BeekmanDJT SumpterFLW Ratnieks2001Phase transition between disordered and ordered foraging in Pharaoh's ants.Proc Natl Acad Sci U S A9897039706

9. Bertschinger N, Natschläger T (2004) Real-time computation at the edge of chaos in recurrent neural networks. Neural Comput 16: 1413–1436.N. BertschingerT. Natschläger2004Real-time computation at the edge of chaos in recurrent neural networks.Neural Comput1614131436

10. Beggs JM, Plenz D (2004) Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures. J Neurosci 24: 5216–5229.JM BeggsD. Plenz2004Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures.J Neurosci2452165229

11. Kinouchi O, Copelli M (2006) Optimal dynamical range of excitable networks at criticality. Nat Phys 2: 348–351.O. KinouchiM. Copelli2006Optimal dynamical range of excitable networks at criticality.Nat Phys2348351

12. Publio R, Oliveira RF, Roque AC (2009) A computational study on the role of gap junctions and rod Ih conductance in the enhancement of the dynamic range of the retina. PLoS One 4(9): e6970.R. PublioRF OliveiraAC Roque2009A computational study on the role of gap junctions and rod Ih conductance in the enhancement of the dynamic range of the retina.PLoS One49e6970

13. Shew W, Yang H, Petermann T, Roy R, Plenz D (2009) Neuronal avalanches imply maximum dynamic range in cortical networks at criticality. J Neurosci 29(49): 15595–15600.W. ShewH. YangT. PetermannR. RoyD. Plenz2009Neuronal avalanches imply maximum dynamic range in cortical networks at criticality.J Neurosci29491559515600

14. Eguíluz VM, Chialvo DR, Cecchi GA, Baliki M, Apkarian AV (2005) Scale-free functional brain networks. Phys Rev Lett 94: 018102.VM EguíluzDR ChialvoGA CecchiM. BalikiAV Apkarian2005Scale-free functional brain networks.Phys Rev Lett94018102

15. Sporns O, Chialvo DR, Kaiser M, Hilgetag CC (2004) Organization, development and function of complex brain networks. Trends Cog Sci 8: 418–425.O. SpornsDR ChialvoM. KaiserCC Hilgetag2004Organization, development and function of complex brain networks.Trends Cog Sci8418425

16. Kitzbichler MG, Smith ML, Christensen SR, Bullmore E (2009) Broadband criticality of human brain network synchronization. PLoS Comput Biol 5: e1000314.MG KitzbichlerML SmithSR ChristensenE. Bullmore2009Broadband criticality of human brain network synchronization.PLoS Comput Biol5e1000314

17. Linkenkaer-Hansen K, Nikouline VV, Palva JM, Ilmoniemi RJ (2001) Long-range temporal correlations and scaling behavior in human brain oscillations. J Neurosci 21: 1370–1377.K. Linkenkaer-HansenVV NikoulineJM PalvaRJ Ilmoniemi2001Long-range temporal correlations and scaling behavior in human brain oscillations.J Neurosci2113701377

18. Mazzoni A, Broccard FD, Garcia-Perez E, Bonifazi P, Ruaro ME, et al. (2007) On the dynamics of the spontaneous activity in neuronal networks. PLoS One 2(5): e439.A. MazzoniFD BroccardE. Garcia-PerezP. BonifaziME Ruaro2007On the dynamics of the spontaneous activity in neuronal networks.PLoS One25e439

19. Pasquale V, Massobrio P, Bologna LL, Chiappalone M, Martinoia S (2008) Self-organization and neuronal avalanches in networks of dissociated cortical neurons. Neurosci 153: 1354–1369.V. PasqualeP. MassobrioLL BolognaM. ChiappaloneS. Martinoia2008Self-organization and neuronal avalanches in networks of dissociated cortical neurons.Neurosci15313541369

20. Bédard C, Kröger H, Destexhe A (2006) Does the 1/f frequency scaling of brain signals reflect self-organized critical states? Phys Rev Lett 97: 118102.C. BédardH. KrögerA. Destexhe2006Does the 1/f frequency scaling of brain signals reflect self-organized critical states?Phys Rev Lett97118102

21. Gervasoni D, Lin SC, Ribeiro S, Soares ES, Pantoja J, et al. (2004) Global forebrain dynamics predict rat behavioral states and their transitions. J Neurosci 24: 11137–11147.D. GervasoniSC LinS. RibeiroES SoaresJ. Pantoja2004Global forebrain dynamics predict rat behavioral states and their transitions.J Neurosci241113711147

22. Paxinos G (2004) The rat nervous system. Elsevier Academic Press. G. Paxinos2004The rat nervous systemElsevier Academic Press

23. Ribeiro S, Shi X, Engelhard M, Zhou Y, Zhang H, et al. (2007) Novel experience induces persistent sleep-dependent plasticity in the cortex but not in the hippocampus. Frontiers in Neurosci 1: 43–55.S. RibeiroX. ShiM. EngelhardY. ZhouH. Zhang2007Novel experience induces persistent sleep-dependent plasticity in the cortex but not in the hippocampus.Frontiers in Neurosci14355

24. Priesemann V, Munk MHJ, Wibral M (2009) Subsampling effects in neuronal avalanche distributions recorded in vivo. BMC Neurosci 10: 40.V. PriesemannMHJ MunkM. Wibral2009Subsampling effects in neuronal avalanche distributions recorded in vivo.BMC Neurosci1040

25. Hennig MH, Adams C, Willshaw D, Sernagor E (2009) Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity. J Neurosci 29: 1077–1086.MH HennigC. AdamsD. WillshawE. Sernagor2009Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity.J Neurosci2910771086

26. Hahn G, Petermann T, Havenith MN, Yu S, Singer W, et al. (2010) Neuronal avalanches in spontaneous activity in vivo. J Neurophysiol. G. HahnT. PetermannMN HavenithS. YuW. Singer2010Neuronal avalanches in spontaneous activity in vivo.J NeurophysiolIn press (doi:10.1152/jn.00953.2009). In press (doi:10.1152/jn.00953.2009).

27. Peng CK, Havlin S, Stanley HE, Goldberger AL (1995) Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos 5(1): 82–87.CK PengS. HavlinHE StanleyAL Goldberger1995Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series.Chaos518287

28. Sanchez-Vives MV, McCormick DA (2000) Cellular and network mechanisms of rhythmic recurrent activity in neocortex. Nat Neurosci 2: 1027–1034.MV Sanchez-VivesDA McCormick2000Cellular and network mechanisms of rhythmic recurrent activity in neocortex.Nat Neurosci210271034

29. Pelt JV, Wolters PS, Corner MA, Rutten WLC, Ramakers GJA (2004) Long-term characterization of firing dynamics of spontaneous bursts in cultured cell networks. IEEE Trans Biomed Eng 51: 2051–2062.JV PeltPS WoltersMA CornerWLC RuttenGJA Ramakers2004Long-term characterization of firing dynamics of spontaneous bursts in cultured cell networks.IEEE Trans Biomed Eng5120512062

30. Steriade M, Nunez A, Amzica F (1993) A novel slow (<1 Hz) oscillation of neo-cortical neurons in vivo: Depolarizing and hyperpolarizing components. J Neurosci 13: 3252–3265.M. SteriadeA. NunezF. Amzica1993A novel slow (<1 Hz) oscillation of neo-cortical neurons in vivo: Depolarizing and hyperpolarizing components.J Neurosci1332523265

31. Amzica F, Steriade M (1998) Electrophysiological correlates of sleep delta waves. Clin Neurophysiol 107(2): 69–83.F. AmzicaM. Steriade1998Electrophysiological correlates of sleep delta waves.Clin Neurophysiol10726983

32. Greenberg DS, Houweling AR, Kerr JN (2008) Population imaging of ongoing neuronal activity in the visual cortex of awake rats. Nat Neurosci 11(7): 749–751.DS GreenbergAR HouwelingJN Kerr2008Population imaging of ongoing neuronal activity in the visual cortex of awake rats.Nat Neurosci117749751

33. Thomson AM, West DC, Lodge D (1985) An N-methylaspartate receptor mediated synapse in rat cerebral cortex: a site of action of ketamine? Nature 313: 479–481.AM ThomsonDC WestD. Lodge1985An N-methylaspartate receptor mediated synapse in rat cerebral cortex: a site of action of ketamine?Nature313479481

34. Harrison N, Simmonds M (1985) Quantitative studies on some antagonists of N-methyl D-aspartate in slices of rat cerebral cortex. Br J Pharmacol 84(2): 381–391.N. HarrisonM. Simmonds1985Quantitative studies on some antagonists of N-methyl D-aspartate in slices of rat cerebral cortex.Br J Pharmacol842381391

35. Leppävuori A, Putkonen PT (1980) Alpha-adrenoceptive influences on the control of the sleep-waking cycle in the cat. Brain Res 193(1): 95–115.A. LeppävuoriPT Putkonen1980Alpha-adrenoceptive influences on the control of the sleep-waking cycle in the cat.Brain Res193195115

36. Hartman BK, Zide D, Udenfriend S (1972) The use of dopamine-hydroxylase as a marker for the central noradrenergic nervous system in rat brain. Proc Natl Acad Sci USA 69(9): 2722–2726.BK HartmanD. ZideS. Udenfriend1972The use of dopamine-hydroxylase as a marker for the central noradrenergic nervous system in rat brain.Proc Natl Acad Sci USA69927222726

37. Hara K, Harris RA (2002) The anesthetic mechanism of urethane: the effects on neurotransmitter-gated ion channels. Anesth Analg 94: 313–318.K. HaraRA Harris2002The anesthetic mechanism of urethane: the effects on neurotransmitter-gated ion channels.Anesth Analg94313318

38. Vahle-Hinz C, Detsch O (2002) What can in vivo electrophysiology in animal models tell us about mechanisms of anaesthesia? Br J Anaesth 89(1): 123–142.C. Vahle-HinzO. Detsch2002What can in vivo electrophysiology in animal models tell us about mechanisms of anaesthesia?Br J Anaesth891123142

39. Jensen HJ (1998) Self-organized criticality: Emergent complex behavior in physical and biological systems. Cambridge: Cambridge University Press. HJ Jensen1998Self-organized criticality: Emergent complex behavior in physical and biological systemsCambridgeCambridge University Press

40. Baiesi M, Paczuski M, Stella AL (2006) Intensity thresholds and the statistics of the temporal occurrence of solar flares. Phys Rev Lett 96(5): 051103.M. BaiesiM. PaczuskiAL Stella2006Intensity thresholds and the statistics of the temporal occurrence of solar flares.Phys Rev Lett965051103

41. Davidsen J, Stanchits S, Dresen G (2007) Scaling and universality in rock fracture. Phys Rev Lett 98(12): 125502.J. DavidsenS. StanchitsG. Dresen2007Scaling and universality in rock fracture.Phys Rev Lett9812125502

42. Corral A, Telesca L, Lasaponara R (2008) Scaling and correlations in the dynamics of forest-fire occurrence. Phys Rev E 77: 016101.A. CorralL. TelescaR. Lasaponara2008Scaling and correlations in the dynamics of forest-fire occurrence.Phys Rev E77016101

43. Christensen K, Danon L, Scanlon T, Bak P (2002) Unified scaling law for earthquakes. Proc Natl Acad Sci U S A 99(1): 2509–2513.K. ChristensenL. DanonT. ScanlonP. Bak2002Unified scaling law for earthquakes.Proc Natl Acad Sci U S A99125092513

44. Corral A (2003) Local distributions and rate fluctuations in a unified scaling law for earthquakes. Phys Rev E 68: 035102.A. Corral2003Local distributions and rate fluctuations in a unified scaling law for earthquakes.Phys Rev E68035102

45. Omori F (1895) On the aftershocks of earthquakes. J Coll Sci Imper Univ Tokyo 7: 111–200.F. Omori1895On the aftershocks of earthquakes.J Coll Sci Imper Univ Tokyo7111200

46. Plenz D, Chialvo DR (2009) Scaling properties of neuronal avalanches are consistent with critical dynamics. D. PlenzDR Chialvo2009Scaling properties of neuronal avalanches are consistent with critical dynamics.arXiv:0912.5369v1. arXiv:0912.5369v1.

47. Churchland MM, Yu BM, Cunningham JP, Sugrue LP, Cohen MR, et al. (2010) Stimulus onset quenches neural variability: a widespread cortical phenomenon. Nat Neurosci 13(3): 369–378.MM ChurchlandBM YuJP CunninghamLP SugrueMR Cohen2010Stimulus onset quenches neural variability: a widespread cortical phenomenon.Nat Neurosci133369378

48. Trapani G, Latrofa A, Franco M, Altomare C, Sanna E, et al. (1998) Propofol analogues. Synthesis, relationships between structure and affinity at GABAA receptor in rat brain, and differential electrophysiological profile at recombinant human GABAA receptors. J Med Chem 41(11): 1846–1854.G. TrapaniA. LatrofaM. FrancoC. AltomareE. Sanna1998Propofol analogues. Synthesis, relationships between structure and affinity at GABAA receptor in rat brain, and differential electrophysiological profile at recombinant human GABAA receptors.J Med Chem411118461854

49. Krasowski MD, Jenkins A, Flood P, Kung AY, Hopfinger AJ, et al. (2001) General anesthetic potencies of a series of propofol analogs correlate with potency for potentiation of gamma-aminobutyric acid (GABA) current at the GABA(A) receptor but not with lipid solubility. J Pharmacol Exp Ther 297(1): 338–351.MD KrasowskiA. JenkinsP. FloodAY KungAJ Hopfinger2001General anesthetic potencies of a series of propofol analogs correlate with potency for potentiation of gamma-aminobutyric acid (GABA) current at the GABA(A) receptor but not with lipid solubility.J Pharmacol Exp Ther2971338351

50. Krasowski MD, Hong X, Hopfinger AJ, Harrison NL (2002) 4D-QSAR analysis of a set of propofol analogues: mapping binding sites for an anesthetic phenol on the GABA(A) receptor. J Med Chem 45(15): 3210–3221.MD KrasowskiX. HongAJ HopfingerNL Harrison20024D-QSAR analysis of a set of propofol analogues: mapping binding sites for an anesthetic phenol on the GABA(A) receptor.J Med Chem451532103221

51. Haeseler G, Leuwer M (2003) High-affinity block of voltage-operated rat IIA neuronal sodium channels by 2,6 di-tert-butylphenol, a propofol analogue. Eur J Anaesthesiol 20(3): 220–224.G. HaeselerM. Leuwer2003High-affinity block of voltage-operated rat IIA neuronal sodium channels by 2,6 di-tert-butylphenol, a propofol analogue.Eur J Anaesthesiol203220224

52. Bekinschtein TA, Shalom DE, Forcato C, Herrera M, Coleman MR, et al. (2009) Classical conditioning in the vegetative and minimally conscious state. Nat Neurosci 12(10): 1343–1349.TA BekinschteinDE ShalomC. ForcatoM. HerreraMR Coleman2009Classical conditioning in the vegetative and minimally conscious state.Nat Neurosci121013431349

53. Sánchez R, Newman DE, Carreras BA (2002) Waiting-time statistics of self-organized-criticality systems. Phys Rev Lett 88: 068302.R. SánchezDE NewmanBA Carreras2002Waiting-time statistics of self-organized-criticality systems.Phys Rev Lett88068302

54. Rothblat LA, Kromer LF (1991) Object recognition memory in the rat: the role of the hippocampus. Behav Brain Res 42(1): 25–32.LA RothblatLF Kromer1991Object recognition memory in the rat: the role of the hippocampus.Behav Brain Res4212532

55. Winters BD, Saksida LM, Bussey TJ (2008) Object recognition memory: neurobiological mechanisms of encoding, consolidation and retrieval. Neurosci Biobehav Rev 32(5): 1055–1070.BD WintersLM SaksidaTJ Bussey2008Object recognition memory: neurobiological mechanisms of encoding, consolidation and retrieval.Neurosci Biobehav Rev32510551070

56. Wallace CS, Withers GS, Weiler IJ, George M, Clayton DF, et al. (1995) Correspondence between sites of NGFI-A induction and sites of morphological plasticity following exposure to environmental complexity. Brain Res Mol Brain Res 32(2): 211–220.CS WallaceGS WithersIJ WeilerM. GeorgeDF Clayton1995Correspondence between sites of NGFI-A induction and sites of morphological plasticity following exposure to environmental complexity.Brain Res Mol Brain Res322211220

57. Nadasdy Z, Hirase H, Czurko A, Csicsvari J, Buzsaki G (1999) Replay and time compression of recurring spike sequences in the hippocampus. J Neurosci 19(21): 9497–9507.Z. NadasdyH. HiraseA. CzurkoJ. CsicsvariG. Buzsaki1999Replay and time compression of recurring spike sequences in the hippocampus.J Neurosci192194979507

58. Ikegaya Y, Aaron G, Cossart R, Aronov D, Lampel I, et al. (2004) Synfire chains and cortical songs: Temporal modules of cortical activity. Science 304: 559–564.Y. IkegayaG. AaronR. CossartD. AronovI. Lampel2004Synfire chains and cortical songs: Temporal modules of cortical activity.Science304559564

59. Ji D, Wilson MA (2007) Coordinated memory replay in the visual cortex and hippocampus during sleep. Nat Neurosci 10(1): 100–107.D. JiMA Wilson2007Coordinated memory replay in the visual cortex and hippocampus during sleep.Nat Neurosci101100107

60. Foote SL, Aston-Jones G, Bloom FE (1980) Impulse activity of locus coeruleus neurons in awake rats and monkeys is a function of sensory stimulation and arousal. Proc Natl Acad Sci U S A 77(5): 3033–3037.SL FooteG. Aston-JonesFE Bloom1980Impulse activity of locus coeruleus neurons in awake rats and monkeys is a function of sensory stimulation and arousal.Proc Natl Acad Sci U S A77530333037

61. Schultz W (1986) Activity of pars reticulata neurons of monkey substantia nigra in relation to motor, sensory, and complex events. J Neurophysiol 55(4): 660–677.W. Schultz1986Activity of pars reticulata neurons of monkey substantia nigra in relation to motor, sensory, and complex events.J Neurophysiol554660677

62. Horvath TL, Peyron C, Diano S, Ivanov A, Aston-Jones G, et al. (1999) Hypocretin (orexin) activation and synaptic innervation of the locus coeruleus noradrenergic system. J Comp Neurol 415(2): 145–159.TL HorvathC. PeyronS. DianoA. IvanovG. Aston-Jones1999Hypocretin (orexin) activation and synaptic innervation of the locus coeruleus noradrenergic system.J Comp Neurol4152145159

63. Lin SC, Nicolelis MAL (2008) Neuronal ensemble bursting in the basal forebrain encodes salience irrespective of valence. Neuron 59(1): 138–149.SC LinMAL Nicolelis2008Neuronal ensemble bursting in the basal forebrain encodes salience irrespective of valence.Neuron591138149

64. Ranade SP, Mainen ZF (2009) Transient firing of dorsal raphe neurons encodes diverse and specific sensory, motor and reward events. J Neurophysiol 102: 3026–3037.SP RanadeZF Mainen2009Transient firing of dorsal raphe neurons encodes diverse and specific sensory, motor and reward events.J Neurophysiol10230263037

65. Lima B, Singer W, Chen NH, Neuenschwander S (2010) Synchronization dynamics in response to plaid stimuli in monkey V1. Cereb Cortex 20(7): 1556–1573.B. LimaW. SingerNH ChenS. Neuenschwander2010Synchronization dynamics in response to plaid stimuli in monkey V1.Cereb Cortex20715561573

66. Canolty RT, Ganguly K, Kennerley SW, C. F. Cadieu CF, K. Koepsell K, et al. (2010) Oscillatory phase coupling coordinates anatomically dispersed functional cell assemblies. Proc Natl Acad Sci U S A. RT CanoltyK. GangulySW KennerleyCF C. F. CadieuK. K. Koepsell2010Oscillatory phase coupling coordinates anatomically dispersed functional cell assemblies.Proc Natl Acad Sci U S AIn press (doi:10.1073/pnas.1008306107). In press (doi:10.1073/pnas.1008306107).

67. Rash JE, Olson CO, Davidson KGV, Yasumura T, Kamasawa N, et al. (2007) Identification of connexin36 in gap junctions between neurons in rodent locus coeruleus. Neurosci 147(4): 938–956.JE RashCO OlsonKGV DavidsonT. YasumuraN. Kamasawa2007Identification of connexin36 in gap junctions between neurons in rodent locus coeruleus.Neurosci1474938956

68. Stezhka VV, Lovick TA (1995) Dye coupling between dorsal raphe neurones. Exp Brain Res 105: 383–390.VV StezhkaTA Lovick1995Dye coupling between dorsal raphe neurones.Exp Brain Res105383390

69. Paxinos G, Watson C (1997) The rat brain: In stereotaxic coordinates. San Diego: Academic Press. G. PaxinosC. Watson1997The rat brain: In stereotaxic coordinatesSan DiegoAcademic Press

70. Ribeiro S, Gervasoni D, Soares ES, Zhou Y, Lin SC, et al. (2004) Long lasting novelty-induced neuronal reverberation during slow-wave sleep in multiple forebrain areas. PLoS Biol 2: 126–137.S. RibeiroD. GervasoniES SoaresY. ZhouSC Lin2004Long lasting novelty-induced neuronal reverberation during slow-wave sleep in multiple forebrain areas.PLoS Biol2126137

71. Dzirasa K, Ribeiro S, Costa R, Santos LM, Lin SC, et al. (2006) Dopaminergic control of sleep-wake states. J Neurosci 26(41): 10577–10589.K. DzirasaS. RibeiroR. CostaLM SantosSC Lin2006Dopaminergic control of sleep-wake states.J Neurosci26411057710589

72. Pereira A, Ribeiro S, Wiest M, Moore LC, Pantoja J, et al. (2007) Processing of tactile information by the hippocampus. Proc Natl Acad Sci U S A 104: 18286–18291.A. PereiraS. RibeiroM. WiestLC MooreJ. Pantoja2007Processing of tactile information by the hippocampus.Proc Natl Acad Sci U S A1041828618291

73. Assis VRV, Copelli M (2008) Dynamic range of hypercubic stochastic excitable media. Phys Rev E 77: 011923.VRV AssisM. Copelli2008Dynamic range of hypercubic stochastic excitable media.Phys Rev E77011923

74. Press W, Teukolsky S, Vetterling W, Flannery B (1992) Numerical Recipes in C. Cambridge: Cambridge University Press. W. PressS. TeukolskyW. VetterlingB. Flannery1992Numerical Recipes in CCambridgeCambridge University Press