Introduction

Wakefulness in healthy subjects is linked to motor behavior and consciousness. Under some pathological or physiological conditions, behavioral response might break down, but consciousness doesn’t necessarily vanish. Anesthetic drugs (Fig. 1a) are used to medically induce comatose states of various depths associated with unresponsiveness. Common anesthetic drugs like Xenon or Propofol lead to an unresponsive and unconscious state resembling NREM sleep(1). Others, such as Ketamine, can cause an unresponsive yet active state similar to REM sleep where conscious experience remains possible even without behavioral response(2). These states of disconnected consciousness can arise following brain injury and pose challenges in clinical practice for consciousness assessment and prognosis evaluation(3, 4). In an attempt to solve this problem, recent studies have developed a technique based on transcranial magnetic stimulation (TMS)(5), in which the brain is stimulated by TMS pulses(6) and evoked responses are recorded with simultaneous electroencephalography (EEG). It is shown that the evoked response can qualitatively differ depending on the patient’s state (awake, REM, anesthesia, etc.), and characterized by a simple scalar metric: the perturbational complexity index (PCI). PCI measures the complexity of the distributed spatiotemporal response in the brain, with values close to 1 when the complexity is high, and closer to 0 for simpler responses. Notably, increasing evidence suggests the PCI index as a reliable marker (Fig. 1b) to differentiate between awake individuals and those in NREM sleep or an anesthetic-induced coma(7, 8).

Fig 1:

Dynamical systems theory is becoming increasingly popular in neuroscience(9) to study the evolution of spontaneous brain activity, i.e. resting state. Resting-state activity is highly organized and displays structured patterns of regional activations in large-scale recordings across imaging modalities. Some co-activation patterns are consistent across individuals and are known as functional brain networks(10). Over time, these activation patterns can emerge and reappear, suggesting a rich temporal structure with unique characteristics across tasks and individuals(11). Functional Connectivity (FC) measures capture the topographic organization of these patterns(12) and are static views of the functional coordination between brain regions. Their temporal organization is estimated by dynamic Functional Connectivity (dFC), generally measured by correlating FC patterns at distinct time windows(13, 14). In this context, data complexity relates to the range of activity patterns seen over time(15, 16), defining a functional repertoire (Fig. 1c). A diverse functional repertoire, supporting complex dynamics, is characterized by the exploration of numerous brain states, found in specific dynamical regimes(13). Further evidence suggests that the brain is self-regulating around the transition between ordered and disordered phases(17, 18) where a high fluidity results in a rich temporal structure. In this study, we aim to examine the connection between perturbational complexity and fluid dynamics near this transition

For that purpose we use a large-scale brain network model The Virtual Brain (TVB)(19) to simulate non-invasive brain stimulation and spontaneous activity (Fig 1. d). Various dynamical regimes are explored by varying a global coupling parameter that scales the influence of the structural connections as compared to the local regional dynamics. Earlier studies indicated that a properly scaled global coupling parameter can produce phenomena of co-activations between brain regions, realistic FC patterns, and non-trivial dFC, as seen in fMRI experiments(20, 21). Here, we quantify the response to simulated TMS across dynamical regimes using the simulation PCI (sPCI) and connect it to dFC metrics. Our hypothesis is that simulations can only produce high sPCI and non-trivial dFC when the system is fine-tuned to a particular fluid regime. We then verify our findings by comparing empirical spontaneous and stimulation data of 15 subjects, each acquired during wakefulness and various states of consciousness (Propofol anesthesia (N=5), Xenon anesthesia (N=5), Ketamine anesthesia (N=5)).

Results

Large-scale brain model and dynamical regimes

