A fast and general method to empirically estimate the complexity of distributed causal interactions in the brain

A novel metric, called Perturbational Complexity Index (PCI), was recently introduced to assess the capacity of thalamocortical circuits to engage in complex patterns of causal interactions. Once validated and calibrated in a large benchmark, PCI showed high accuracy in detecting consciousness in brain injured patients. In its original formulation, the index was computed by stimulating the brain with transcranial magnetic stimulation (TMS) and using the Lempel-Ziv algorithm to quantify the spatiotemporal complexity of cortical responses, as derived by the inverse solution of high-density EEG (hd-EEG). Despite its accuracy in assessing consciousness, the original formulation of PCI was limited by its reliance on source estimation, which not only renders its calculation timeconsuming, dependent on elaborate experimental setups and offline processing, but also restricts its extension to other types of brain signals beyond TMS/hd-EEG recordings. Here, we address these limitations and introduce PCIST, a fast and simple method for estimating perturbational complexity of any given brain response signal. Our approach combines dimensionality reduction with a novel metric of complexity derived from recurrence quantification analysis, in which state transitions (ST) present in the signal’s principal components are quantified in order to estimate spatiotemporal complexity directly at the sensors level. PCIST was tested on a large dataset of TMS/hd-EEG recordings obtained from 108 healthy subjects and 108 patients with brain injuries, and also applied to sparse intracranial recordings of 9 patients undergoing intra-cerebral single-pulse electrical stimulation. The new method performs with the same accuracy as the original formulation, but the index can be computed in less than a second, requires a simpler setup and can be generalized to sparse intracranial recordings. PCIST advances towards the development of a clinical bedside index of consciousness and can be used as a general tool to study the brain’s complexity across scales, experimental models and species.


