Dynamical criticality during induction of anesthesia in human ECoG recordings

1. Introduction

It has been suggested that neural systems operate in a critical regime similar to phase transitions in physics, given several computational desirable features of such states represented by the statistics of the thermodynamic variables (Chris, 1990). Evidence for statistical criticality is based on the observation that various aspects of neuronal activity such as avalanches observed in local field potentials and action potentials in tissue preparations and in animal models (Gireesh and Plenz, 2008; Ribeiro et al., 2010), as well as magneto-encephalography (MEG) and electro-corticography (ECoG) in human subjects (He et al., 2010; Shriki et al., 2013), exhibit long tailed-distributions well approximated by power laws. The critical regime provides important functional benefits; quantities such as dynamic range and information transmission are optimized near criticality (Shew and Plenz, 2013).

More recently, the dynamical aspect of criticality has been brought into focus, as a similarly desirable feature not fully captured by steady-state statistics such as avalanche size distributions (Magnasco et al., 2009; Chialvo, 2010; Mora and Bialek, 2011; Beggs and Timme, 2012); a perturbation in an extended dynamical system that is close to a critical point will neither decay nor explode, thus allowing for long range communication across the entire system. In contrast, if the system is far from criticality (therefore stable), perturbations damp out and no information integration takes place beyond the characteristic damping time scale (Tononi, 2008).

While models of self-organized criticality exhibit both dynamically and statistically critical behavior (Bak et al., 1987; Gil and Sornette, 1996), the two aspects of criticality are not necessarily related. The winnerless network provides an illuminating example: under very generic conditions, neural systems can display a phase space determined by heteroclynic orbits connecting saddle nodes (i.e., at least one unstable manifold), such that the resulting dynamics are quasi-periodic cycles over the nodes, without necessarily exhibiting statistically critical distributions (Rabinovich et al., 2001; Aguiar et al., 2011; Ashwin et al., 2011). A model connecting statistical and dynamical criticality in neural systems was proposed recently by Magnasco et al. (2009). They consider an abstract model in which the activity of a set of neurons is encoded in a N-dimensional vector $\overrightarrow{x}$ which evolves in time according to a N × N connectivity matrix A, characterized by its set of N eigenvalues {λ_n}. By assuming anti-Hebbian dynamics for the connectivity matrix a very rich dynamical scenario emerges. The eigenvalues of the matrix A evolve toward the dynamically critical point Re(λn) ≈ 0 ∀ n and the solutions of the model exhibit complex spatio-temporal dynamics, as well as long tailed avalanche distributions and other signatures of statistical criticality. Consistent with this observation, experimental evidence of both statistical and dynamical criticality was reported in human ECoG recordings; however, the precise mechanism by which critical dynamics occur has not been investigated. The analysis showed that the eigenvalues crowd near the critical line, and moreover that task performance (finger tapping) implies a subtle but significant decrease in dynamical criticality, presumably because the modes related to motor execution impose higher stability (Solovey et al., 2012). Of note, signatures of statistical criticality were not strongly affected by task performance.

If indeed dynamical criticality is a useful feature of brain activity rather than an epiphenomenon, stability of neuronal dynamics ought to be modulated by the behavioral state of the subject. Here, we hypothesized that a particularly dramatic change in stability accompanies changes in the level of wakefulness (consciousness). When the brain is awake and displaying complex behavior its dynamical state ought to be close to a bifurcation point; marginally stable modes contribute to long range interactions across the system. Conversely when higher-order functions associated with wakefulness have been diminished and eventually completely disrupted by anesthetics, brain dynamics should exhibit more stability. In other words, anesthesia induction should lead to stabilization of brain dynamics.

Changes in the level of arousal (wakefulness) have been historically quantified using spectral analysis of neuronal activity. In this view, decrease in the level of wakefulness is reflected in the increase and prevalence of low frequency oscillations and the concurrent decrease in the high frequency oscillations reviewed in Brown et al. (2010). While this is true for some states of decreased arousal such as slow wave sleep, this association breaks down during other states in which arousal is similarly depressed such as rapid eye movement (REM) sleep for instance. Furthermore, state of general anesthesia can be characterized by different spectral signatures depending on the specific choice of anesthetic agent (Maksimow et al., 2006). This makes current modes of detecting the “depth of anesthesia” unreliable (Avidan et al., 2011).