To explore the relationship between resting state dynamics and brain responsiveness, we used a large-scale brain network model(22). This model consists of 84 cortical and subcortical regions, each of them being conceptualized by a neural mass model (NMM). The NMM is composed of a set of stochastic nonlinear differential equations, which describes the evolution of both the average firing rate r(t) and the average membrane potential V(t) of an ensemble of quadratic integrate and fire (QIF) neurons(23). The equations ascribe the dynamics of the node to two fixed points. One is an oscillatory, high-firing state, the upstate (inward spiral), and the other a low-firing rate downstate (stable fixed point). Every node is influenced by a stochastic noise component, allowing it to transition between these two states. When embedded into a brain network, the state of each node is also influenced by the afferent input from its neighbors through the connectome(24). The intensity of input any region receives is weighted by the size of the bundle of white fiber tracts, as estimated by diffusion weighted imaging (DWI). For our study, due to lack of DWI data, we employed a generic connectome from a healthy subject obtained from the Human Connectome Project (HCP)(25). Examining the system at the network level, spontaneous activity can present a variety of dynamics depending on two global parameters: the global coupling G (which scales the connectome) and noise intensity σ (white noise with zero mean and standard deviation σ) (Fig 2a). In the low coupling regime (Fig 2a, left side), regions of the network are weakly connected and mostly behave independently, mostly driven by noise. Hence, when the noise is low (Fig 2a, bottom left corner), the activity is generally low and sparse, as most of the nodes are in the down state, with some rare and transient activations. As noise increases (Fig 2a, top left corner), activations are more frequent, yet stochastic and non-coherent. In contrast, when global coupling is set to a high value (Fig 2a, right side), brain regions are tightly connected, and the activity of each node is heavily entrained by the rest of the network. In this regime, most nodes tend to stay in the upstate due to high average input (proportional to G, since I(t) = f(G, W)). Only nodes that are weakly connected to the rest of the network will stay in the downstate despite high coupling. This behavior gives rise to a stable pattern of regional coactivations that involve characteristic subnetwork mirroring the structural connectivity(20, 26) and resulting in stereotyped activity on a slow timescale. However, in-between these two regimes, lies an intermediate regime, around a working point, (Fig 2a, middle) where cascades of coactivations of different sizes and durations occur in a complex and unpredictable fashion. Yet these coactivations patterns are not purely stochastic(27). They are often repetitive, revisiting the same neighborhood in state space(28).

Fig 2:

Modeling stimulation and spontaneous activity

Brain stimulation was modeled by adding a transient external current to the membrane potential, mimicking a single pulse stimulus to a cortical brain region. When perturbed, the whole system shifts from its current state as the stimulus propagates through the network until the system settles to a new stable state. In line with the PCI experimental protocol, we sampled from multiple initial conditions and regions, presenting the maximum sPCI for each regime (i.e., each {G, sigma}). For each simulation, we measured the complexity of activity over 10 seconds post-stimulus. We then compared this to the complexity of the baseline (i.e. the same simulation without a stimulus) using a ratio (see Methods). This helped us determine the relationship between the current dynamical regime and the system’s ability to integrate a perturbation flexibly. We found that the maximum increase in complexity is reached when global coupling is set around G = 0.55, and complexity itself can be increased five-fold with respect to baseline (Fig. 2b). If the system operates far from this working point, the stimulus has minimal impact on global activity, meaning no notable change in complexity. The visual inspection of stimulation time series (Fig. S1) across different regimes indicates that stimulation can trigger (or disrupt) coactivation cascades when the system is near maximum fluidity (G=0.55). This behavior is behind a symmetry between minimum and maximum sPCI as one can see that they share the same distribution in the parameter space. If sPCI is greater than one (lower than one) the trajectory of the system moves towards (away from) an attractor thus triggering (destroying) a coactivation cascade. When sPCI is close to one, either stable attractors are absent, or it’s challenging to deviate from the attractor, but the resulting response won’t exceed the baseline activity’s complexity.

To compare responsiveness to resting state activity we examined the fluidity of the system, defined by the variance of the upper triangular part of the dFC. Fluidity captures the recurring patterns of co-activation across time and, more specifically, the switches between highly correlated and non-correlated patterns. Even though direct evidence is sparse due to challenges in probing human brain activity’s manifold, fluidity is believed to reflect the complexity of the energy landscape and the number of attractors(13, 28). The distribution of fluidity in the model’s parameter space (Fig 2c top-left) shows a peak near the working point value, reflecting recurring cascades during simulation, and decreases monotonously as the regime is drawn away from it.