Introduction
Measures of brain complexity have recently begun to move from the realm of theoretical neuroscience [1][2][3][4][5][6] into the field of experimental neurophysiology to study differences between global brain states, from wakefulness to sleep and anesthesia [7][8][9][10]. Further, measures of brain complexity have been considered as useful paraclinical indices to assess consciousness at the bedside of brain-injured patients [11][12][13][14][15]. In this spirit, a novel strategy based on quantifying the global effects of direct cortical perturbations was recently introduced [16]. This approach is motivated by the general theoretical principle that a brain's capacity for conscious experience relies on its ability to integrate information [5]. In this view, a key mechanism of consciousness is the ability of different neural elements to engage in complex patterns of causal interactions such that the whole system generates information over and above its parts.
Practically, in order to estimate the amount of causal, irreducible information that a system can generate, a general procedure was implemented based on two steps: (i) locally perturbing the system in a controlled and reproducible way to trigger a cause-effect chain and (ii) quantifying the spatiotemporal complexity of the ensuing deterministic response to estimate information. The original implementation of this perturb-and-measure approach [16] involved (i) stimulating the brain with transcranial magnetic stimulation (TMS) and (ii) computing the algorithmic (Lempel-Ziv) complexity of the resulting patterns of activations at the level of cortical sources derived from the inverse solution of high-density electroencephalographic (hd-EEG) responses; this metric will be henceforth referred to as Lempel-Ziv Perturbational Complexity Index (PCI LZ ).
Albeit macroscopic and coarse, PCI LZ provided maximum (100%) accuracy in detecting consciousness in a large (n=150) benchmark population of subjects who could confirm the presence or absence of conscious experience through immediate or delayed reports. PCI LZ was lower in all unresponsive subjects who did not report any conscious experience upon awakening from non-rapid eye movement (NREM) sleep or midazolam, xenon, and propofol anesthesia, and was invariably higher in conditions in which consciousness was present, including awake controls, conscious brain-damaged patients and subjects who were disconnected and unresponsive during dreaming and ketamine anesthesia but retrospectively reported having had vivid conscious experiences upon awakening [17,18]. Once calibrated on the gold-standard of subjective reports, PCI LZ measurements performed at the bedside of non-communicating subjects with brain injuries offered high sensitivity (94%) in detecting minimally conscious patients and allowed identifying a significant percentage (about 20%) of vegetative state/unresponsive wakefulness syndrome (UWS) cases with high brain complexity, who had a higher chance of eventually recovering consciousness [17].
While PCI LZ performs with unprecedented accuracy, it also has practical drawbacks and limitations. First, PCI LZ can only be computed on spatiotemporal matrices of cortical activations that are obtained after an intensive processing of TMS/hd-EEG data, including forward modeling [19], source estimation [20] and permutation-based statistics at the single-trial level. All these steps imply a complicated and lengthy off-line analysis pipeline that hinders the dissemination of the method and its application as a routine clinical bedside tool. Clearly, the possibility of estimating perturbational complexity directly at the level of EEG sensors, may have critical advantages: not only it would render the analysis process faster (ideally, on-line), easier to standardize and immune to the technical caveats of source modeling, but it would also allow the use of simplified and cheaper set-ups (i.e. not requiring hd-EEG and subject-specific MRI scans).
A second important drawback of PCI LZ is its limited application to signals other than TMS/hd-EEG evoked potentials. Intracranial stimulations/recordings in humans [21,22] and in animal models [23][24][25][26] as well as intra and extracellular responses recorded from cortical slices [27,28] offer an unprecedented range of opportunities to elucidate the relationships between neuronal dynamics, network complexity and consciousness [29]. However, because PCI LZ relies on EEG source estimation, its extension to other types of recordings, such as sparse matrices of intra-cerebral stereo EEG recordings and in vivo/in vitro local field potentials, is not straightforward [27].
Here we propose a novel measure of perturbational complexity that bears conceptual similarities with PCI LZ but is much faster to compute and in principle generalizable to any type of brain signal evoked by perturbations or event related. Conceptually, we started from the notion that the binary sequences of activation and deactivations which are compressed by PCI LZ can be considered as sequences of transitions between different states: a "response state" and a "non-response" or "baseline state". Thus, one should expect to find high values of perturbational complexity in systems that react to the initial perturbation by exhibiting multiple and irreducible patterns of transitions between response and non-response states. Following this intuition, we here introduce PCI ST , an index that combines dimensionality reduction and a novel metric of recurrence quantification analysis (RQA) to empirically quantify perturbational complexity as the overall number of non-redundant state transitions (ST) caused by the perturbation.
We validated PCI ST on a large dataset of 719 TMS/hd-EEG sessions recorded from 108 healthy subjects during wakefulness, NREM sleep and anesthesia (propofol, midazolam and xenon) as well as 108 brain-injured patients with disorders of consciousness (DOC).
Furthermore, we applied the same index to 84 stereotactic EEG recordings (SEEG) obtained in 9 epileptic patients undergoing intra-cerebral single-pulse electrical stimulation (SPES) during both wakefulness and sleep for clinical evaluation. The index was able to approximate the performance of PCI LZ in healthy subjects and patients, while improving on the previous method by dispensing with both source estimation and permutation statistics, being computed in less than a second; further, its application generalizes from the scalp EEG to sparse intracranial recordings thus potentially enabling the exploration of the complexity of the brain's causal structure across scales, experimental models and species.

