## Abstract

What aspects of neuronal activity distinguish the conscious from the unconscious brain? This has been a subject of intense interest and debate since the early days of neurophysiology. However, as any practicing anesthesiologist can attest, it is currently not possible to reliably distinguish a conscious state from an unconscious one on the basis of brain activity. Here we approach this problem from the perspective of dynamical systems theory. We argue that the brain, as a dynamical system, is self-regulated at the boundary between stable and unstable regimes, allowing it in particular to maintain high susceptibility to stimuli. To test this hypothesis, we performed stability analysis of high-density electrocorticography recordings covering an entire cerebral hemisphere in monkeys during reversible loss of consciousness. We show that, during loss of consciousness, the number of eigenmodes at the edge of instability decreases smoothly, independently of the type of anesthetic and specific features of brain activity. The eigenmodes drift back toward the unstable line during recovery of consciousness. Furthermore, we show that stability is an emergent phenomenon dependent on the correlations among activity in different cortical regions rather than signals taken in isolation. These findings support the conclusion that dynamics at the edge of instability are essential for maintaining consciousness and provide a novel and principled measure that distinguishes between the conscious and the unconscious brain.

**SIGNIFICANCE STATEMENT** What distinguishes brain activity during consciousness from that observed during unconsciousness? Answering this question has proven difficult because neither consciousness nor lack thereof have universal signatures in terms of most specific features of brain activity. For instance, different anesthetics induce different patterns of brain activity. We demonstrate that loss of consciousness is universally and reliably associated with stabilization of cortical dynamics regardless of the specific activity characteristics. To give an analogy, our analysis suggests that loss of consciousness is akin to depressing the damper pedal on the piano, which makes the sounds dissipate quicker regardless of the specific melody being played. This approach may prove useful in detecting consciousness on the basis of brain activity under anesthesia and other settings.

## Introduction