To highlight that fluidity is a property of the dynamics and is not just tied to a single metric, we investigate and present three additional metrics characterizing spontaneous activity. We first calculate the raw Lempel-Ziv (LZ) complexity of each regime after binarization of the time series. LZ complexity counts the number of unique patterns (sequences of bits) found in a binary string, allowing for patterns of different sizes. A simple fixed threshold was used to binarize each time series into a sequence of upstates and downstates. LZ complexity (Fig 2c bottom-right) is found to increase abruptly as global coupling crosses the working point separating a low complexity-low coupling regime from a high complexity-high coupling regime. Second, we extracted the size of the functional repertoire of the model, by counting the number of unique configurations that the network takes over time. As for LZ complexity, there is a separation of the space around the working point by a sharp increase of the functional repertoire (Fig 2c bottom-left), although this separation shifts towards lower coupling when noise increases and is lost for regimes with high noise. Lastly, we measured the bursting potential of the system. Bursting potential aims at capturing the fastest change in global activity of the network and it is computed as the maximum of the gradient of the residual sum of squares between sources (see Methods). Again, a similar pattern is seen with an increase around the working point (Fig 2c top-right). These tight relationships between simulation features (on spontaneous and stimulated activity) and the tuning of global coupling suggest that the same fluid dynamical regime underpins both optimal brain responsiveness and optimal spontaneous activity. Under the assumption that the global coupling models a free parameter in the brain which is subject to change upon internal and/or external states, global dynamical changes should be found for different brain states. Based on this reasoning, we sought to retrace the signatures of these same dynamical regimes in real data in spontaneous EEG recordings during wakefulness and anesthesia.

Data analysis

To do so, we used previously published experimental data, provided by(8). 15 healthy subjects underwent a TMS-EEG protocol and a spontaneous EEG during anesthesia (Xenon, Propofol, Ketamine) and wakefulness. These anesthetics were posited as empirical models of loss of consciousness. Reports upon awakening were used as a ground truth for conscious experience and it was found that subjects under Xenon or Propofol had no recollection of conscious experience whereas Ketamine induced anesthesia was always associated with reports of consciousness (i.e. dreaming). During the TMS-EEG protocol multiple sessions of stimulation each consisting of multiple trials per stimulation site were performed. A PCI value (across trials average) was obtained for each session and the maximum PCI across sessions was retained. Only the maximum value of PCI for each subject and condition was available to us. A PCI at 0.31 was benchmarked as a threshold for consciousness (7). As published previously, PCI successfully distinguishes conscious from unconscious states. Complexity measured with PCI is systematically higher during wakefulness than under anesthesia for Xenon and Propofol drugs, but not for Ketamine induced anesthesia which is grouped with wakefulness (PCI > 0.31). Based on modeling results linking perturbational complexity and spontaneous activity metrics to the same underlying dynamical regimes modulated by global coupling, we hypothesized that variations of complexity should also be found at rest and associated with consciousness.

We then studied the four resting state metrics described previously and present their association with consciousness to demonstrate that consciousness is associated with more complex activity at rest. In fact, for all subjects, at least one metric can perfectly distinguish wakefulness from anesthesia. We first extracted fluidity of the recordings on a broad band spectrum (highpass filter at 0,5Hz and lowpass at 60Hz). Dynamic functional connectivity was calculated using circular correlation between electrodes with a sliding window approach (1 second window, 100 ms step size). Then, for each subject and recording we used a bootstrap procedure to sample uniformly with replacement 50 segments of 1 minute. Complexity, the size of the functional repertoire and bursting potential were also computed on each segment. When binarization was applied, we z-scored each electrode independently and thresholded at 3 standard deviations. This yielded a distribution of 50 points for each subject and condition (anesthesia or wakefulness). We present in Fig. 3 the results, subject-wise, for the 4 resting state metrics we studied: fluidity, complexity, functional repertoire and bursting potential between anesthesia and wakefulness conditions. In Fig. 4, we show classification results when grouping participants either on anesthesia (vs wakefulness) or conscious report (vs no report). We cross plotted each metrics against the PCI on the top and show classification accuracies of all metrics after fitting an SVM classifier with a linear kernel. For these last results we computed complexity and the size of the functional repertoire on the full recordings and adjusted by the length of the signal.

Fig 3:

Fig 4:

According to our hypothesis and modeling results, we expected higher values of fluidity during wakefulness than anesthesia (Fig. 3a). The separation between conditions was perfect for Propofol and Xenon anesthesia. Ketamine results were inconsistent, with subject 2 showing a reversed pattern (higher fluidity during anesthesia). The classification by conscious report was perfect (100% accuracy), just like with the PCI (Fig. 4a). Based on the results obtained on perturbational complexity, we expected to find systematically more spontaneous complexity during wakefulness than anesthesia (Fig. 3c). For Propofol anesthesia all subjects had a good separation and for Xenon, 3 subjects were well classified but two had overlapping distributions of complexity. For Ketamine, subject 4 had overlapping distributions when the other 4 showed good separation though the relationship was inverted for subjects 3 and 5 with less complexity found in wakefulness. In line with(20), we hypothesized that bursting potential was associated with the flexibility of the brain to switch between distant functional states and expected this ability to be lost during anesthesia. The results for bursting potential were satisfying for the Xenon and Propofol conditions but poor for all subjects under Ketamine (Fig. 3b). As expected, when the separation occurs between wakefulness and anesthesia, the bursting potential is consistently higher during wakefulness. It is likely hinting at different mechanisms of action of the drugs on neural dynamics, in which bursting potential is not altered by Ketamine. Lastly, as the size of the functional repertoire reflects the extent of explored functional patterns, we were expecting to find it associated with wakefulness, when the brain is more solicited and needs to reconfigure more often. Indeed (Fig. 3d), it is the metric that permits a perfect separation between wakefulness and anesthesia for all subjects and conditions. The number of states explored during the recordings was consistently larger during wakefulness than anesthesia for all drugs and subjects.