Lack of clear association between changes in the spectral content of brain signals and level of arousal is not entirely surprising. It is likely that the overall level of wakefulness is a consequence of the interactions among many brain regions rather than any specific feature of neuronal activity observed at any one region taken in isolation. Therefore, more recent efforts have been aimed at detecting decreases in arousal using connectivity measures based on spectral coherence as well as mutual information and phase relationships among brain activity recorded simultaneously at multiple locations (Imas et al., 2005; Cimenser et al., 2011; Lee et al., 2012). While this connectivity analysis does suggest that integration of information between different brain regions may be decreased when the level of wakefulness is reduced, it is not trivial to relate changes in connectivity to the changes in global dynamics of the brain.

To address the dynamics, we fitted vector autoregressive (VAR) models to ECoG signals collected directly from the cortex of human subjects as they were gradually induced into the state of general anesthesia. These models were independently fit to short temporal windows with an arbitrarily large overlap. Thus, while we assume that the dynamics are locally linear and stationary over a short temporal window, on a longer time scale the dynamics are expected to be arbitrarily non-linear and non-stationary. This locally linear approximation allows us to quantify the changes in stability of brain activity in terms of temporal evolution of the distribution of eigenmodes of the fitted models. As previously reported (Solovey et al., 2012), we found a prevalence of critical eigenmodes across the entire recordings. However, the stability of the models shows statistically significant differences across different stages of induction. While the distribution of eigenvalues changes in non-trivial ways, high frequency modes become more damped as anesthesia is induced. Moreover, modes closer to criticality, regardless of frequency, show a gradual shift to stability spanning several drug volleys over approximately 20 min.

This work is organized as follows. In the next section we describe the induction protocol and the analysis method. We present our results in section 3. In section 4 we summarize and discuss our findings.

2. Methods

All experimental protocols were approved by the IRB at the Weill Cornell Medical College (protocol number 1106011763). After obtaining informed consent, three subjects undergoing surgical treatment for intractable epilepsy were enrolled in this study. Subdural electrode grids and strips (Ad-tech, Medical Instruments Corp., Racine, WI) were implanted for the purposes of localization of the epileptogenic loci. The location and the number of electrodes were determined by the clinical considerations (temporal lobe for all subjects in this study). After the initial implantation of the subdural electrodes, the subjects underwent video and EEG monitoring, duration of which was dictated solely by clinical considerations (1–2 weeks in these subjects). The recordings analyzed in this work were obtained during induction of anesthesia for the second craniotomy performed after completion of this observation period. During induction of anesthesia (see below), blood pressure, ECG, heart rate, pulse oxymetry, and end tidal carbon dioxide were monitored and maintained within normal limits. Patients were given supplemental oxygen via nasal cannula.

After obtaining baseline recordings (without any pre-medication) anesthesia was gradually induced using target controlled infusions of propofol using pharmacokinetic parameters derived by Schnider et al. (1999), administered using STANPUMP. Target propofol concentration was increased slowly while the level of sedation was accessed using responses to simple verbal commands. Propofol infusion continued until subjects lost the ability to respond to verbal commands. At this point additional propofol, opioids, and neuromuscular blockers were administered (at the discretion of the anesthesia provider) and trachea was intubated. Recordings were terminated at this point.

Recordings were obtained using SynAmps² (Neuroscan) using DC coupled recording. Data were acquired at 10 KHz. 64 channels of ECoG signals were acquired in each subject. While both conventional EEG and ECoG are thought to primarily reflect the sum of synchronized postsynaptic potentials of neurons in the vicinity of the electrode, the invasive nature of the ECoG signals allows for much greater signal to noise ratio and significantly improves spatial and temporal resolution of the signals.