Overview of the method
In the original formulation [16], perturbational complexity was calculated by binarizing TMS-evoked potentials (TEPs) at the EEG sources level using a fixed threshold derived from non-parametric statistics with respect to the baseline (pre-stimulus) and subsequently compressing the binary spatiotemporal patterns with the Lempel-Ziv algorithm. An underlying assumption of this binarization strategy is that complex activations engaged by the perturbation appear on the evoked signals as patterns of oscillations around a fixed amplitude scale. Although proven successful when applied to the sources level, this approach is less sensitive in detecting complexity when calculated directly at the EEG-scalp level, where fast oscillations can appear riding on top of larger and slower envelopes as result of volume conduction and signal mixing (S1 Fig). Furthermore, complex neuronal oscillations occurring in amplitude scales that are not determined by a fixed threshold can also be observed in microscopic and mesoscopic recordings due to cross-frequency couplings [30][31][32][33][34], and a binarized measure applied to such scales would have limited applicability to detect complex physiological activations.
Aiming at a general index of perturbational complexity that can be fast and efficiently calculated directly at the EEG-sensors level, we here took a non-binary, geometrical approach and tried to quantify the spatiotemporal complexity of evoked potentials by exploring multiple amplitude fluctuations present in the response. At the scalp level, TEPs appear as time-varying voltage configurations, with local maxima shifting over time across different EEG sensors [16].
When voltages of all electrodes are put together in a vector space, TEPs can be geometrically represented by a trajectory in the so called EEG space [35]. projecting the data in a lower dimensional space spanned by the first three principal components (PCs) of the evoked response. Contrary to the trajectories observed during NREM sleep, which are either local and short-lasting (Fig 1B) or widespread and non-recurrent (Fig 1C), complex TEPs recorded during wakefulness appeared in the principal components' space as entangled, interlaced paths ( Fig 1A). In this geometrical approach, estimating perturbational complexity becomes a matter of quantifying the complexity of these trajectories.
The state transition based perturbational complexity index (PCI ST ) estimates complexity by quantifying transitions between states within the trajectory spanned by the system's response to a perturbation. In the following, we summarize the steps involved in calculating PCI ST (see Materials and Methods for further details). Starting from trial-averaged signals recorded in response to a perturbation, principal component analysis (PCA) of the response is performed in order to effectively reduce the dimension of the data (Fig 2A). The complexity of each resulting principal component is then evaluated by quantifying what we call state transitions, using a method derived from RQA [36,37]. In order to do so, distance matrices, defined by the voltageamplitude distances between all time-points of the signal, are calculated separately for pre-stimulus and post-stimulus samples (Fig 2B). These 2D distance matrices provide a summary of the geometrical relationships of the baseline and the response signals. By thresholding both distances matrices (baseline and response, Fig 2C) at a given scale ε, we obtain a contour plot, called the transition matrix, that depicts the temporal transitions between states -roughly, the ups and downs in the signal -, for both the baseline and the response (Fig 2D). By varying the threshold ε and comparing the number of state transitions (NST) in the matrices of the response with that of the baseline, we can look for the scales at which there are state transitions in the signal's response over and above transitions present in the baseline activity ( Fig 2E). The complexity of each component (red dot in Fig 2E) thus becomes the weighted difference of NST (ΔNST) between baseline and response signals calculated at the threshold ε* where this value is maximum. Finally, PCI ST is defined as the sum of these maximized significant state transitions across all principal components of the evoked signal.
By its definition, PCI ST is high when there are multiple linearly independent components in a spatially distributed response (spatial complexity), each one of them contributing with significant amounts of state transitions (temporal complexity). Conversely, PCI ST is expected to be low either if the perturbation evokes a strongly correlated response across different spatial recordings or if the independent components carry few temporal transitions in the response as compared to the baseline.

PCI ST is reliable and fast in benchmark conditions
PCI ST was calculated on a benchmark of 382 TEPs obtained in a group of 108 subjects during conscious (alert wakefulness) and unconscious (NREM sleep and anesthesia) conditions ( Fig 3A). The wakefulness group presented significantly higher and more variable PCI ST values (mean ± SD, 47.89 ± 12.65) than the NREM sleep/anesthesia group (mean ± SD, 14.19 ± 5.26, p = 4.7x10 -40 , Wilcoxon-ranksum test). In terms of classification performance between conscious and unconscious conditions, PCI ST showed a high classification power that was equivalent to the performance of the original version of PCI in this same dataset [17] (area under the curve (AUC) of the receiver-operator characteristics (ROC) for PCI ST =0.998, AUC for PCI LZ = 0.995).
Indeed, when values for each TMS/hd-EEG session were compared, we found a significant linear correlation between the metrics (r = 0.82, p<10 -95 , Fig 3B). On the other hand, because PCI ST estimates perturbational complexity without employing source localization and surrogate techniques, a key difference found was that PCI ST computation on the whole benchmark dataset was approximately 380 times faster than with PCI LZ . While PCI LZ took about 300 seconds per session to compute (mean ± std, 270s ± 99), PCI ST was calculated in less than one second (mean ± std, 0.71 ± 0.20, p< 10 -127 ) ( Fig 3C).