Discussion

Understanding the neural substrates of consciousness stands as one of the most significant scientific challenges encompassing both philosophical and theoretical dimensions. This challenge also stems from clinical practice, such as the need to characterize disorders of consciousness and their prognosis(3, 4, 29). The detection of consciousness independently of sensory processing and motor behavior has recently improved with the development of the Perturbational Complexity Index that relies on the analysis of TMS-evoked potentials (7). Despite its remarkable performance TMS-EEG can be difficult to implement in clinical practice and it requires complex equipment. In this work, we explored both the dynamics of spontaneous activity and responsiveness in a large-scale brain model, and found strong evidence that common network-level factors govern the two. We demonstrate using a large-scale brain model that a similar dynamical regime of activity can cause perturbational complexity and complex (fluid) spontaneous activity. We then showed using empirical data that spontaneous activity analyzed through appropriate metrics can be effective in assessing consciousness in wakefulness and anesthesia.

A large body of evidence explores a relationship between the complexity of brain signals and conscious experience, as recently reviewed(30). The first theoretical argument started from a phenomenological description of conscious experience that led to theorizing the expected properties of the underlying neural activity(31, 32). In this view, the two building blocks of consciousness are functional integration (each experience is unified as a whole) and functional differentiation (each experience is unique and separated from others). These theoretical findings lead to the development of metrics rooted in information theory to estimate these functional properties on neural data. Of all these metrics, the PCI (relying on Lempel-Ziv complexity and entropy) revealed an interpretable and practical link between the spatio-temporal complexity of the propagation of a perturbation in the brain and the level of consciousness. In parallel, a large literature is dedicated to understanding brain function from the perspective of dynamical systems theory. In this field, the brain is viewed as a complex open system far from equilibrium composed of a large number of coupled components (e.g., neurons at the microscale or neural masses at the mesoscale) whose interactions are responsible for the emergence of macroscopic patterns. In such systems, drastic changes such as phase transitions or bifurcations can occur spontaneously or under the influence of a global parameter, with major effects on the dynamics. Brain simulations like the one presented in this work support the idea that a rich dynamics is possible when the system’s parameters are fine-tuned around specific values(33). In fact, although debated, one prominent theory in neurosciences posits that the brain self-regulates near a critical regime i.e. in the vicinity of a phase transition(17, 18). Simply put, this property allows the brain to be flexible enough to reconfigure and adapt dynamically to a changing environment all while being stable enough to engage in complex sustained activities. Among other things, the critical regime is known to maximize information processing(34) or related to cognitive processes(35), and criticality in the brain is also possibly linked to consciousness(36–38). Dynamical systems theory in neurosciences has also found a strong paradigm with the concept of manifolds. It’s been shown that a high dimensional nonlinear dynamical system can be described in a low dimensional space. This concept was first developed as synergetics(39), but restricted to systems living close to a bifurcation, and later generalized to symmetry breaking systems(40). Recent applications of this concept can be found in the domain of motor control(41) or to sleep/wake studies(42, 43). The link between dynamics and complexity is still under theoretical exploration(16) and, here, we bring an empirical argument in that direction. Indeed, by using a large-scale brain model, we were able to control the dynamical regimes of the network and show that fluidity and responsiveness are maximal within the same parameter range. Fluidity captures the complexity of the manifold of the brain activity (i.e., number of attractors) and responsiveness the complexity upon perturbation. It could be argued that the noise level is not similar in both settings but noise is only present to ensure continuous exploration of the manifold without affecting its topology(15). Complexity had already been shown to vary across conscious states(44) and the 1/f slope of the EEG spectrum was also correlated with PCI on the same data we used(45). Our validation results are in line with previous work on spontaneous activity and consciousness. In addition, we demonstrate that it’s possible to distinguish ketamine induced unresponsiveness from wakefulness at the individual level with a single metric.