Consciousness is the most prominent feature of the mind, yet its neurophysiological underpinnings remain clouded by a web of disconnected, circumstantial observations. One source of conceptual obfuscation is that, according to a growing consensus (Thompson and Varela, 2001), consciousness is an emergent phenomenon that arises out of the interactions between components of the nervous system but is not reducible to them (O'Connor and Wong, 2012). What, then, are the emergent features of brain activity that are associated with consciousness?

We begin with a simple intuition—the sine qua non of consciousness is responsiveness. Responsiveness to sensory stimuli is the cornerstone of assessment of consciousness in the settings of brain injury, neurologic disorders, and anesthesia. However, responsiveness need not relate specifically to sensory perturbations. It seems likely that, to exhibit consciousness, activity arising in one brain network must be able to perturb activity in a different network. This property of responsiveness is related closely to the concept of stability in dynamical systems theory. Stability governs whether the trajectory of the system in phase space is going to change drastically after a small perturbation or remain essentially unchanged. If the brain dynamics were too stable, then any perturbation would quickly dampen and brain activity would remain essentially unchanged. Conversely, if the dynamics were unstable, then even infinitesimal amounts of noise would grossly disrupt ongoing activity. Therefore, we propose that “dynamical criticality”—the behavior exhibited by systems in the vicinity of a bifurcation between dynamically stable and unstable behaviors (Magnasco et al., 2009)—is essential for consciousness.

General anesthetics present an ideal set of tools to test this hypothesis. Anesthetics act on a variety of receptors (Franks, 2008) distributed among many brain regions (Concas et al., 1990; Devor and Zalkind, 2001; Chen et al., 2005; Jia et al., 2005). Although some, such as propofol, induce the canonical large-amplitude low-frequency oscillations in the EEG, resembling slow-wave sleep (Brown et al., 2010), others, such as ketamine, produce “awake-like” EEG characterized by oscillations in the ∼40 Hz γ range (Maksimow et al., 2006). In fact, despite routine use of EEG-based monitors, some patients awaken during surgery, have postoperative recall (Avidan et al., 2011), and develop posttraumatic stress disorder. Thus, it is currently impossible to reliably distinguish conscious from unconscious brain even in the relatively controlled setting of anesthesia.

To determine whether dynamical criticality is associated with consciousness, we performed stability analysis of high-density electrocorticography (ECoG) in nonhuman primates during reversible transitions into unconsciousness induced with different anesthetic regimens. An array of ECoG electrodes each recording the field potentials of ∼10^{5} neurons (Miller et al., 2009; Ritaccio et al., 2010) was constructed to cover an entire cerebral hemisphere. This allowed us to address simultaneously the macroscopic stability of cortical dynamics and the dependence of stability on the interactions among activities in different cortical regions.

Thus, our overarching hypothesis that dynamical criticality is essential for the responsiveness characteristic of consciousness suggests more specifically that loss of consciousness would be accompanied by a loss of criticality, an implication first preliminarily tested by Alonso et al. (2014). Here we show that dynamical criticality is disrupted universally during loss of consciousness regardless of the anesthetic regimen or specific microscopic features of neuronal activity. Furthermore, stabilization is an emergent macroscopic phenomenon dependent on the correlations in activity among different cortical regions rather than on any specific feature of each individual local field potential.

## Materials and Methods

##### Subjects and data acquisition.

Data from four male monkeys were collected at the Laboratory for Adaptive Intelligence, Brain Science Institute, RIKEN. The datasets are shared in the public server Neurotycho (http://neurotycho.org/). The ECoG array consists of 128 electrodes covering the visual, auditory, somatosensory, and motor areas and the parietal and frontal association cortices (Fig. 1). The ECoG potentials were digitized at 1 kHz. Additional details of the procedure have been described previously (Nagasaka et al., 2011; Yanagawa et al., 2013). ECoG recordings were obtained during the induction of anesthesia starting from the awake state. In this study, we analyze a total of 16 experiments each consisting of reversible induction of anesthesia starting from the awake state. The drug doses used to induce anesthesia are shown in Table 1.

A total of 12 anesthetic inductions were performed using one of the ketamine–medetomidine doses. Four anesthetic inductions were performed with propofol. Ketamine–medetomidine inductions were performed by injecting the drugs intramuscularly, whereas propofol was administered intravenously. Each monkey received more than one anesthetic induction that was separated by at least 1 d.

##### Data processing.

Before the stability analysis (see below), the power spectra of the ECoG channels were examined to identify the potential for contamination with electrical line noise. One channel in monkey M2 was found to have significant noise contamination and was removed from subsequent analysis. All other channels for all other monkeys were included in the autoregressive models (see below). All channels were filtered using a noncausal filter to remove harmonics of the 50 Hz electrical line noise and bandpass filtered between 5 and 500 Hz. Both notch and bandpass filters were implemented using the idealfilter function in MATLAB (MathWorks) to avoid phase shifts.

##### Stability analysis.

ECoG is a multivariate time series whose dynamical properties can be inferred using autoregressive models fitted independently to short time segments as described previously (Solovey et al., 2012). Note that, for our purposes here, we are not interested in reconstructing the dynamics per se. Our focus here is linear stability analysis. This is accomplished by performing a locally linear approximation of the system dynamics in short temporal windows using autoregressive modeling. For every non-overlapping 500 ms window centered at time *T*, a first-order autoregressive model,
was fit independently to the ECoG potentials using a least-squares approach after casting the autoregressive model in the form of an ordinary regression (Neumaier and Schneider, 2001). In Equation 1, *x⃗*(*t*) is a 128-dimensional (number of electrodes in the ECoG array) vector of random variables representing ECoG potentials at time *t*, *A*(*T*) is the 128 × 128 evolution matrix fit to the 500 ms temporal window centered at time *T*, τ = 1 ms (the sampling period), and ε⃗(*t*)is the 128-dimensional error vector at time *t* (*T* − 250 ms ≤ *t* < *T* + 250 ms).

An assumption implicit in this analysis is that the data are both linear and stationary over the time window for which the model is estimated. Autoregressive models have been used previously to fit different neuronal signals including ECoG (Leuthardt et al., 2004). The usual approach is to fit a single autoregressive model to the entire time series (usually several seconds or minutes long). In this case, it is common to use an order *p* autoregressive model *x⇑* = ∑_{i=1}^{p}*A _{i}* ·

*x⇑*(

*t*−

*i*· τ) + ε⇑, in which the

*p*is left as a free parameter. An alternative approach is to fit a complex, biologically inspired nonlinear model to the entire dataset (Boly et al., 2012). However, this approach requires many free parameters and models of the mechanisms responsible for the generation of the signal.

Our approach is purely phenomenological in that it does not make any assumptions about the underlying neural mechanisms and does not require that the signal remain globally stationary. Rather than fitting a single linear model to the entire dataset, we divide the dataset into short non-overlapping time segments and estimate the evolution matrix independently for each time segment. This relies on a much weaker assumption that, although the signal is globally nonstationary, it may still, to a good approximation, be locally stationary over a short time segment.

One significant advantage of autoregressive models is that the stability of the system can be understood in terms of the eigenvalues of the evolution matrix *A* (Neumaier and Schneider, 2001). For a 128-electrode array, using an order 1 autoregressive matrix (128 × 128), we obtain 128 eigenvalues of the system in a short temporal window. Complex eigenvalues define a frequency of oscillation along the corresponding eigenvector. More importantly, for our analysis, the stability of the dynamics (i.e., exponential growth or decay along a corresponding eigenvector) is given by the absolute value of the eigenvalue λ. This can be understood if we consider the formal solution to Equation 1, *x⃗*(*n* × τ) = *A ^{n}x⃗*(0), where τ is the time step (1 ms in our case). Writing the eigenvalue decomposition A =

*U*Λ

*U**, where

*U*is orthonormal, Λ is the diagonal eigenvalue matrix, and * represents the transpose conjugate, we obtain

*x⃗*(

*n*× τ) =

*U*Λ

***

^{n}U*x⃗*(0). The norm of the solution can be expressed as ‖

*x⃗*(

*n*· τ) ‖ = Π

_{k=1}

^{D}

*c*‖

_{k}^{n}*x⃗*(0) ‖, where

*D*= 128 and which we will call the criticality index. If

*c*< 1, the mode is stable: a small perturbation along the eigenvector will decay and return the system to its original trajectory. Conversely, if

_{k}*c*> 1, the mode is unstable: any perturbation along the eigenvector will grow exponentially.

_{k}*c*= 1 corresponds to a critical value at the transition between the stable and unstable modes. In other words, rather than reconstructing system dynamics, in this work, autoregressive models are used to perform linear stability analysis of the ECoG potentials.

_{k}##### Methodology controls: order of the autoregressive models.

Assuming that Equation 1 is a good local approximation, the evolution matrix *A*(*T*) is representative of ECoG dynamics in the corresponding short window around *T*. Thus, complexity emerges from the time evolution of the model parameters, which allows us to follow abrupt changes in neuronal dynamics that accompany induction of anesthesia.

Therefore, it is critical that the model actually captures the dynamics of the ECoG signals within the short windows. Because autoregressive models are stochastic, the goodness of fit can be estimated by the fraction of the total covariance of the dataset that is captured by the model. For 128 channels sampled at 1 kHz, the highest-order model that can be fit to a window of 500 ms is *p* = 3. We evaluated the goodness of fit for orders 1, 2, and 3 by computing the ratio ||*C*_{ε}||/||*C _{x}*||, where C

_{ε and}

*C*are the covariance matrices for the error (Eq. 1) and the ECoG time series, where ||… || stands for the matrix norm defined as the largest singular value of the matrix. Using data from monkey M1, we obtained mean ratio values of 0.011, 0.0088, and 0.0011 for orders 1, 2, and 3, respectively. Using data from monkey M2, the mean over the rest and anesthesia conditions were 0.0082, 0.0044, and 0.00066; however, during the recovery condition, the model performed considerably worse, with means of 0.069, 0.063, and 0.004. As expected, although order 3 yields a better fit (at the expense of three times the number of regressors), order 1 models capture ∼99% of the covariance of the ECoG time series and are therefore a good approximation of the ECoG dynamics. Thus, we use order 1 models for the rest of the analysis of both real and surrogate datasets (see below). Alonso et al. (2014) has shown that the choice of the order of the model does not strongly affect the results. This was confirmed with this dataset as well (data not shown).

_{x}##### Methodology controls: robustness of the conclusions with respect to filtering and window size used to estimate the autoregressive models.

Although the choice of the window for model fitting (500 ms in this case) is somewhat arbitrary, there are some constraints on this choice. If the time segment is too long, the system is likely to deviate significantly from the assumed stationary and linear regimes. Conversely, if the window is made too short, then less data are available to fit the model with confidence. This latter constraint depends on the sampling rate of the data (1000 Hz in this case). Given a 128-electrode grid, the shortest time segment that can be used to fit a first-order 128 × 128 autoregressive matrix is 128 ms. This corresponds to having a single data point for each coefficient. Thus, 500 ms is a reasonable compromise between the two constraints. Note that, in our previous work with human ECoG, the data were sampled at 10 kHz using a 64-channel ECoG array. This allowed us to use 200 ms windows to estimate the models (Alonso et al., 2014).

The choice of the window size limits the bandwidth of frequencies that can be fit by the model. Furthermore, filtering is a nontrivial operation that can, in principle, affect the results. Thus, we varied the size of the window (350–750 ms) and studied the effect of the high-pass filter on the stability analysis using the following procedure.

For each choice of window size and filter, we estimated autoregressive models and computed the distribution of the criticality indices. To study how this distribution evolves in time as the monkey is induced into unconsciousness and subsequently recovers, we computed the *z* statistic using the nonparametric Mann–Whitney *U* test between the distribution in the initial window and each subsequent window. Thus, every experiment is represented by a sequence of *z* statistics. To quantify the degree of consistency between different windows and filter conditions, we computed the correlation between the sequences of *z* statistic. Table 2 shows the correlation coefficients averaged across all 16 experiments. Thus, although we used high-pass-filtered data and window size of 500 ms in this work, within limits, the choice of the window length or high-pass filtering has only minimal effect on the stability analysis.

##### Spectral analysis.

Power spectra were estimated using Thomson's multitaper method (Thomson, 1982) implemented in Chronux (http://chronux.org/; Mitra and Bokil, 2008) with time bandwidth product of 5 (nine tapers). Spectra for each channel were estimated in sliding windows (window duration, 3 s; window step, 1 s). The spectrum in each window was normalized by total power such that all power estimates are expressed as fraction of total power in the corresponding window. Frequency bands were defined as follows: δ, 0–4 Hz; θ, 4–7 Hz; α, 7–14 Hz; β, 15- 30 Hz; γ, 30–100 Hz; high γ, 100–250 Hz. Fraction of power contained in each frequency band was computed by summing the power across all frequency elements in each band.

For the purposes of the “mean spectrum” (see Fig. 6*A*,*B*), the spectrum for each channel was estimated using 10 s windows using a window step of 1 s. Then spectra from different channels within each experiment were averaged. Spectra were then expressed as deviation from the mean awake spectrum computed by averaging spectral windows from 4 min culminating in drug injection.

##### Connectivity.

We did not require the autoregressive matrices to be sparse, and, thus, each electrode is connected to all others. For the purpose of simplifying the content of the matrices, we binarized each one according to the following procedure: we estimated the mean and SD of the distribution of values, *Ā* = 〈*A _{ij}*〉, σ

*= 〈*

_{A}*A*

_{ij}

^{2}〉 −

*Ā*

^{2}(where

*A*is the autoregressive matrix in Eq. 1), respectively, so that only elements 2 SDs above the mean were considered. Formally, the summary matrix

*B*is defined as: For each experiment, we took 200 s of data and binarized each autoregressive matrix using the above procedure.

Then we constructed a summary matrix *S* out of these binarized matrices as follows: *S _{ij}*

_{(i ≠ j)}= 1 if

*B*= 1 for any of the matrices within the 200 s time window. This was performed for control conditions before induction of anesthesia over a 200 s window. Because the kinetics of ketamine–medetomidine and propofol effects were different (see Figs. 2⇓–4), for the ketamine–medetomidine experiments, we started the 200 s “anesthesia” period at 500 s after drug injection, and for the propofol experiment, the “anesthesia” period started at 150 s.

_{ij}The diagonal and off-diagonal elements of the matrix capture two different aspects of neuronal dynamics. Diagonal elements reflect the dependence of activity in each electrode on its own past history, the off-diagonal elements represent the effect of signal at one electrode on future signal in a different electrode. Off-diagonal elements *i*,*j* are shown as arrows originating at *j* and pointing toward *i*. The arrow color changes from blue to red (origin → target).

Note that *S* reflects only the off-diagonal elements. For diagonal elements, we created a summary vector *D*_{i=}∑_{t=1}^{T}*B _{ii}*. This was estimated for the window before,

*D*

_{i}

^{b}, and after,

*D*

_{i}

^{a}, anesthesia. We show

*Ď*=

_{i}*D*

_{i}

^{b}−

*D*

_{i}

^{a}normalized for the purposes of visualization to lie between 0 (red) and 1 (blue). That is, if the diagonal element was higher while awake, it appears as blue; otherwise, it appears as red.

## Results

### Induction of anesthesia results in stabilization of cortical dynamics

Our main finding is that cortical dynamics become stabilized during induction of anesthesia, with both ketamine–medetomidine and propofol (Figs. 2*A*,*D* for representative examples, *B*,*C*,*E*,*F* for group data; 3, 4). In the awake state before drug injection (*t* < 0; Figs. 2 and timing of drug administration is indicated by vertical lines in 3, 4), the criticality indices given by Equation 2 (see Materials and Methods) crowd in the immediate vicinity of the critical value (criticality index = 1). During the induction of the anesthetic state, the number of critical modes decreases: the number of modes with an associated criticality index larger than 0.98 decreased from ∼17.5 to ∼ 13.4% and from ∼15 to ∼4.5% within 4 min after the injection of ketamine–medetomidine and propofol, respectively (Fig. 2*C*,*F*). Although 0.98 is a somewhat arbitrary threshold, as can be seen in Figure 2, *A* and *B* (as well as in Figs. 3, 4) induction of the anesthetized state is accompanied by the decrease in the density of the criticality indices near 1, and, thus, the observed stabilization does not depend strongly on the threshold. Thus, consistent with our hypothesis, dynamical criticality normally present in the awake brain (Solovey et al., 2012) is abolished after the administration of anesthetics.

Although both anesthetic regimens result in statistically significant stabilization of cortical activity, the onset of effect of propofol is faster than that for the ketamine–medetomidine mixture and is less persistent. Note that ketamine and medetomidine were administered intramuscularly, while propofol was given intravenously (Yanagawa et al., 2013). Pharmacokinetic models (Schnider et al., 1999) suggest that the peak in the “effect site” concentration occurs at 1.7 min after intravenous bolus of propofol. Consistent with this, we observed robust stabilization of cortical activity at ∼2 min. The pharmacokinetics of ketamine and medetomidine are more complex, but the peak in the plasma ketamine concentration after intramuscular injection in children occurs at ∼15 min (Grant et al., 1983). Consistent with this, we observed sustained decrease in the number of critical modes 15 min after administration of ketamine and medetomidine.

### Stability is an emergent phenomenon dependent on correlations in cortical activity

Induction of anesthesia is associated with complex changes in the power spectrum (John et al., 2001), phase relationships among different oscillations (Mukamel et al., 2011), and changes in coherence (Cimenser et al., 2011) and correlations (Ku et al., 2011). To address specifically the role of correlations among cortical sites to the global stability, we constructed time-shuffled surrogates produced by independent random constant shift of the time stamps for each electrode, with values drawn from a Gaussian distribution with μ = 0 and σ = 50 ms (Fig. 5*A*). Note that this is a very subtle perturbation to the signal (Fig. 5*B*); it selectively disrupts higher-order correlation but preserves all features of each individual channel taken in isolation, and all of the pairwise correlations between signals are maintained but arbitrarily time shifted by a small amount. However, in every experiment on every monkey, we observed statistically significant differences in the distributions of the eigenvalues obtained in corresponding time windows for the real and surrogate datasets (*p* < 0.05, Mann–Whitney *U* test).

Dynamical criticality is a highly unlikely state of the system and is not generically expected. Thus, the presence of dynamical criticality implies some self-tuning. It has been proposed (for review, see Chialvo, 2010) that many aspects of neuronal activity are found at a critical point of a second-order phase transition. These critical points are characterized by power law distributions of events such as sizes of avalanches in models of sand piles (Bak et al., 1987) and activity avalanches in cultured neurons (Beggs and Plenz, 2003), as well as other biological systems (for review, see Mora and Bialek, 2011). One appealing aspect of this theory is that it has been shown that some systems can spontaneously exhibit such critical behavior without the need for fine tuning of parameters (Bak et al., 1987).

However, the analysis of statistical properties of neuronal activity, such as avalanche size distribution, glosses over the dynamical properties, such as stability. It is possible that the same kind of tuning process that produces power law distributions of states also yields marginally stable dynamics (Magnasco et al., 2009), but currently the relationship between dynamical and statistical aspects of criticality is not well understood (Mora and Bialek, 2011).

Although statistically significant differences between the stability of real and time-shuffled surrogates essentially rule out the possibility that stability depends on any set of features of each signal in isolation, this analysis does not directly address self-tuning, which we address by comparing the internal consistency of the real and time-shuffled surrogates.

It is not easily knowable when, exactly, the monkey loses consciousness after the administration of the anesthetic agent, but the monkey is known to be awake before drug administration. If the awake brain was self-tuned to exhibit dynamical criticality, the distribution of the criticality indices obtained for different time windows during the awake state ought to be stationary, and hence the distribution of eigenvalues for any two distinct windows should remain consistent. In other words, the eigenvalues obtained in any given epoch while the monkey is awake should be drawn from a single probability distribution that remains stable in time, and therefore we should not be able to reject the null hypothesis “the monkey is awake” based on the distribution of the criticality indices. To test this prediction, we computed a reference distribution of the criticality indices over the first 5 min for real and time-shifted surrogate datasets (using 500 ms windows to estimate each autoregressive matrix independently) and compared this distribution with those obtained subsequently during wakefulness using Mann–Whitney *U* test. This procedure is illustrated in Figure 5*C*. Figure 5*D* shows cumulative distributions of *p* values for real (blue) and time-shifted (red) datasets. The abscissa shows the *p* value for rejecting the null hypothesis that the monkey is awake, and the ordinate shows the probability of rejecting the null hypothesis at a given level of statistical significance. The results indicate that the eigenvalues for the real data appear to be drawn from one single (consistent) distribution (probability of falsely rejecting the null hypothesis is low). In contrast, the eigenvalues for the surrogate data fluctuate in distribution; in other words, the surrogate data appear to have been drawn from a distribution that fluctuates haphazardly in time. Thus, given the surrogate dataset, the null hypothesis “monkey is awake” can be readily falsely rejected with high statistical confidence. This result is consistent with the notion that the brain is self-tuned to exhibit dynamical criticality and that this self-tuning involves global correlations. That is, dynamical criticality is an emergent phenomenon.

Note that we chose 50 ms as the width of the Gaussian distribution from which the time shifts were drawn. This preferentially disrupts eigenmodes with higher frequencies. Our objective here was to produce the most subtle perturbation to the data that can nonetheless disrupt the results of the stability analysis. Broadening the distribution of time shifts is expected to further alter the data and exacerbate the differences between real and time-shifted datasets.

### Power spectrum is uninformative of the action of anesthesia

It is well known that induction of anesthesia is associated with changes in the spectral signatures of cortical signals (Brown et al., 2011). Note, however, that in contrast to the uniform increase in stability observed with both ketamine–medetomidine and with propofol, changes in spectral content induced by these drugs are distinct (Fig. 6*A*). Changes observed with propofol are, on average, much larger than with ketamine–medetomidine. In addition to the increase in the slow oscillations <5 Hz, complex changes in other frequency ranges variably accompany induction of anesthesia.

Although it is true that, on average, spectral signatures of cortical signals change with induction of anesthesia, these average quantities conceal both the temporal and spatial heterogeneity of changes in the power spectrum (Fig. 6*C–F*). For instance, although, on average, we observe the expected increase in the low frequency (δ) power with both ketamine–medetomidine and propofol (Fig. 6*C*,*D*, thick lines), each individual channel can deviate quite significantly from the mean behavior (thin lines). Furthermore, the spatial distribution of changes in power is different between the two anesthetic regimens (Fig. 6*D*,*E*). In many frequency bands, profound fluctuations in the spectral content are observed such that the changes in the power spectrum depend strongly on time, the location of the electrode, and the anesthetic regimen. The complexity of changes in the spectral content observed here is consistent with that observed by Breshears et al. (2010) in human ECoG recordings obtained during induction and recovery from propofol anesthesia.

It is not entirely surprising that the distribution of the stability parameter behaves differently from the power spectrum. In a steady-state autoregressive process, there is a complex relationship between the power spectrum of the process (the absolute value squared of its Fourier transform) and the distribution of the eigenvalues in the complex plane. Summarily, each eigenvalue contributes to the power spectrum a Lorentzian distribution, whose center is the frequency given by the imaginary part of the eigenvalue and whose width is given by the (negative) real part of the eigenvalue. The height is given by the projection of the noise term onto the corresponding eigenvector, meaning that each channel can have in principle a sharply different spectrum. Therefore, very similar power spectra can be achieved by very different eigenvalue distributions, and, conversely, the distribution of eigenvalue frequencies by itself does not determine the spectrum. The objects in which we are most interested in this study, the criticality indices, do not determine the spectrum at all. For instance, it is easy to conceive of a system in which a high-frequency mode decays either faster or slower than a low-frequency mode and thus can be either more or less stable.

Non-steady-state autoregressive processes, such as in this work, add additional complexity to the relationship between stability and spectral characteristics because the dependence on the initial conditions cannot be brushed aside, as the system never decays into its asymptotic state before the matrix is changed. The potential presence of unstable modes generate growing oscillations that can only truly be taken into account by remembering that the window is finite and that, at later times, the process will move that eigenvalue onto the stable regime.

### Global reorganization of “functional connectivity”

Although the eigenvalues of the autoregressive matrix (see Materials and Methods) define the stability of the neuronal dynamics, the autoregressive coefficients forming these matrices define a network of interactions between individual cortical areas. Thus, characterization of the structure of the autoregressive matrices could be used to understand changes in directed “functional connectivity” during transition from consciousness into the anesthetized state (Fig. 7).

Although the fitting procedure was blind to the anatomy, note the consistency in the overall flow of connections in the awake state from the primary sensory (occipital) cortices toward the frontal association and inferotemporal cortical areas. Thus, at least to some degree, the autoregressive models fitted to the ECoG capture the patterns of flow of information in the primate cortex.

Loss of consciousness induced with both propofol and the ketamine–medetomidine mixture results in the disruption of the global connectivity pattern. Propofol anesthesia results primarily in the disappearance of connections. In contrast, anesthesia induced with ketamine–medetomidine results in the reversal of direction such that the inferotemporal cortex became the origin of the bulk of the connections. This is consistent with the decrease in the strength of self-connections observed in the occipital cortex (Fig. 7*B*). These results were consistent across experiments and monkeys. Thus, although loss of consciousness results in the disruption of functional connectivity patterns, the specific nature of the disruption depends quite strongly on the identity of the anesthetic agent.

### Frequency dependence of stabilization

The focus of our analysis thus far has been the criticality index given by the |λ|, but generally, the eigenvalues (λ values) of the autoregressive matrices are complex: λ* _{j}* = ρ

_{j}e^{iφj}. Thus, the frequency of the

*j*th eigenmode is given by

*f*=

_{j}*=*

_{j}*t*is the sampling rate (1 ms). Note that ω ≈ 0 corresponds to critical dynamics, ω < 0 corresponds to damped modes, and ω > 0 corresponds to unstable modes. Thus, a more complete analysis of changes in stability examines the changes in the distribution of the eigenvalues induced by anesthetics in the plane spanned by the frequency and damping rate (Fig. 8). Although propofol anesthesia resulted in the increased damping of eigenmodes in the broad range of frequencies from 5 to ∼100 Hz, stabilization induced by ketamine–medetomidine was limited to the high-frequency range greater than ∼50 Hz.

## Discussion

Our main finding is that reversible loss of consciousness is accompanied by reversible stabilization of brain dynamics such that critical oscillations in the awake brain become damped during unconsciousness. The consistency of the stabilization is remarkable given the clear inconsistencies in the spectral characteristics (Fig. 6), directed functional connectivity (Fig. 7), and frequency dependence of damping (Fig. 8) observed with different anesthetics.

In dynamical systems theory, “stability analysis” refers to an evaluation of whether the trajectory of the system would be similar or vastly different after a perturbation. It is performed by evaluating a linear approximation to the full nonlinear dynamics centered at the observed trajectory. This is what our analysis attempts to do: to capture coarse-grained linear components of the dynamics valid for a short time window. As the system moves in its phase space, the linear approximation changes, which is observed as changes to our matrix. An established result of dynamics systems theory is that the dynamics is strongly dominated by the linear components, except when the linear components are critical, in which case the higher-order nonlinear terms decide the ultimate fate of the system. Thus, within our analysis to find that a large number of eigenvalues are critical means that the system is not long-term predictable by this very analysis. What our analysis is capable of demonstrating is that the system hovers over a regime with many critical eigenvalues, the central thesis of this study.

What are the functional consequences of this? Critical dynamics render the system exquisitely sensitive to perturbations (e.g., sensory stimuli), not in an exponentially unstable “butterfly effect” manner but in a more subtle “propagate the perturbation across the system” manner. Consistent with this intuition, our results suggest that the decrease in the responsiveness of an anesthetized subject is heralded by a shift of the dynamics away from the critical regime.

Several seemingly disparate results can be explained by the stabilization hypothesis and integrated into a deeper conceptual framework. Imas et al. (2005) demonstrate that the long-latency component of an evoked potential is suppressed under anesthesia and that the evoked potentials exhibit damping. Perturbation need not be sensory. The effect of transcranial magnetic stimulation (TMS) in the awake subject persists longer and is more complex than in the anesthetized subject (Ferrarelli et al., 2010). This is exactly what one expects when the system moves from a dynamically critical to a stable regime (Yan and Magnasco, 2012). Imas et al. used volatile anesthetics and Ferrarelli et al. used midazolam, neither one of which shares significant similarities in terms of molecular mechanisms or changes in the cortical power spectrum with the ketamine–medetomidine mixture. Thus, loss of dynamical criticality maybe a universal feature of loss of consciousness induced by anesthetics.

Anesthetics may at least in part impinge on the same neuronal mechanisms as sleep (Brown et al., 2010). There is indeed evidence that sleep and anesthesia share essential features in terms of the qualitative dynamics. For example, the effect of TMS is almost identical in the naturally sleeping (Massimini et al., 2005) and anesthetized (Ferrarelli et al., 2010) subjects. Furthermore, the effect of TMS-induced perturbation during disorders of consciousness is similar to sleep and anesthesia (Casali et al., 2013). This leads us to make a bold and tantalizing suggestion that dynamical criticality is an essential requirement for wakefulness in the brain.

### Relationship to functional connectivity

Many methods have been applied previously to show changes in the functional connectivity in anesthesia—correlation between fMRI voxels (Peltier et al., 2005) and EEG channels (Lee et al., 2011), α band coherence (Cimenser et al., 2011), phase lags between EEG recordings (Ku et al., 2011), nonlinear generative models of EEG (Boly et al., 2012), and frequency-specific modes of large-scale communications (Yanagawa et al., 2013) revealed using spectral Granger causality—all show changes with anesthetic-induced unconsciousness.

Here, we use autoregressive matrices to define “directed functional connectivity.” Although our results with propofol are consistent with the notions of “disconnected” brain, under ketamine–medetomidine anesthesia, we observe reversal of the direction of many connections rather than disappearance of connections. One manifestation of this reversal is that, in the brain anesthetized with ketamine–medetomidine, signals from the visual cortex become less “self-driven” and more influenced by the activity of the inferotemporal cortex. Ketamine is known to elicit visual hallucinations and has been used as a pharmacological model of positive symptoms of schizophrenia (Krystal et al., 2003). Visual hallucinations are associated with increased activity in the occipital cortex (Ffytche et al., 1998). Our results suggest that, during visual hallucination, this activity in the occipital cortex is driven by other brain areas rather than sensory experience to generate the illusion or hallucination.

For propofol, stabilization can be related directly to the decrease in functional connectivity. Elimination of the off-diagonal elements in the autoregressive matrices has the effect of decreasing the number of feedback loops that can destabilize the system. The precise nature of the relationship between functional connectivity and stability is less clear for the brain anesthetized with ketamine–medetomidine.

However, regardless of the specifics, our methodology offers significant advantages over other methods. While autoregressive modeling is conceptually closer to the notion of causation proposed by Granger than correlation-based connectivity (Garg et al., 2011), the main advantage in the present context is that it allows for direct interrogation of the global brain dynamics that arise out of this correlation structure through eigenmode decomposition. By fitting autoregressive models to short temporal windows, our method allows for globally nonstationary and nonlinear dynamics and yet takes advantage of the well tractable linear stability analysis.

### Neurophysiological measures of the depth of anesthesia

In clinical and most research settings, depth of anesthesia (DOA) is estimated through changes in the various statistical measures of EEG channels taken in isolation, such as the bispectral index and spectral entropy (Palanca et al., 2009). Despite years of research and use in clinical practice, these DOA monitors fail to prevent intra-operative awareness and postoperative recall (Avidan et al., 2011). One reason for this failure is that, in contrast to the slow progression measured by the criticality index, there is no uniform signature of DOA in the statistical structure of individual recordings. This is consistent with the notion that consciousness is likely an emergent phenomenon that involves correlated activity distributed among multiple brain regions. Here, we suggest that stability is one such emergent property exhibiting robust changes with onset of anesthesia and is independent of individual differences in the spectral properties of individual recordings and of the specific identity of anesthetic agents. Thus, stability analysis may prove useful in characterizing DOA in clinical settings.

Our recordings were confined to a single hemisphere. Although administration of anesthetic to a single hemisphere (Wada test) can trigger confusion, disorientation, and amnesia, it does not render one unconscious (Posner et al., 2007). Furthermore, en bloc resection of an entire cerebral hemisphere does not routinely produce a comatose or a stuporous state (Austin and Grant, 1955).

For systemic drug administration, our stabilization hypothesis predicts that the combined system including both hemispheres should exhibit stabilization during loss of consciousness. Furthermore, we would expect that spatial projections of some eigenmodes would span both hemispheres. These predictions can be tested experimentally in the future. What is less clear is what should happen with direct administration of the anesthetic agents to a single hemisphere. In light of the observations with the Wada test and hemispherectomy, we would hypothesize that modes localized to the ipsilateral hemisphere and those that span both hemispheres may exhibit stabilization, whereas the modes localized only to the contralateral hemisphere may be spared.

### Outlook

One common interpretation of dynamical criticality in the brain is that a system can only achieve optimal information processing in the vicinity of a phase transition (Langton, 1990), yet most physical systems enjoy “physical stability”: they preserve their qualitative behavior under perturbations, i.e., the defining parameters of most systems are found away from a bifurcation. In stark contrast to this, our results and those of others (Solovey et al., 2012; Alonso et al., 2014) suggest that a significant fraction of cortical dynamics is found close to a bifurcation. This conclusion is not unique to cortical dynamics. The most compelling case for dynamical criticality has been made for the auditory periphery in which Hopf bifurcation is thought to give rise to some of the fundamental aspects of hearing, such as compressive gain, frequency tuning, and spontaneous otoacoustic emissions (Hudspeth et al., 2010).

The very fact that a system is found close to a bifurcation implies self-tuning because bifurcations occupy a small region in the parameter space. Although our analysis is phenomenological and thus does not allow for the direct interrogation of underlying mechanisms, in theoretical models (Magnasco et al., 2009), simple and biologically plausible anti-Hebbian synaptic plasticity leads robustly to a dynamically critical regime. It is possible that the slow increase in stability seen with induction of anesthesia is a manifestation of the changes in the synaptic plasticity.

## Notes

Supplemental material for this article is available at http://sur.rockefeller.edu/Plone/lab-members/leandro-alonso/. Animation showing evolution of stability analysis during reversible loss of consciousness. This material has not been peer reviewed.

## Footnotes

This work was supported by National Science Foundation Grant EF-0928723 (M.O.M.) and National Institute of General Medical Sciences Grant 1K08GM106144-01 (A.P.).

- Correspondence should be addressed to Dr. Alex Proekt at his present address: Department of Anesthesiology and Critical Care, Perelman School of Medicine at the University of Pennsylvania, Philadelphia, PA 19104. proekt{at}gmail.com