PCI ST captures spatiotemporal complexity as the joint presence of spatial independence and temporal complexity
PCI ST is the product between the number of selected principal components (N C ) of the evoked response and the average number of state transitions across all its components ( NST ).
The first quantity can be regarded as a measure of the underlying dimensionality of the data and is calculated as the number of linearly independent components that span most of the evoked response in terms of its square amplitude and signal-to-noise ratio. In the case of the TMS/hd-EEG signals, N C is high when different EEG sensors record different evoked potentials and, therefore, it can be viewed as a measure of the spatial differentiation of the brain's response to the perturbation. The second quantity, NST , corresponds to the average temporal complexity present in the individual principal components as measured by the quantification of state transitions. Thus, NST is high when the signal's response displays a variety of patterns of activations across time over and above the patterns occurring before the stimulation onset.
At the group level, TEPs recorded during alert wakefulness presented both significantly higher number of principal components N C (mean ± SD, 7. AUC=0.998, FDR = 6.07 , Fig 3A).
Each nth principal component contributes to PCI ST with its correspondent ΔNST n value.

PCI ST allows a simple and fast set-up at the bedside of patients
After ascertaining the effectiveness of PCI ST to discriminate consciousness from unconsciousness in healthy subjects while dispensing both source estimation and hd-EEG, we tested the performance of the novel index in 108 brain-injured patients with disorders of consciousness (DOC). For each EEG setup (60, 19 and 8 channels, Figure 5C When calculated on the original high-density setup, the sensitivity of PCI ST in detecting signs of consciousness in brain-injured patients was comparable to PCI LZ [17] (Fig 6A, top): In UWS patients, the absence of behavioral signs of consciousness per se cannot be considered a proof of the absence of consciousness. Thus, brain-based measures that do not require subject's interaction with the external environment can be useful to detect a covert capacity for consciousness. In a previous study, PCI LZ detected conscious-like complexity in 20.9% (9/43) of UWS patients, who also had a higher chance of recovery at 6 months [17]. Here, we evaluated whether these patients could also be identified by PCI ST . The novel index calculated on both high-density and standard 10-20 EEG setups detected all (n=9) the patients with high PCI LZ , whereas more than 82% of patients classified as low-complexity by PCI LZ were also below threshold for PCI ST (hd-EEG: 88.2%, 10-20 setup: 82.3%, Fig 6B-top and middle).

PCI ST can be generalized to intracranial responses
Given that PCI ST estimates perturbational complexity directly at the level of sensors, its extension to other types of recordings is straightforward. As a proof of concept, we calculated PCI ST on sparse intracranial SEEG recordings in 9 patients undergoing presurgical mapping for epilepsy. In each patient, SEEG responses evoked by single-pulse electrical stimulation (SPES) of different contacts were recorded both during wakefulness and NREM sleep, resulting in a total of 84 sessions (see S1 Table).

Discussion
In this paper we have developed and tested a new method of estimating pertubational complexity based on dimensionality reduction and state-transitions quantification. The novel index, called PCI ST , represents a fundamental advancement towards the implementation of a reliable and fast clinical tool as well as a general measure than can be employed to explore the neuronal mechanisms of loss/recovery of brain complexity across scales and models.
PCI combines a "perturb and measure" approach with a complexity metric in order to empirically estimate the amount of irreducible information generated through causal interactions within a system [16,39]. This approach is unique because it permits, at least in principle, to explore the structure of causal interactions from the intrinsic perspective of the system under study, above and beyond the patterns of correlations that can be appreciated by a purely extrinsic, observational perspective [5,40]. Further, by measuring the average response evoked by repeated stimulations of the same cortical target, only the effects causally related to the perturbation are amplified, whereas unrelated patterns of activity are averaged out. In this way, it is possible to explore effective interactions while minimizing the effect of noise, common drivers and elements that are causally segregated. Finally, to estimate the complexity of intrinsic causal interactions, only the irreducible (i.e. non-redundant) patterns of spatiotemporal significant activity are considered. Perturbational complexity is both spatial and temporal in nature inasmuch as it is only produced when different elements interact causally in an integrated and differentiated manner. In practice, perturbational complexity is high only when a perturbation propagates reliably over a set of integrated elements that react differently, giving rise to a transient response encompassing a variety of irreducible spatiotemporal patterns of activation.
In the case of PCI ST , this transient response is interpreted as a trajectory in the multidimensional space spanned by the sensors. For example, the TMS-evoked trajectory rests in the sensor space defined by the EEG channels (Fig 1). Mutanen et al. [35] has recently studied the ability of TMS to modulate brain activity by considering single-trial EEG responses to TMS as trajectories in the EEG signal space. Here we extended this geometrical approach and looked at the trajectories spanned by the evoked (trial-averaged) potentials in the rotated basis defined by the principal components of the response (Fig 1). PCI ST quantifies perturbational complexity by estimating the spatiotemporal complexity of these evoked trajectories.
As a first step in this procedure, singular value decomposition (SVD) is used to rotate the trajectories into an orthonormal basis and optimally project them onto a space of lower dimension. This approach is similar to the concept of spatial PCA [41,42], a data-driven technique that has been widely applied to evoked potentials [43][44][45][46][47][48][49] and that can be used to estimate the dimensionality of the trajectories in the EEG space in terms of the number of principal components accounting for the variance of the evoked response [50]. However, dimensionality reduction as operationalized in PCI ST differs from spatial PCA in two important ways: (i) PCI ST calculates SVD exclusively from the post-stimulus signals; (ii) the projection is based not on the covariance but on the correlation across sensors. This procedure guarantees that the principal components, onto which the evoked trajectories are projected, optimally explain the data in terms of the strength (mean field power) of the responses to the perturbation. As expected, when projected onto the space spanned by the components that accounts for most of the response strength, TEPs recorded from healthy awake subjects, which invariably present high value of perturbational complexity as assessed by PCI LZ , were found to have significantly higher-dimensional trajectories as compared to responses obtained in unconscious conditions (Fig 3).
While the number of selected principal components captures the spatial linear independence present in the response, the estimation of spatiotemporal complexity must also take into account the temporal information content present in each component. Since the complexity metric employed in the Lempel-Ziv formulation of perturbational complexity was based on binarizing temporal patterns with respect to a fixed significance threshold, PCI LZ shows a reduced performance when calculated in the absence of source modeling, such as in scalp EEG, when fluctuations of different frequencies relying on distinct neural mechanisms may appear linked due to volume conduction. In this case, complex patterns riding on top of larger and slower components may be easily missed (see S1 Fig). Crucially, sensitivity to cross-frequency patterns of fast oscillations that are modulated by slower envelopes are particularly important because they have been reported at multiple levels, from microscopic to macroscopic recordings [30][31][32][33][34]51], and because they may serve as the mechanism to coordinate neural dynamics and transfer information across spatial and temporal scales [52][53][54][55].
One way to explore the multiple amplitude scales of the activations in a time-series is by calculating its distance matrix [36,56], which provides a summary of the geometrical relationships of the signal, not only in respect to the baseline but to every time point of the signal. Here, we proposed to quantify temporal complexity by binarizing not the original timeseries but its corresponding distance matrix: a technique known as recurrence quantification analysis [36,37]. Recurrences are typically present in complex dynamical systems [57]: eventually every system will return to a state which is sufficiently close to a previous one and quantifying these recurrent states have been shown to be an effective way to detect rhythms and state changes within complex, possibly non-stationary, time-series [58][59][60][61][62][63]. Although inferring recurrences in the "true" state space of a dynamical system from a given time-series depends on properly choosing the embedding parameters to reconstruct the system's dynamical attractors [64] (see Materials and Methods), recurrence techniques are employed in calculating PCI ST simply to quantify patterns in the time-series, without assuming that the underlying attractors are faithfully reconstructed. Indeed, it has been shown that, when applied to experimental data, recurrence plots and the measures derived from them are able to quantify patterns in the dynamics even if no embedding is performed [65]. In the case of PCI ST , independence on embedding dimension has been substantiated by the fact that the performance of the method in discriminating between brain states was found to be stable across different embedding parameters (see S4 Fig).
In the context of RQA, several metrics were proposed to quantify the complexity of the dynamics [36], such as the diagonal-line based entropy [63], the recurrence probability density entropy [66] or the Kolmogorov-Sinai entropy [67]. Although detecting recurrent states does not require signals to be stationary, these metrics are based on statistical properties of densities of recurrent states and may be unreliable in estimating complexity of very short and highly nonstationary time series, such as those evoked by a perturbation. We here took a different perspective by simply counting the number of transitions between recurrent and non-recurrent states (NST) in the principal components space as a measure of the information content of the response to a perturbation. Notably, as this approach might be sensitive to patterns with a regular temporal structure that do not contribute to complexity as measured by compression-based metrics, PCI ST can be expected to accurately estimate complexity only when applied to transient short-lasting signals, such as evoked potentials, in which regular patterns repeating over a long time are unlikely to occur. Future studies in simulated dynamical systems [68,69] and brain network models [2,70,71] should investigate to which extent counting state transitions with PCI ST is enough to capture the complexity of the causal structure of a system. Empirically, PCI ST performed well in the case of the brain, which generates large-scale complex activations in response to a local perturbation in the form of short and transient oscillations with highly non-stationary frequencies. Indeed, PCI ST was strongly correlated to PCI LZ and reliably discriminated consciousness from unconsciousness in the large benchmark population of healthy individuals (Fig 3). Furthermore, and similarly to PCI LZ , the absolute values of PCI ST could be used to compare healthy subjects to brain-injured patients who emerged from coma. When applying the threshold derived from the benchmark population to individual brain-injured patients, PCI ST exhibited high sensitivity in detecting consciousness (Fig 6), with a performance comparable to PCI LZ . This performance was achieved by dispensing with source estimation, forward modeling and non-parametric statistics, thus entailing several practical advantages. First and foremost, PCI ST calculation takes less than one second per session (Fig 3), thus being more than three hundred times faster than PCI LZ . Further, PCI ST can afford an accurate estimation of perturbational complexity even on a reduced (standard 10-20) EEG set-up in both the healthy and brain-injured populations (Figs 5 and 6). Taken together, these advances pave the way to the implementation of a practical, fast and potentially online method to be applied at the bedside in the routine clinical setting.
Beside these practical applications, estimating perturbational complexity based on state transitions enables the exploration of the brain's causal structure across different recordings scales, from macroscopic EEG signals, to mesoscopic local field potentials and, in principle, to microscopic multisite electrophysiological/optical recordings. In the present study we investigated this possibility at the mesoscale level by computing PCI ST on intracranial SPESevoked potentials. This novel metrics was able to detect systematic relative differences in brain complexity between wakefulness and sleep despite the substantial differences between the intracranial (SPES/SEEG) and non-invasive (TMS/EEG) recordings. First, the type of perturbation essentially differs from TMS in that bipolar intracranial SPES not only activates superficial fibers but also deep white matter as well as cell bodies [72] and it does so at a scale of millimeters. Second, the spatial distribution of the recording SEEG electrodes is sparse rather than arranged along a regular spatial matrix. Third, the local field potentials recorded by SEEG electrodes reflect the activity of neuronal assemblies than are substantially smaller than the ones sampled by scalp electrodes [73,74]. Finally, SEEG can reach deep structures, such as hippocampus, amygdala and cingulate cortex [75] that cannot be directly accessed with TMS/EEG experiments.
In addition to the differences listed above, SPES presents key advantage over TMS in that the first is not associated with concurrent auditory and/or somatosensory stimulation [76].
Indeed, when not appropriately controlled for, these additional peripheral inputs induced by TMS may add to the EEG response, thus confounding the genuine effects of direct cortical stimulation [77,78]. In this respect, the present intracranial results serve as a definite confirmation of the fundamental interpretation of perturbational complexity, as originally assessed through TMS/hd- Finally, mesoscale assessment of perturbational complexity is in an ideal position to link microscale explorations at the bench to the macroscale measurements performed at the bedside of brain-injured patients. This connection is fundamental as experiments in both cortical slices [27] and unresponsiveness wakefulness syndrome patients [79] suggest that the fundamental mechanism of loss of complexity is the tendency of neurons to enter a silent period upon an initial activation (OFF-period).
In conclusion, PCI ST , a novel brain complexity index based on dimensionality reduction and the quantification of state transitions, may not only provide a reliable, fast and potentially online option for the assessment of consciousness in the clinical setting, but also serve as a general translational tool for exploring the mechanisms of loss and recovery of brain complexity across species, scales, and models.