Of important consideration, global coupling is a phenomenological parameter without any biophysical meaning. Several studies have focused on biological mechanisms of consciousness, for instance on the role oscillations in consciousness or the E/I balance. Here, we don’t address these questions but rather give a proof of concept that large-scale brain models can help understanding the dynamics related to brain function. We used a model derived from QIF neurons(23) that lacks biological parameters such as ion concentration or synaptic adaptation. Related work has shown that adaptation successfully reproduces dynamical regimes coherent with NREM and wakefulness (Cattani et al 2023) with corresponding realistic PCI values(46). Similarly, using anatomically constrained parameters for brain regions has already been shown to increase the predictive value of brain network models (Wang et. al. Sci Adv 2021, Kong et al. Nat Comm 2022). Nevertheless, we demonstrate that even the symmetry breaking caused by the connectome is sufficient for setting the global working point of the brain, which then links the brain’s capacity for generating complex behavior in the different paradigms, that is, rest and stimulation.

One major limitation of our study lies in the algorithm we used to assess complexity in the model. The empirical PCI calculation is done on the average evoked response following stimulation that only lasts a few hundred milliseconds before returning to baseline activity. In the model we spanned this calculation over ten seconds thus capturing a slower dynamic than in real data. Nevertheless, we believe it doesn’t affect our statement and in theory the neural mass model we used doesn’t have a specific timescale. In future work, some caveats of our work could also be addressed and remedied. First, the imperfect separation of groups for some of the metrics could be improved by personalized brain modeling. The size of the functional repertoire, for instance, distinguishes wakefulness from anesthesia (including ketamine) but lacks predictive power at the group level. This could be solved using a personalized model including structural information and parameter inference. Second and last, more realistic parameters could be included in the models such as neuromodulatory pathways(47, 48) to improve explanatory power.

Materials and Methods

Connectome based modeling

The connectome used in this work was obtained from the Human Connectome Project (HCP)(25), it was extracted through Diffusion Tensor Imaging (DTI) and belongs to a healthy subject. This connectome captures the tracts weights and lengths between regions of interest (ROI) according to the Desikan Kiliani Atlas(49). In the end 84 cortical and subcortical regions are connected through an undirected graph. Each region of interest (ROI) of the brain was modeled by a neural mass model developed by Montrbrió et al(23). This model is a derivation of the mean firing rate r and membrane potential V of an all-to-all coupled network of quadratic integrate and fire (QIF) neurons. It is exact in the thermodynamic limit (N 0.55 0.55, with N the number of neurons) and is described by a set of two differential equations:

η is a neuronal excitability parameter that follows a Lorentzian distribution for which the mode appears in the mean field equation and Δ is the scale (inter-quartile interval) of the same distribution. Parameter J is the weight of synaptic integration of the QIF neurons. And I(t) is an external current All neural masses were constructed equal with parameter values: , Δ = 1.10 and j = 1.0 and J = 15.0. Regions models are connected through their membrane potential by setting:

Where, W_nm is the weight (white fiber tract width) and τ_nm is the tract length connecting the two regions n and m.

Data

We retrieved data from studies(8, 45), where 15 healthy subjects between 18 and 28 years old (5 males) were randomly assigned to one of three groups (5 subjects per group) to undergo anesthesia induced by either Xenon, Propofol or Ketamine. For each subject, a 60 channel EEG recorded spontaneous activity of the brain followed by a TMS-EEG session for PCI assessment, once before drug administration (wakefulness) and one under anesthesia. All subjects were unresponsive during anesthesia as measured by the Ramsay(50) scale score 6. Participants were asked to report any previous conscious experience upon awakening. Xenon and Propofol anesthesia led to no report of consciousness while Ketamine allowed for retrospective conscious reports (vivid ketamine-dreams). We retrieved the spontaneous EEG recordings for every subject during wakefulness and during anesthesia, as well as the maximum PCI value from the TMS-EEG session and the absence/presence of conscious experience. For more details in the experimental setup and preprocessing steps, see (8).

Complexity

For complexity assessment in spontaneous activity, we calculated the Lempel-Ziv complexity(51) denoted as LZ_L on binarized activity. For synthetic data, the bistable nature of the Montbrió model allowed us to binarize firing rate time series (both spontaneous activity and after stimulation) by simple thresholding. The threshold (r_thresh = 0.7) was determined empirically by visual inspection with the intent to limit false positive and false negative rates of upstate activity. For empirical data, the binarization of spontaneous EEG done after z-scoring each electrode and thresholding at 3 standard deviations.