No online filtering was performed. ECoG data was collected from three human subjects as they were induced into general anesthesia. For all subjects, the infusion started 60 s into the recording and the concentration of anesthetics was increased every 300 s. For Subject 1, propofol infusion started 60 s into the recording. 360 s into the recording the subject reports being awake. 510 s into the recording the subject no longer responds. 960 s into the recording the subject is given additional drugs and intubated. For Subject 2, propofol was incremented at 300, 600, and 900 s. At 660 s the subject opened eyes. 720 s into the recording the subject no longer responded. 1140 s into the recording subject was given additional drugs and was intubated. For Subject 3, propofol infusion started 60 s into the recording. The concentration was increased every 300 s and maintained constant before and after. 900 s into the recording the subject no longer responded to verbal commands or light taps on the shoulder. 1200 s into the recording subject was given additional drug and was intubated.

Data was bandpass filtered at 0.1 − 500 Hz and detrended in segments of 10 s. We applied notch filters at 60 ± 2, 120 ± 2, and 180 ± 2 Hz. Finally, the amplitude of the signal in each channel was normalized by its standard deviation. For our analysis data was partitioned in equally sized windows of τ = 200 ms centered every t_{j + 1} − t_j = Δ = 100 ms. In each window, we assumed that the dynamics is locally linear and fitted a vector auto regressive model (VAR) of order p = 1.

$\begin{array}{cc}{y}_{n+1}=A{y}_n+u_n& \left(1\right)\end{array}$

where y_n are the fitted values, A ∈ ℝ^{N × N} is the matrix to be estimated and u_n is assumed to be white noise. Here, y_n ∈ ℝ^N is a multivariate time series that represents the recorded activity in all channels at time t_n and A corresponds to the LAG 1 correlation between channels. A comprehensive treatment of this model and its estimation can be found in Lütkepohl (2006). In this work we used a python implementation of Schnider's et al. algorithm to estimate A (Schneider and Arnold, 2001). This procedure yields a set of matrices A_j which govern the stability properties of the VAR model at time t_j. In order to address changes in the stability of the fitted models we considered the distribution of the modulus of the eigenvalues at each time step. Also, since our underlying hypothesis corresponds to a continuum model we performed a transformation in order to obtain a correspondence between the eigenvalues of A_j and the timescales of the dynamics. Let λ_j = ρ_j e^iϕ be the eigenvalue corresponding to the j-th mode, the frequency of the mode is given by

$\begin{array}{cc}{f}_j=\frac{{\varphi }_j}{2\pi dt}& \left(2\right)\end{array}$

while the growth rate (timescale) of the mode is given by

$\begin{array}{cc}{\tau }_j=\frac{log({\rho }_j)}{dt}& \left(3\right)\end{array}$

Here dt = $\frac{1}{S_f}$ = 0.0001 s, where S_f is the sampling frequency of the recordings. A mode with eigenvalue λ is critical if

$\begin{array}{cc}\Vert \lambda \Vert =1& \left(4\right)\end{array}$

In practice however, we call a mode critical if ǁλǁ ≈ 1 (thus τ ≈ 0 s). These are modes which are close to alternate their behavior between damping and growth (Strogatz, 2006).

The distributions so obtained were compared to the initial distribution (prior to induction) by means of two statistical tests. Kolmogorov-Smirnov (KS) tests the null hypothesis that the distributions are the same and yields the maximal difference of the cumulative distributions to quantify for the changes. Wilcoxon rank-sum (W) tests the null hypothesis that the distributions are the same against the alternative hypothesis that they are shifted and returns a z-value to account for the magnitude of the shift. If the values of the subsequent distributions are smaller than those of the reference distribution (awake state) then z > 0, therefore, an increase of the z-value indicates an increase of the stability.

We settled on a VAR-1 model because the main results related to the effect of anesthesia are robust for VAR-2 and VAR-3 models. We have explored window sizes ranging from 100 ms to 1 s and found no significant changes. Our method was tested against surrogate data obtained by phase randomization of the signal; for each channel we computed the Fourier transform of the signal, changed the phase value by a random number [drawn from a flat distribution in (0, 2π)] and transformed back to obtain the surrogated signals. Note that by construction this procedure preserves the power spectrum of each signal.

3. Results