Participants
Healthy subjects. The benchmark dataset consisted of 382 TMS/hd-EEG sessions reported in previous works [16,17,79] . Data were recorded from 108 healthy subjects (female, n = 63; age range = 18-80 years) in two conditions: (1) while they were unresponsive and did not provide any subjective report upon awakening (NREM sleep, n = 19; midazolam sedation at anesthetic concentrations, n = 6; anesthesia with xenon, n = 6; anesthesia with propofol, n = 6) and (2) while they were awake and able to provide an immediate subjective report (n = 103, including 32 subjects also recorded in the previously described unresponsive conditions). The experimental protocols applied during sleep and anesthesia have been detailed elsewhere [16,17]. In short, healthy volunteers with history or presence of major medical/neurological disorders and of drug/alcohol abuse were excluded. All participants underwent neurological screening to exclude those at risk of potential adverse effects of TMS. Protocols and informed consents were approved by the local ethical committees (Ospedale "L. Sacco" in Milan, Italy; University of Wisconsin, Madison; Medical School of the University of Liège, Belgium).
Brain-injured patients. TMS/hd-EEG data were also obtained in a population of 108 brain-injured patients (95 reported in a previous work [17]) with newly acquired data recorded following the same previously reported protocol [17]. Besides screening for potential adverse effects of TMS, exclusion criteria were medical instability, refractory generalized seizures, and history of neurodegenerative or psychiatric disease. Sixteen brain-injured patients were Epileptic patients. Data included in the present study derived from a dataset collected during the pre-surgical evaluation of nine (eight previously reported [22]) neurosurgical patients with a history of drug-resistant, focal epilepsy. All subjects were candidates for surgical removal of the epileptic focus. During the pre-surgical evaluation all patients underwent individual investigation with SPES and simultaneous SEEG recordings for mapping eloquent areas and for precisely identifying the epileptogenic cortical network [75]. The investigated hemisphere, the duration of implantation and the location and number of stimulation sites were determined based on the non-invasive clinical assessment. The stimulation, recording and data treatment procedures were approved by the local ethical committee (Niguarda Hospital, Milan, Italy). All patients provided written informed consent.