To capture the complexity of the activity after stimulation we computed an index strongly inspired by the Lempel Ziv version of the perturbational complexity index(5). We define the simulation PCI (sPCI) over a time series of length L as:

where LZ_L is the Lempel-Ziv complexity of a series of length L and 0.55_L is the asymptotic Lempel-Ziv complexity for a random series of the same length L, which is proportional to the source entropy that would have generated the data with the same average activity p. For each {G, σ} couple, we stimulated a given node node_i at a given time (or trial) t to obtain a sPCI_stim value over a window of ten seconds of activity:

To distinguish the deterministic response of the system from the spontaneous activity, each of these pci_stim value was matched with a baseline value:

which corresponds to the complexity associated with the spontaneous activity of the realization of the exact same simulation (same parameters, same noise seed) but without the stimulus. The ratio of those two quantities was considered as the change in complexity due to the stimulus:

A value of sPCI > 1 (respectively sPCI < 1) indicates an increase (respectively a decrease) in complexity in response to the stimulus. 17 nodes were randomly selected (for computational limitations) among the 68 cortical nodes and 140 trials were done for each node near the working point G ∈ [0.52, 0.57] and σ ∈ [0.022, 0.039], and 6 trials elsewhere. We then selected the maximum sPCI across time and nodes:

Empirical maximum PCI for every participant was provided by Sarasso et al(8) and was also based on Lempel-Ziv complexity. Details about the implementation and calculation can be found in the original publication.

Fluidity

The extent of the repertoire of functional patterns explored spontaneously by a network is measured by fluidity defined as the variance of the upper triangle of the dynamic functional connectivity matrix. A sliding window (of length τ_win ) of functional connectivity (FC) is applied:

This yields a stream of FC matrices across the length of the signal s(t). Dynamic functional connectivity is then computed as the correlation between the upper triangular part of those FC matrices to extract recurrences of functional patterns across time:

To finish, fluidity is the variance of the upper triangle of the dFC matrix, after removal of overlapping windows:

In practice, the diagonal offset to compute fluidity is given by the formula:

With n_overlap the number of overlapping time points between two adjacent windows.

In synthetic data functional connectivity was based on the zero-lag Pearson correlation between the firing rate activity of every pair of nodes of the network:

In empirical EEG recordings functional connectivity was assessed by the circular correlation coefficient among all pairs of electrodes (N_con = 1770 for 60 channels) over a sliding window of length τ_win 0.55 We computed the Hilbert transform of each electrode to extract the instantaneous phase φ_i(t) and mean phase within each window on a broad ban spectrum (0.5-60Hz). Then between for each pair of electrodes i, j we computed functional connectivity:

Functional repertoire

To measure the size of the functional repertoire in empirical data we first extracted avalanches on z-scored and binarized EEG recordings (threshold=3s.d). We used a bin size of 1 and verified that branching ratio values were all close to one. The size of the functional repertoire was then defined by the count of unique avalanche configurations during a recording(52). In synthetic data, avalanche extraction was not always possible depending on the dynamical regime. We simplified the procedure by counting the number of unique configurations across all time points after binarization.

Bursting potential

Various modalities have shown that structured brain patterns emerge in short-lived bursts of neuronal activations(53–56). To capture these events in our data, we first z-score the data along time, and then perform the root-sum-square (RSS) across electrodes. The RSS profile has values described by a fat-tail distribution and presents a bursting profile characterized by rare sudden jumps. Similar to the PCI, which measures the potential for the system to express complex responses by selecting the maximum of perturbational complexity across an array of stimulations, here we define the potential of spontaneous bursting for the system by selecting the maximum variation of the RSS, which we dub Bursting Potential (BP):

Classification

To test for classification accuracy of the resting state metrics we used a Support Vector machine algorithm with a linear kernel. Two different classifiers were evaluated. The first, grouping all participants in the wakefulness condition (N=15) against anesthesia condition (Ketamine N=5, Xenon N=5, and Propofol N=5). The second, grouping participants in conscious report (wakefulness N=15, and Ketamine anesthesia N=5) against no report (Xenon N=5, and Propofol N=5 anesthesia). Each classifier was trained using one feature at a time to find the optimal threshold for classification: PCI, fluidity, complexity, size of the functional repertoire, and bursting potential. All features were calculated on the full duration of each recording. Complexity and Functional repertoire features were normalized by the length in minutes.