We performed VAR analysis on three human subjects as they were induced into general anesthesia. Our primary focus was to detect changes in the distribution of the stability parameters ǁλ_jǁ during induction of anesthesia. To quantify changes in the stability of the models we used two non-parametric statistical tests [Kolmogorov-Smirnov (KS) and Wilcoxon rank-sums(W)]. The results of this analysis are shown in Figure 1A [for each subject top row shows (KS) and bottom row shows (W)]. To improve visualization the results were smoothed using moving average windows of 10 s. The distribution of eigenmodes computed over different windows during the awake state fluctuates. To scale the observed differences in stability during induction of anesthesia by these spontaneous fluctuations, we computed the time average of both KS and W statistics over the awake period (1 min) and subtracted this value from the curves shown in Figure 1A. In all cases, the temporal average of the p-values behaves similarly to the KS-Z values. During the first minutes of the procedure we find that p ≈ 0.75, thus, the null hypothesis that the distributions are the same cannot be safely rejected. However, we find a drastic drop of the p-value concomitant with changes in KS-Z values. For the regions indicated in blue and green (Figure 1A), the average p-value of both tests are in the range of 0.2−0.3 suggesting that the distributions have changed. While the KS test simply indicates that the distributions of stability parameters during awake and anesthetized states are different, the increase in the z-values of the Wilcoxon test implies that ǁλ_jǁ tends to decrease with induction of anesthesia, i.e., the dynamics is becoming more stable. Note that the change in the distribution of the stability parameter is not observed in phase randomized surrogates (red curves in Figure 1A). Thus, the observed changes in stability are not given by the spectral properties of neuronal activity.

FIGURE 1

Figure 1. ECoG signals were recorded from three human subjects as they were induced into general anesthesia. Data was locally fitted with VAR(1) models in windows of 200 ms every 100 ms (see methods). The linear stability of each model is compared to the awake state by means of two statistical tests. (A) top rows: Kolmogorov Smirnov test. For each model we plot the KS statistics of comparing the fitted distribution of eigenvalues against the awake state. (A) bottom rows: The distribution of time scales is compared using Wilcoxon test. Both quantities were averaged in time intervals of 10 s. The stability properties of locally fitted VAR(1) models change as the subjects undergo anesthesia. We defined three different segments (color rectangles) which were used for subsequent figures. (B) Changes in the frequency and stability of the eigenmodes. We compared the distributions of frequencies and time scales using a Wilcoxon test. In each figure, the vertical axis shows the z-value of comparing the frequency distributions whereas the horizontal axis shows the same test for the stability parameters distributions. The color code represents time elapsed since the beginning of the recording. In this representation all realizations yield qualitatively similar results: as the subjects are induced, the fitted frequencies shift to higher values at the same time they become more damped.

Note that in general the eigenvalues of the autoregressive matrices fitted to the ECoG signals are complex numbers whose real and imaginary parts give rise to the timescale τ and frequency f of the corresponding eigenmode (see Equations 2, 3). While Figure 1A focused just on changes in the distribution of the stability parameters, Figure 1B shows changes in both the distribution of timescales (abscissa) and bulk frequencies (ordinate) treated independently. Time elapsed since the onset of experiment is color coded from red (awake) to blue (anesthetized). The bulk evolution of the eigenmodes is consistent in all subjects: as induction progresses, modes shift to higher frequencies while they become more stable. To validate that the results obtained with the VAR-1 model are robust, we show in Figure A1 (included as Appendix) the same analysis as in Figure 1A implemented with a VAR-3 model. As it can be seen, the changes in the distribution of eigenmodes are almost identical to those for VAR-1.

While Figure 1B suggests an increase in the bulk frequency and decrease in the time constant, this does not fully characterize the way in which anesthetics change the distribution of eigenvalues in the plane spanned by timescale and frequency. Figure 2 shows how we represent the distributions of the eigenvalues of A_j. The vertical axes corresponds to frequencies plotted on a logarithmic (base 2) scale. Horizontal axes indicate the modes damping/growth timescale. The sign indicates whether the mode's amplitude is growing (positive) or decaying (negative). Histograms are color coded with blue indicating low occupancy to red indicating high occupancy. Note that the damping time and frequency are not independent and modes with lower frequencies tend to have longer damping times, with slow oscillations found near the critical point (τ ≈ 0). Traces on the margin of the figure illustrate the dynamics for particular pairs of damping time and frequency. Note that the traces are plotted on the timescale commensurate to the damping time rather than on an absolute time scale. The inter-relationship between damping time and frequency assures that most modes located along the most densely populated ridge go through several complete cycles before being damped out, while the modes located to the left of the dominant ridge are damped out earlier and are thus less likely to carry out meaningful computations performed by the brain.