TMS-EEG measurements and data analysis
Specific protocols for acquiring TMS/hd-EEG potentials were described in [16] and [17].
In brief, data were recorded with a 60-channel TMS-compatible EEG amplifier and MRI-guided TMS pulses were delivered with a focal biphasic stimulator. In each subject, multiple sessions of 200 stimuli were collected with TMS targeted to different areas at different intensities accordingly to the specific protocol. Precision and reproducibility of stimulation were controlled using a Navigated Brain Stimulation (NBS) system (Nexstim Ltd., Finland). The maximum electrical field at the cortex varied within the range of 80 to 160 V/m and targets were selected based on the individual anatomical MRI from within the middle-caudal portion of the superior frontal gyrus, superior parietal gyrus, superior occipital gyrus, and the midline sensorimotor cortex. In brain-injured patients, the stimulation of targets affected by cortical lesions identified on individual magnetic resonance images was deliberately avoided because the EEG response of these areas may be absent or unreliable.
Data was analyzed as previously described [17]. EEG responses to TMS were visually inspected to reject single trials and channels with bad signal quality. Independent component analysis was applied to remove residual artifactual components resulting from eye movements and muscle activations. Bad channels were then interpolated using spherical interpolation and data were bandpass filtered (0.1-45Hz), downsampled to 725 Hz, segmented between -400 and 400ms, re-referenced to the average, baseline corrected (-400 to -5 ms) and averaged across trials.