FIGURE 2

Figure 2. Qualitative behavior of eigenmodes. The histogram corresponds to the eigenmodes of VAR-1 processes fitted to ECoG signals. The count in each bin is color-coded and the number of samples is N > 10⁶. The frequency axis is in logarithmic scale (base-2). The arrows indicate points in the stability plane for which the qualitative dynamics of the corresponding mode is illustrated. The dynamics of each mode can be expressed as an oscillation of frequency f whose amplitude (red curves) is modulated by an exponential decay/growth (blue curve). Each solution is shown for 2/τ s. Note that for the points in the plane with non-zero count, a number of oscillations occur before the mode is damped out. For the case labeled with ⋆, the mode grows exponentially (i.e., it is super-critical).

Figure 3A shows the distribution of eigenvalues in the plane introduced in Figure 2 during three stages of the induction process (100 s segments shown in Figure 1A). In order to better resolve the distributions we performed a moving VAR analysis with t_{j + 1} − t_j = 1 ms of spacing between adjacent windows. In order to visualize changes in the eigenvalue distributions we normalized the count value of each histogram by its maximum. Then, we used the normalized values in each bin to code for color in RGB space as indicated in the filled circles. Figure 3B, correspond to the superposition of such images. In this way, regions of the stability space that are similarly occupied in the three stages are coded in gray scale [with white corresponding to maximal occupancy (1,1,1)] and pure colors RGB correspond to values that are exclusive to the first, second and third stage respectively. A prominent feature shown by these panels is the shift of high-frequency eigenmodes toward increased stability. While the full worm-like distribution of eigenvalues changes in subtle ways, the left-ward shift in these frequencies is ubiquitous in all subjects. Figures 4A,B correspond to vertical and horizontal “slices” respectively of the histogram shown in Figure 3A for subject 1. Figure 4C shows details of how these distributions change for subject 1. We performed the same comparison as before but choosing a smaller frequency range for computing the histograms. The shift to damped states is more pronounced for modes with frequencies that are greater than 64 Hz.

FIGURE 3

Figure 3. Statistics of the models stability at each stage indicated in Figure 1. Each stage is sampled in 100 s intervals and VAR(1) models were fitted every 1 ms. For each stage we computed the frequency and growth \decay timescale of each mode. Figures correspond to 2D histograms of this quantities. The count in each bin is color coded and the number of samples is N ≈ 10⁶. The frequency axis is on logarithmic (base 2) scale. (A) Distributions of eigenvalues. (B) Differences across segments. We normalized each bin in (A) histograms by their maximum. The figure is constructed by superposing the three histograms, each coding for a color in RGB space. In this way pure red, pure green and pure blue correspond to eigenvalues that are only present in the first, second and third stage. The rainbow-like pattern indicates a shift of high frequency modes as they become more damped.

FIGURE 4

Figure 4. Detailed view of eigenmodes distributions for subject 1. Colors indicate the segments indicated in Figure 1A. (A) Histogram in Figure 3A is sliced by restricting the analysis to frequency bands of 5 Hz centered around the indicated values. Note that the blue histograms are always to the left of the red ones indicating increased stability as anesthesia is induced. (B) Same analysis as (A) restricted to slices of 50 $\frac{1}{s}$. Note the emergence of highly damped high frequency oscillations in the anesthetized condition (blue). In both cases we plot count number on the scale (0.46000) (A) and (0.120000) for (B). (C) Similarly to Figure 3B, we computed the superposition of histograms in a smaller frequency range for better visualization. The histograms were done with logarithmic binning and the frequency axis is logarithmic (base 2). The rainbow indicates a shift toward more damped states. The organization of stability undergoes non-trivial changes in low frequency bands (4–128 Hz).

Finally, we investigated how the distribution of the most critical modes is affected by induction. This was partially inspired by results previously reported in human ECoG, showing that differences between task and resting conditions can be detected precisely by changes in these populations (Solovey et al., 2012). We show in Figure 5 the result of comparing the distribution of modes truncated to eigenvalues with damping constant above a given threshold close to criticality. For all subjects, the distributions show a gradual change in the stability of near-critical modes along the entire span of the induction process Figure 5A. This is somewhat surprising, as the induction process is controlled by discrete events of drug increase which notably affect the full eigenmode distribution.

FIGURE 5

Figure 5. Stability analysis for truncated distributions of eigenvalues. We performed the same analysis as in Figure 1A for truncated distributions of eigenvalues; we kept all the eigenvalues with damping constant greater than a threshold. Each color corresponds to different threshold as indicated in the figure labels. (A) KS statistics of comparing the truncated distribution of eigenvalues. (B) Same as (A) using Wilcoxon z-value statistics.

4. Conclusions

Dynamical systems theory indicates that systems that are capable of performing computations should have a large number of modes with marginal stability. In such a scenario an arbitrary perturbation will not decay or explode, thus allowing for information integration across the entire system. Previous work suggest that the brain might operate in a dynamically critical regime (Magnasco et al., 2009; Solovey et al., 2012). A simple model exhibiting complex spatio-temporal dynamics was recently proposed, in which statistically critical behavior emerges due to dynamical instabilities. Within this framework we tested the hypothesis that the stability properties of the system change as anesthesia is induced; specifically, we hypothesized that wakefulness is related to dynamical criticality while the anesthetized state corresponds to increased damping of the dynamics. To test this hypothesis we assumed locally linear dynamics estimated in short segments of the recordings using eigenmode decomposition of VAR models.

We found that as the subjects become anesthetized the linear stability of the ECoG recordings show significant changes which are efficiently tracked by non-parametric statistical methods. These markers are remarkably robust to changes in the way data is normalized (choice of filters, amplitude normalization, sampling frequency). Moreover, changes in this quantities were found to be consistent with the subjects behavior as reported by the medical team. This suggests that our indicators could be used to monitor depth of anesthesia.

Our results are also consistent with the criticality hypothesis: we found a prevalence of modes close to criticality across the whole induction procedure. However, as the subjects became anesthetized there were significant changes in the stability properties of the fitted dynamics. These changes were examined closely in selected stages of the procedure and are visualized by the superposed histograms in Figures 3B, 4C. This analysis revealed that changes in the stability exhibit much richer structure than a simple block shift to damping across all frequencies. Yet, we observe a consistent pattern in all three subjects; the eigenvalues of the fitted models shift toward higher frequencies and increased damping. This should be interpreted carefully; it is not necessarily the case that there is an increase of high frequency spectral content of the ECoG signals. Although there ought to be a relationship between a moving spectral analysis and the eigenmodes of a moving VAR analysis, this relationship may be complex.

The increase in the prevalence of eigenmodes characterized by high frequency (high gamma) may be seen as surprising given the well-known observation that the power of high frequency oscillations tends to decrease with some anesthetics including propofol. This result, however, ought to be interpreted carefully. The increase in the number of eigenmodes does not equate to the increase in power. For instance, there could be fractionation of a single correlated pattern of high frequency oscillations in the awake state into multiple mutually independent patterns of high frequency oscillations.

The finding that high frequency modes become more damped as the subject is anesthetized is to some extent reassuring. If we adopt the traditional view that high frequency activity is associated to cognitive processes our results are consistent with an appealing interpretation. The effect of the anesthetic procedure is to damp out high frequency activity while still allowing for low frequency modes to perform a function. Low frequency activity can then presumably be associated to the maintenance tasks which keep the subject alive.

A number of recent reports have been aimed at characterizing criticality as a universal feature in ECoG recordings (He et al., 2010), and as particularly relevant to differentiate wakefulness from sleep (Meisel et al., 2013; Tagliazucchi et al., 2013) (see also Ribeiro et al., 2010 for comparable results with action potential recordings). In this context, our results provide support for a consistent and theoretically founded interpretation of the relationship between criticality and wakefulness. While the theoretical model is not the focus of the present publication, it is interesting to note that it implies a specific and falsifiable prediction: the model achieves self-tuned criticality by means of plastic synaptic adaptation. It follows that blocking synaptic changes should result in a breakdown of criticality; similarly, the model should also be able to explain changes in criticality during the sleep cycle, given the concomitant changes in plasticity patterns (Ribeiro et al., 2007). This will be the subject of future publications, along with further validation of the stabilizing effect of anesthesia in animal models, effects of different anesthetic agents, larger number of subjects, recovery from anesthesia, and application of the methods to EEG recordings.