Intra-cerebral measurements and data analysis
The procedures for SEEG data acquisition are described in [22]. Briefly, intracerebral activity was recorded from platinum-iridium semiflexible multi-contact intracerebral electrodes, with a diameter of 0.8 mm, a contact length of 2 mm, an inter-contact distance of 1.5 mm and a maximum of 18 contacts per electrode (Dixi Medical, Besancon France - Fig 1A-B). SEEG signals were recorded using a 192-channel recording system (NIHON-KOHDEN NEUROFAX-110) with a sampling rate of 1000 Hz. Data were recorded and exported in EEG Nihon-Kohden format. Recordings were referenced to a contact located entirely in the white matter. Electrical stimulation was applied through one pair of adjacent contacts (amplitude 5 mA, frequency=1 Hz, # of stimuli per session=30), while SEEG activity was simultaneously recorded from all other bipolar contacts. The stimulation, recording and data treatment procedures were approved by the local Ethical Committee (protocol number: ID 939, Niguarda Hospital, Milan, Italy). All patients provided written informed consent. SEEG data recorded during both wakefulness and NREM were imported from EEG Nihon Kohden format into Matlab and converted using a customized Matlab-based script. Data were subjected to linear detrend and bandpass filtering (0.5 -300 Hz), using a third order Butterworth filter. Bipolar montages were calculated by subtracting the signals from adjacent contacts of the same depth-electrode to minimize common electrical noise and to maximize spatial resolution [80,81]. Single trials were obtained by using a digital trigger simultaneous to each SPES delivery. Stimulation artifact was reduced by applying a Tukeywindowed median filtering, as in [82], between -5 and 5 ms. Finally, trials and contacts showing pathological activity [83] were detected by visual inspection excluded from the analysis.

Perturbational complexity index based on state transitions
The state transition based perturbational complexity index (PCI ST ) is composed of two steps: dimensionality reduction and transition quantification. In the first step, singular value decomposition (SVD) of the response is performed in order to effectively reduce the dimension of the data. More specifically, let X be the N T  matrix containing trial-averaged signals Then let a n be the normalized eigenvectors of the autocorrelation matrix A calculated from the response signal, In the new orthogonal basis,   n a , X can be represented by a rotated trajectory ( ) t Y : The N T  matrix Y is thus composed by projecting the data X in the space spanned by the principal components (PC) obtained from SVD of the response signal. The PCs are then selected in terms of the eigenvalues of A sorted in decreasing order, so as to account for at least 99% of the response strength measured in terms of the square mean field power, Next, components with low signal-to-noise ratio (SNR ≤ 1.1) are removed, and the remaining N C responsive principal components are used for transition quantification in the following (Fig 2A).
As the resulting components are linearly independent, we regard each component ( ) where L and  are the dimension and the delay of the embedding, respectively and Given an embedded time-series y n (t), its distance matrix D n is then defined as the Euclidean distance between all timepoints of the signal ( Fig 2B): Both matrices are then thresholded at several scales ( Fig 2C): for each threshold  we calculate the transition matrix defined as the contour plot of the correspondent distance matrix at the scale  ( Fig 2D): 1, ( , ) ( , , ) 0, The transition matrix can be thought as the unidirectional derivative of the recurrence matrix R n , also known as recurrence plot [56]: Transition matrices can thus be interpreted as measuring not recurrent states, but transitions between states -between recurrent states, R n (,t i ,t j )=1, and non-recurrent states, R n (,t i ,t j )=0.
Next, we simply cumulate the transition matrix across T C time samples and calculate the mean number of state transitions (NST) for a component n at a given threshold  as The parameter k is introduced to control the relative weight between pre and poststimulus state transitions. By setting the value of k, the crucial scale  n * becomes a data-driven parameter determined by the pre-post comparison for each principal component. Therefore, the temporal complexity of the nth component (ΔNST n ) is the maximized weighted difference at the optimal scale  n * scaled by the number of response samples (T R ) (Fig 2E). The scaling factor yields a quantity that is extensive with the length of the signal's response and largely independent on the sampling rate ( PCI ST was calculated on TMS/EEG signals using three different EEG setups: hd-EEG (all 60 channels), 10-20 system (channels numbered 1,2,3,7,9,11,13,19,25,27,29,31,33,45,47,49,51,53,57,59) and a reduced setup with 8 channels (numbered 9,11,27,29,31,49,57,59). For each setup, signals were referenced to the corresponding common average reference.

Statistical analysis
Statistical analysis was performed with IBM SPSS v. 22 A linear support vector classifier [38], independently trained on all TMS/EEG sessions of the benchmark dataset, was employed to predict the level of consciousness of brain-injured patients who had recovered from coma (see S2 Fig