5. Author Contributions

Leandro M. Alonso, performed data analysis; Leandro M. Alonso, Alex Proekt, Guillermo A. Cecchi, and Marcelo O. Magnasco wrote the manuscript. Alex Proekt, Guillermo A. Cecchi, Marcelo O. Magnasco designed the experiments; Alex Proekt, Kane O. Pryor, and Theodore H. Schwartz performed the experiments; Leandro M. Alonso and Alex Proekt contributed equally to this work.

Conflict of Interest Statement

The Guest Associate Editor A. Ravishankar Rao declares that, despite being affiliated to the same institution as the author Guillermo A. Cecchi, the review process was handled objectively and no conflict of interest exists. The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was funded by NSF grant EF-0928723 awarded to Marcelo O. Magnasco and by NIGMS (1K08GM106144-01) awarded to Alex Proekt

References

Aguiar, M., Ashwin, P., Dias, A., and Field, M. (2011). Dynamics of coupled cell networks: synchrony, heteroclinic cycles and inflation. J. Nonlin. Sci. 21, 271–323. doi: 10.1007/s00332-010-9083-9

CrossRef Full Text

Ashwin, P., Karabacak, O., and Nowotny, T. (2011). Criteria for robustness of heteroclinic cycles in neural microcircuits. J. Math. Neurosci. 1, 1–18. doi: 10.1186/2190-8567-1-13

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Avidan, M. S., Jacobsohn, E., Glick, D., Burnside, B. A., Zhang, L., Villafranca, A., et al. (2011). Prevention of intraoperative awareness in a high-risk surgical population. N. Engl. J. Med. 365, 591–600. doi: 10.1056/NEJMoa1100403

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bak, P., Tang, C., and Wiesenfeld, K. (1987). Self-organized criticality: an explanation of 1/f noise. Phys. Rev. Lett. 59, 381–384. doi: 10.1103/PhysRevLett.59.381

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Beggs, J. M., and Timme, N. (2012). Being critical of criticality in the brain. Front. Physiol. 3:163. doi: 10.3389/fphys.2012.00163

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Brown, E. N., Lydic, R., and Schiff, N. D. (2010). General anesthesia, sleep, and coma. N. Engl. J. Med. 363, 2638–2650. doi: 10.1056/NEJMra0808281

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Chialvo, D. R. (2010). Emergent complex neural dynamics. Nat. Phys. 6, 744–750. doi: 10.1038/nphys1803

CrossRef Full Text

Chris G. L. (1990). Computation at the edge of chaos: phase transitions and emergent computation. Phys. D 42, 12–37. doi: 10.1016/0167-2789(90)90064-V

CrossRef Full Text

Cimenser, A., Purdon, P. L., Pierce, E. T., Walsh, J. L., Salazar-Gomez, A. F., Harrell, P. G., et al. (2011). Tracking brain states under general anesthesia by using global coherence analysis. Proc. Natl. Acad. Sci. U.S.A. 108, 8832–8837. doi: 10.1073/pnas.1017041108

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gil, L., and Sornette, D. (1996). Landau-Ginzburg theory of self-organized criticality. Phys. Rev. Lett. 76, 3991. doi: 10.1103/PhysRevLett.76.3991

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gireesh, E. D., and 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. doi: 10.1073/pnas.0800537105

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

He, B. J., Zempel, J. M., Snyder, A. Z., and Raichle, M. E. (2010). The temporal structures and functional significance of scale-free brain activity. Neuron 66, 353–369. doi: 10.1016/j.neuron.2010.04.020

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Imas, O. A., Ropella, K. M., Ward, B. D., Wood, J. D., and Hudetz, A. G. (2005). Volatile anesthetics disrupt frontal-posterior recurrent information transfer at gamma frequencies in rat. Neurosci. Lett. 387, 145–150. doi: 10.1016/j.neulet.2005.06.018

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lee, U., Lee, H., Müller, M., Noh, G. J., and Mashour, G. A. (2012). Genuine and spurious phase synchronization strengths during consciousness and general anesthesia. PLoS ONE 7:e46313. doi: 10.1371/journal.pone.0046313

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lütkepohl, H. (2006). New Introduction to Multiple Time Series Analysis. Berlin; Heidelberg: Springer-Verlag.

Magnasco, M. O., Piro, O., and Cecchi, G. A. (2009). Self-tuned critical anti-Hebbian networks. Phys. Rev. Lett. 102, 258102. doi: 10.1103/PhysRevLett.102.258102

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Maksimow, A., Särkelä, M., Långsjö, J. W., Salmi, E., Kaisti, K. K., Yli-Hankala, A., et al. (2006). Increase in high frequency EEG activity explains the poor performance of EEG spectral entropy monitor during S-ketamine anesthesia. Clin. Neurophysiol. 117, 1660–1668. doi: 10.1016/j.clinph.2006.05.011

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Meisel, C., Olbrich, E., Shriki, O., and Achermann, P. (2013). Fading signatures of critical brain dynamics during sustained wakefulness in humans. J. Neurosci. 33, 17363–17372. doi: 10.1523/JNEUROSCI.1516-13.2013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mora, T., and Bialek, W. (2011). Are biological systems poised at criticality? J. Stat. Phys. 144, 268–302. doi: 10.1007/s10955-011-0229-4

CrossRef Full Text

Rabinovich, M., Volkovskii, A., Lecanda, P., Huerta, R., Abarbanel, H. D. I., and Laurent, G. (2001). Dynamical encoding by networks of competing neuron groups: winnerless competition. Phys. Rev. Lett. 87:068102. doi: 10.1103/PhysRevLett.87.068102

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ribeiro, S., Shi, X., Engelhard, M., Zhou, Y., Zhang, H., Gervasoni, D., et al. (2007). Novel experience induces persistent sleep-dependent plasticity in the cortex but not in the hippocampus. Front. Neurosci. 1:43. doi: 10.3389/neuro.01.1.1.003.2007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ribeiro, T. L., Copelli, M., Caixeta, F., Belchior, H., Chialvo, D. R., Nicolelis, M. A., et al. (2010). Spike avalanches exhibit universal dynamics across the sleep-wake cycle. PLoS ONE 5:e14129. doi: 10.1371/journal.pone.0014129

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Schneider, T., and Arnold N. (2001). Algorithm 808: ARfitA Matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models. ACM Trans. Math. Softw. 27, 58–65. doi: 10.1145/382043.382316

CrossRef Full Text

Schnider, T. W., Minto, C. F., Shafer, S. L., Gambus, P. L., Andresen, C., Goodale, D. B., et al. (1999). The influence of age on propofol pharmacodynamics. Anesthesiology 90, 1502–1516.

Pubmed Abstract | Pubmed Full Text

Shew, W. L., and Plenz, D. (2013). The functional benefits of criticality in the cortex. Neuroscientist 19, 88–100. doi: 10.1177/1073858412445487

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shriki, O., Alstott, J., Carver, F., Holroyd, T., Henson, R. N. A., Smith, M. L., et al. (2013). Neuronal avalanches in the resting MEG of the human brain. J. Neurosci. 33, 7079–7090. doi: 10.1523/JNEUROSCI.4286-12.2013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Solovey, G., Miller, K. J., Ojemann, J., Magnasco, M. O., and Cecchi, G. A. (2012). Self-regulated dynamical criticality in human ECoG. Front. Integr. Neurosci. 6:44. doi: 10.3389/fnint.2012.00044

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Strogatz, S. (2006). Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering. Cambridge, MA: Perseus Publishing.

Tagliazucchi, E., von Wegner, F., Morzelewski, A., Brodbeck, V., Jahnke, K., and Laufs, H. (2013). Breakdown of long-range temporal dependence in default mode and attention networks during deep sleep. Proc. Natl. Acad. Sci. U.S.A. 110, 1541915424. doi: 10.1073/pnas.1312848110

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tononi, G. (2008). Consciousness as integrated information: a provisional manifesto. Biol. Bull. 215, 216–242. doi: 10.2307/25470707

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Appendix

FIGURE A1

Figure A1. Model robustness. The figure shows the same analysis presented in Figure 1A (2), but for a VAR-3 model. The results are almost identical to those for VAR-1.