Nonlinear dynamics captures brain states at different levels of consciousness in patients anesthetized with propofol

Summary

The information processing capability of the brain decreases during unconscious states. Capturing this decrease during anesthesia-induced unconsciousness has been attempted using standard spectral analyses as these correlate relatively well with breakdowns in corticothalamic networks. Much of this work has involved the use of propofol to perturb brain activity, as it is one of the most widely used anesthetics for routine surgical anesthesia. Propofol administration alone produces EEG spectral characteristics similar to most hypnotics; however, inter-individual and drug variation render spectral measures inconsistent. Complexity measures of EEG signals could offer better measures to distinguish brain states, because brain activity exhibits nonlinear behavior at several scales during transitions of consciousness. We tested the potential of complexity analyses from nonlinear dynamics to identify loss and recovery of consciousness at clinically relevant timepoints. Patients undergoing propofol general anesthesia for various surgical procedures were identified as having changes in states of consciousness by the loss and recovery of response to verbal stimuli after induction and upon cessation of anesthesia, respectively. We demonstrate that nonlinear dynamics analyses showed more significant differences between consciousness states than spectral measures. Notably, attractors in conscious and anesthesia-induced unconscious states exhibited significantly different shapes. These shapes have implications for network connectivity, information processing, and the total number of states available to the brain at these different levels. They also reflect some of our general understanding of the network effects of consciousness in a way that spectral measures cannot. Thus, complexity measures could provide a universal means for reliably capturing depth of consciousness based on EEG changes at the beginning and end of anesthesia administration.

Keywords:

Surgical and invasive medical procedures – Electroencephalography – Nonlinear dynamics – Consciousness – Ellipses – Anesthesia – Anesthetics – Syncope

Introduction

The brain performs computations in both linear and nonlinear fashions, and both cognitive and arousal states have been characterized using analyses based on nonlinear dynamics [1–5]. Measures capturing the complexity of EEG signals from nonlinear dynamics positively correlate with changes in brain functional connectivity [6] and thus can be an indicator of the brain’s ability to process information. In particular, electroencephalogram (EEG) signals derived from the frontal cortex show stereotypic responses to GABAergic anesthetics (e.g. propofol and sevoflurane). EEG frequency domain measures are particularly sensitive to anesthetic levels [7–12] and anesthetics produce frequency slowing, partially due to prolonged GABA-mediated inhibition [13–15]. For most anesthetics, loss of response (equated with loss of consciousness) in humans occurs when low amplitude, higher frequency EEG waveforms are replaced by higher amplitude slow-wave patterns, similar to the delta rhythms seen during slow wave sleep [10,16]. These delta rhythms are gradually replaced by burst suppression patterns as patients transition to deeper surgical planes of anesthesia [17]. Interestingly, recovery of consciousness is associated with much less delta activity and higher amounts of low amplitude, fast activity on awakening, compared to loss of consciousness, indicating a hysteresis of brain states for loss and recovery from anesthesia [16,18], since a more wakeful looking pattern of higher frequency low amplitude activity is seen on awakening. This has also been termed ‘neural inertia’ and appears to reflect differential activation of neural networks during these state transitions [19,20]. Emergence from general anesthesia also does not typically follow a common pattern, instead a more heterogenous spectral pattern appears during re-establishment of conscious awareness and responsivity [16]. Importantly, frequency domain measures can differ for various classes of anesthetics and adjuvant agents (GABA versus non-GABAergic), in response to surgical stimulation, across age groups, and between people, making spectral domain measurements non-invariant.

Anesthesiologists have noted that nonlinear dynamic measures of EEG signals could distinguish anesthetic depth [1,3] as well as recovery of consciousness in surgical patients [2]. One such method is the creation of time-delayed embeddings from EEG signals. Sometimes, time-delayed embedding signals are observed to explore only a limited region of all possible states in a multidimensional space in which they are embedded; they settle onto attractors in state space, and can thus be seen as probes of the underlying dynamics. We have previously shown that attractors generated from microEEG signals are sensitive to isoflurane-induced loss of righting responses (surrogate loss of consciousness measure; [21] in rodents [22]. Specifically, we observed that initially spherical attractors became flattened and more ellipsoidal when moving from waking to loss of righting reflex with both anesthesia administration and during sleep [22]. We observed a similar attractor shape change when characterizing EEGs from human subjects anesthetized with a combination of remifentanil and nitrous oxide [23], and another geriatric patient population anesthetized with propofol and fentanyl [24]. We quantified the shape changes we observed using a geometric phase-space analysis; we fit the 3D attractor with an ellipsoidal solid of revolution [25] and then calculated the ratio between the major and minor axes. This ellipse radius ratio (ERR) changes significantly and consistently before and after loss and recovery of response during several anesthetic regimes [23,24]. Additionally, the ERR also shows a significant difference in adult patients across a wide range of ages [23,24]. Characterization of attractors using correlation dimension provides a measure of complexity or information content of EEG signals [2,3]. Currently, it is unknown whether correlation dimension or our ERR analysis is sensitive to recovery of consciousness for propofol anesthesia. It is also unknown how loss and recovery of response properties of EEG signals differ, and how the dynamics change with anesthetic administration and time as patients lose and regain consciousness. Understanding the dynamical complexity changes that occur during loss and recovery of consciousness would support theories demonstrating breakdowns in connectivity and information processing during unconsciousness [26–29] and enhance anesthetic depth monitoring.

To determine whether complexity of brain activity changes with depth of consciousness, we tested the ability of correlation dimension and our ERR analysis to distinguish between EEG activity before and after the loss and recovery of response during propofol anesthesia in surgical patients. We compared our analyses from nonlinear dynamics to frequency-derived measures, and characterize the dynamics that occur during these transitions.

Materials and methods

Patient selection

A retrospective analysis of EEG data as part of a medical chart review was approved by the Stanford University School of Medicine Institutional Review Board Protocol 28130. Twenty patients were chosen from a population of surgical all-comers to the operating room at Stanford Hospital and Clinics from March 1^st, 2015 until July 1^st, 2016. Patients were included if they were to be anesthetized with propofol as their maintenance anesthetic (a total intravenous anesthetic or TIVA). Subjects were between the ages of 23 and 87 years old (median age 51), comprised American Society of Anesthesiology (ASA) classes 1–3, and were undergoing various surgical procedures (breast, gynecological, head and neck, and neurosurgery). Subjects included 17 females and 3 males. We deliberately did not restrict the population studied to a single surgical cohort, in order to better reflect the diversity of patients that would undergo surgical anesthesia, and to demonstrate that nonlinear dynamics tools would be sensitive to brain state.

Anesthetic administration

Patients were administered an infusion of intravenous propofol over minutes in order to achieve a slowly increasing plasma concentration of anesthetic for hemodynamic stability on induction of anesthesia, and to more accurately determine the precise time of the loss of response (LOR). LOR was defined at the moment the patient lost response to verbal stimulus (eye opening or verbal response to name), assessed every 5 seconds after the infusion was started. If there was doubt regarding the LOR, an additional verbal stimulus was delivered to the patient. Recovery of response (ROR) was similarly assessed by verbal stimulus (stating the patient's name and asking them to open their eyes) after the maintenance anesthetic was turned off. Timing of LOR and ROR were marked to the nearest second for off-line analysis of EEG signals. We estimate an error range for the observer to be less than 5 seconds based on continuous vigilance of the patient between stimuli. Raw EEG data, processed EEG data (patient state index or PSI), and time stamps were downloaded for off-line data analysis.

Choice of airway device (laryngeal mask airway or endotracheal tube) for general anesthesia was determined by the anesthesiologist based on the type of surgical procedure and patient needs. LOR measurement was unaffected since it occurred prior to instrumentation of the airway. Patients were maintained on propofol infusions as the maintenance anesthetic (TIVA), with analgesic adjuncts as needed according to the anesthesiologist’s clinical judgment. Some patients received pre-operative medication (midazolam and/or fentanyl) prior to initiation of the propofol infusion as clinically necessary. Intraoperative analgesia was supplied by fentanyl, hydromorphone, ketorolac, acetaminophen, or ketamine IV boluses, and supplemented by lidocaine, fentanyl, or remifentanil infusions according to the anesthesiologist's discretion and the degree of surgical stimulation. Timing of ROR would have been influenced by maintenance anesthesia and analgesic adjuncts, but was equally compared across all analyses.

EEG data acquisition and preprocessing

Standard American Society of Anesthesiologists (ASA) monitors (pulse oximeter, non-invasive blood pressure cuff, capnograph, electrocardiogram) [30] were utilized for each patient in order to measure intraoperative oxygen saturation, hemodynamic parameters, and end-tidal gas concentrations (O₂, CO₂). A self-adhesive, 5-lead frontal EEG electrode array was placed according to the manufacturer's instructions (SedLine Legacy, Masimo, Irvine, CA) on the forehead of each patient prior to induction of anesthesia. Electrodes were applied at approximately F7, F8, Fp1, and Fp2 referenced to AFz (Fig 1A) in the standard 10–20 electrode EEG montage. Appropriate electrode impedance was checked according to an automated SedLine routine that sends a ~78 Hz impedance pulse to each electrode. Data was digitized at 250 Hz. EEG recordings from F7 were de-trended and then filtered using a 50 Hz low-pass Butterworth Infinite Impulse Response (IIR) filter, and a 1 Hz high-pass Butterworth IIR filter for subsequent analysis. Note the 50 Hz low-pass filter does not remove dynamics at short timescale (i.e. that below 50Hz). Dynamics remaining at short time-scales stem from phase differences between lower frequency components.

Identification of clips for analysis has been previously described[23,24]. Briefly, 20 second artifact-free clips occurring within 2 minutes before to 2 min after loss and recovery of response to verbal commands (LOR and ROR respectively) were selected manually for subsequent analysis (Fig 1A). Two of the authors visually inspected the EEG traces as well as the EEG spectrograms to ensure the clips did not contain burst suppression or artifacts (SLE, MBM). From the original 20 patients, 19 were selected for LOR analysis and 16 selected for ROR analysis given these were the patients who had artifact-free EEGs for at least 20 continuous seconds before and after the events of interest.

Spectral analysis

EEG Data was analyzed using custom scripts written in MATLAB (Mathworks, Natick, Massachusetts). Power spectral density was calculated on the 20 second EEG clips during the pre- and post-LOR and pre- and post-ROR periods using the multitaper method as part of the ‘mtspectrumc’ function of the Chronux toolbox [31] (http://chronux.org/). Specifically, we used a time-bandwidth product of 5, 9 tapers, and limited the frequency range to 0 to 50 Hz. The error range for power spectral density plots was computed using the theoretical estimate method at 0.05 significance level. Power values were expressed in decibels.

To visualize the temporal profiles of the spectral changes occurring during the 4 minute windows surrounding the LOR and ROR transitions (i.e. 2 minutes before to 2 minutes after), we computed normalized spectrograms with custom MATLAB scripts (Fig 1B). A Fourier transform (using ‘fft’ function in MATLAB) was performed using Hann windows with half window overlap. The magnitude was converted to decibels (dB) and the spectrogram was normalized by its maximum magnitude.

Total power and spectral edge frequency (SEF, the frequency below which 95% of the total EEG power resides) were also calculated using multitaper spectral analysis [31] on the filtered EEG. In addition, we calculated the percentage of total power for individual frequency bands per condition (pre- and post-LOR, pre- and post-ROR). The percentage of total power was used because significant changes in EEG total power with exposure to propofol have been previously reported. The frequency band ranges were defined as follows: delta: 1 to 4 Hz; theta: 4 to 8 Hz; alpha: 8 to 12 Hz; beta: 12 to 25 Hz; and gamma: 25 to 50 Hz [32,33].

Characterization of dynamical attractors

Geometric phase-space analysis

We constructed three-dimensional, time-delayed embeddings (attractors) of the EEG signal before and after LOR and ROR using a 4 ms delay for correlation dimension analysis (Fig 1C). We chose this delay as it was the smallest delay possible for our sampling frequency and showed the most significant effects based on previous work [23,24]. We observed shape changes in attractors before and after our time points of interest when plotted at this timescale. We quantified this shape change by fitting the three-dimensional attractor to an ellipsoidal solid of revolution [23–25] (Fig 1D). The lengths of the symmetry axes of the ellipsoid were calculated and the ratio of the minimum and maximum axes (which we term the ellipsoid radius ratio, ERR) was used to quantify the shape change. A radius ratio of 1 implies a sphere, while smaller ratios imply more strongly ellipsoidal shapes.

Correlation dimension

Using the same three-dimensional attractor with 4 ms embedding delay (Fig 1C), we tested whether the correlation dimension captured the EEG changes that occurred before and after LOR and ROR. We used a previously reported and commonly used algorithm to compute the non-integer (fractal) correlation dimension that works for irregularly sampled objects (e.g., a point cloud in our case) [1,2,34]. We also tested whether significant differences in correlation dimension could be observed when we increased the dimensionality of the embeddings from 3 to 5 dimensions.

Specifically, we use the classic Grassberger and Procaccia method to estimate the correlation dimension [34]. Grassberger and Procaccia defined the correlation integral as:

They demonstrated that, for small r, C(r) behaves as a power law in r [34], so that

Complexity measures by embedding delay

We tested also whether the ellipse radius ratio and correlation dimension were changed by varying the embedding delay time. To estimate the optimal delay for creating the attractors we calculated the first zero-crossing of the autocorrelations of the pre-LOR and post-LOR signals. We used this value to set our maximum range, and tested multiple embedding delays (4, 8, 12, 52, 100, 500, 1000, 1500, 2000, and 2500 ms) between the shortest time window (shifting the EEG by 1 point) to the largest (set by the autocorrelation zero crossing). Spearman’s correlation was used to determine whether a trend existed between the embedding delays and ellipse radius ratio or correlation dimension. We tested whether a trend existed for both 3D and 5D time-delayed embeddings and correlation dimension.

Dynamics of complexity measures

To explore the dynamics of our measures around our two clinically relevant time periods (LOR and ROR) we calculated ERR and correlation dimension in a 4 min window surrounding these events. Specifically, 20 second windows of EEG activity starting 2 min before an event to 2 min after the event were used to create attractors. We shifted these 20 second windows every 5 seconds so windows had 75% overlap. We created the attractors using the smallest embedding delay (4 ms) as this delay showed the most significant differences before and after LOR and ROR. We then calculated the ERR and correlation dimension (both 3D and 5D) for each 20 second segment of EEG activity. The average and standard error of the mean of all EEG data from 19 LOR and 16 ROR patients were calculated to show the dynamics of our measures across this period.

Effect size and correlations of EEG measures

We calculated a paired-data Cohen’s D on our EEG measures before and after LOR and before and after ROR for the complexity (ERR and correlation dimension) measures. We also calculated the Cohen’s D for the spectral edge frequency, total power, and percentage of power in each of the frequency bands (delta, theta, alpha, beta, and gamma) for comparison. We used the 4 ms embedding delay complexity measures as these showed the most significant results. Additionally, we calculated Spearman correlations between the ERR, 3D, and 5D correlation dimension measures and age to determine if these patient characteristics influenced our results. We also calculated the Spearman correlation between the ERR and 3D correlation dimension (again at the shortest 4 ms delay) with each other and with the percentage of power in the band-limited frequency ranges. We calculated these correlations independently for LOR and ROR periods.

Since we had a large age range in our study, we tested whether patient age may have influenced our measures. To determine whether complexity measures change with age we calculated the Spearman correlation between age and change in pre- to post-ERR and 3D CD (at 4 ms delay) measures individually for LOR and ROR. There were no significant correlations in the change in ERR or 3D CD and age for LOR (ERR and age r = 0.17, p = 0.48; 3D CD and age r = -0.21, p = 0.38, p values are uncorrected) and ROR (ERR and age r = 0.006, p = 0.98; 3D CD and age r = 0.04, p = 0.88, p values are uncorrected). Given the lack of correlation in these measures we reasoned that our complexity EEG measures were not confounded by age.

Statistics

Significance values for multiple comparisons within a particular analysis type were corrected using the Holm-Bonferroni method—a sequentially-rejective procedure [35]. Specifically, we corrected p values for pre- vs post-metrics for LOR and ROR separately within each of the following analyses: (1) the power percentage across all 5 frequency bands (delta, theta, alpha, beta, gamma), (2) the ellipse radius ratio (ERR) when comparing across all 10 embedding delays, (3) the correlation dimension when comparing across all 10 embedding delays, and (4) the Spearman correlation of ERR, correlation dimension, and frequency band power. We report our results as medians (25, 75 percentiles), and corrected significance values (p) are calculated from Wilcoxon Signed Rank Tests.

Results

Spectral differences between pre- and post-LOR and ROR transitions are best captured by spectral edge frequency and gamma

We measured the EEG spectrum collected from F7 referenced to AFz in 20 s artifact-free clips identified pre- and post-loss (LOR) and recovery (ROR) of response (Fig 1A). Visualization of the EEG spectrum before and after loss of response (LOR) shows similar changes to what has been previously reported: notably, increases in lower frequencies and decreases in higher frequencies, including a visually apparent increase in the alpha range (8–12 Hz) [16,32,36] (Fig 2). In contrast, emergence from propofol anesthesia is marked by a more uniform distribution of power across frequency bands (Fig 2), which has also been reported previously [16]. When we look at the spectral changes through time we observe the slower dynamic changes in the spectral content throughout LOR (Figs 1B and 3). It appears that frequencies below 12 Hz gradually increase, while frequencies above 12 Hz gradually decrease and the dominant frequency, near the alpha (8–12 Hz) range, emerges near to and following the LOR transition (Figs 1B and 3). However, this alpha increase does not reach statistical significance in this sample (Table 1). In this cohort of patients, the ROR dynamics appear more abrupt, with sharp changes in frequency content near the ROR transition (Fig 3). The more uniform power distribution across frequency bands occurs within seconds of the marked ROR event, and the previously dominant alpha rhythm is replaced by a more pronounced beta rhythm (Fig 3), though it does not reach statistical significance in this sample (Table 2).

We quantified these spectral changes in 20 second, artifact-free clips selected before and after LOR and ROR using spectral edge frequency, total power, and the percentage of power at delta (1-4Hz), theta (4 – 8Hz), alpha (8 – 12Hz), beta (12- 25Hz), and gamma (25 – 50Hz) frequency bands. We observed a significant decrease in spectral edge frequency from pre- to post-LOR (Fig 4A), and a significant increase in spectral edge frequency from pre- to post-ROR (Fig 4B). Additionally, we observed a significant increase in total power from pre- to post-LOR (Fig 4C), but no significant change in total power from before to after ROR (Fig 4D). This is reflected in a more uniform distribution of power across frequency bands (Fig 3), which recapitulates the emergence patterns of other GABAergic agents, like sevoflurane, reported in other studies [16]. We summarize the percentage of power in band limited frequency ranges for LOR in Table 1 and for ROR in Table 2. Our pre- and post-LOR clips showed significant differences in the percentage of gamma power (Table 1). The percentage of power in the alpha range was close to significance but did not pass the correction for multiple comparisons (Table 1). Significant differences in pre- versus post-ROR clips were observed in delta (decrease), theta (decrease), alpha (decrease), beta (increase), and gamma (increase) percentages of power (Table 2). There was also an increase in beta following ROR, but it did not reach significance after being corrected for multiple comparisons (Table 2).

Attractor characterization using ellipse radius ratio and correlation dimension distinguishes between pre- and post-response transitions

To assess complexity differences before and after LOR and ROR we created three-dimensional time-delayed embeddings of the 20 second clips using a 4 ms delay (shown in Figs 1C, 1D and 5 flattened in 2D). These embeddings appear to explore a limited amount of the possible 3D space, suggesting that they reveal attractors in the dynamics. As in previous reports, we observed flatter, more ellipsoidal attractors during periods when responses were absent and patients were anesthetized (Figs 1C, 1D and 5). Conversely, broader, more spheroidal shapes were observed when patients were responsive (Figs 1C, 1D and 5).

To quantify the changes in the attractors observed before and after LOR and ROR, we measured the ellipse radius ratio (ERR) of the 3D attractor point clouds fit with ellipsoidal solids of revolution as previously described [23,24]. Consistent with our previous work, a significant decrease in ERR was observed in the clips preceding and following LOR, while a significant increase in ERR was observed in the clips preceding and following ROR (Fig 6A). The attractor shape change was also quantified using the correlation dimension (CD). A nonsignificant decrease was observed in the correlation dimension computed using a three-dimensional embedding (3D CD) following LOR; however, a significant increase in 3D CD from before to after ROR was observed (Fig 6B). To determine if this effect occurred if the embedding dimensionality was increased, we used a five-dimensional embedding and recalculated the correlation dimension (5D CD). A significant change in 5D CD was observed in both conditions–a significant decrease pre- and post-LOR, and a significant increase pre- and post-ROR (Fig 6C).

Ellipse radius ratio is embedding delay dependent, but the correlation dimension is not

We explored whether these differences in ERR and correlation dimension were dependent on the embedding delay used to create the attractors. To test this, we characterized attractors using 10 successively increasing delays (4, 8, 12, 52, 100, 500, 1000, 1500, 2000, and 2500 ms, max set using autocorrelation of EEG signals). Visual inspection of the attractors shows that the obvious shape differences seen preceding and following both LOR (Fig 7A) and ROR (Fig 7B) at smaller delays become less pronounced as the delay increases, successively fading as attractors are plotted at larger delays. We assessed the significance of the measured values by delay. The ellipse radius ratio (ERR) was significant for the first three delays (p, corrected = 0.0015, 0.0015, 0.011 for 4, 8, and 12 ms respectively) for loss of response, but not for the rest of the delays (Fig 8A, top). Additionally, the ellipse radius ratio (ERR) was significant for the first three delays (p, corrected = 0.0039, 0.0039, 0.009 for 4, 8, and 12 ms respectively) for recovery of response, but not for the rest of the delays (Fig 8A, bottom). The 3D correlation dimension (3D CD) for loss of response was not significant for any of the delays (Fig 8B, top). However, for recovery of response, the 3D CD was significant for all delays (p, corrected = 0.0013, 0.0045, 0.011, 0.0032, 0.0013, 0.0019, 0.0045, 0.011, 0.0013, 0.0032 for 4, 8, 12, 52, 100, 500, 1000, 1500, 2000, and 2500 ms, respectively, Fig 8B, bottom). The 5D CD was not significant for any of the delays during LOR (Fig 8C, top). Note this result is different from Fig 6C because here we are controlling for multiple comparisons, since we are now comparing 10 embedding delays. The significance value needs to be adjusted in this case. Unlike the 3D CD, the 5D CD was significant for only the first delay (p, corrected = 0.022 for 4 ms) for recovery of response (Fig 8C, bottom).

To reveal the scale of the dynamics of our measures and to determine whether there was a relationship between the embedding delay and the difference in our complexity measures, the Spearman correlation was calculated between the 10 embedding delays we used and the ellipse radius ratio (ERR), 3D, and 5D correlation dimensions. There were no significant correlations between the delays and ERR for LOR (r = 0.24, p = 0.51) or ROR (r = -0.66, p = 0.04). Additionally, there were no significant correlations between embedding delay and 3D correlation dimension for LOR (r = 0.62, p = 0.06) or ROR (r = -0.42, p = 0.23). The Spearman correlation was not significant (when controlling for multiple comparisons), but was strongest between the embedding delay and the 5D correlation dimension for LOR (r = 0.72, p = 0.02) and ROR (r = -0.70, p = 0.03). Given the lack of significant relationships between embedding delays and our attractor characterizations, p-values reported here have not been corrected for multiple comparisons.

Attractor differences track spectral changes and patient behavioral state

To explore the dynamics of our complexity measures around our events of interest, we evaluated the two complexity measures in the 2-minute epochs preceding and following either LOR or ROR. LOR spectrograms and their corresponding attractors, calculated in 20 second adjacent windows, are shown in Fig 9. One can clearly see the attractors flattening and becoming more ellipsoidal when approaching LOR and continuing to flatten after LOR in both example patients. ERR, 3D, and 5D correlation dimensions are also shown for both patient examples in Fig 9A and 9B. In both individuals, the LOR trend appears to be a decline in these measures across this window, though not consistently. In these patients, ROR shows the opposite trend, with attractors inflating and becoming more spheroidal as a patient becomes responsive (Fig 10A and 10B). This is captured by all three attractor measures, the ERR, 3D, and 5D correlation dimensions, which all increase in the period spanning ROR. There is, however, a difference in the two patient examples shown. The patient in Fig 10A appears to be moving to a more unresponsive state at approximately minute +1.75, indicated by the decrease in higher frequencies and reappearance of a slow-wave dominated spectrogram. The attractors follow this trend during this period and become more ellipsoidal and flat. In the patient in Fig 10B, however, we do not see this reversion, reflecting a behavioral state in which the subject remained responsive. When we explore the dynamics of these measures for all patients, we can see a generally decreasing trend across the LOR time window and an increasing trend across the ROR time window (Fig 11). The LOR dynamics appear to change more gradually and consistently, while the ROR appears to have an abrupt jump in the complexity measures once the patient regains their ability to respond (Fig 11).

Complexity measure effect sizes are comparable to spectral measures and strongly correlated with gamma

To determine the effect size of our EEG measures we calculated pairwise Cohen’s D for spectral edge frequency, total power, band-limited frequency ranges, and complexity measures (Table 3). This measure gives a description of the difference of means, with +/-1 reflecting a large effect size. The ERR phase-space analysis, spectral edge frequency, total power, and gamma power reflected the largest effect sizes for LOR, of approximately 1 (Table 3, middle column). Alpha power, 5D correlation dimension, and 3D correlation dimension had medium effect sizes of 0.6–0.7. For ROR the largest effects sizes were ERR (1.4) and spectral edge frequency (1.3), followed closely by 3D (1.2) and 5D (1.0) correlation dimension estimates (Table 3, last column). While these effect sizes for ROR were greater than LOR, the total power effect size was much smaller in the ROR condition (Table 3). Effect sizes of changes in power in all the frequency bands were relatively larger for ROR, with the exception of the gamma band (Table 3). To determine whether there was a relationship between our EEG measures we calculated the Spearman correlation. Specifically, the Spearman correlation between the pre- and post-event complexity measures (ERR, 3D CD at 4 ms delay) and percentage of power differences in individual frequency bands was computed separately for LOR (Table 4) and ROR (Table 5). During LOR, there was a significant relationship between the change in ERR and gamma, as well as between 3D CD and gamma (Table 4). Similarly, during recovery of response (ROR), significant correlations were observed between ERR and gamma as well as 3D CD and gamma (Table 5). Interestingly, the correlation between ERR and 3D CD reached significance for the ROR condition, even when corrected for multiple comparisons; it did not rise to significance for the LOR condition.

Discussion

Our work supports other neuroscience studies showing that complexity changes occur in neural activity with consciousness transitions [2,3,23,24,26–29]. We demonstrate here that these approaches provide equivalent discrimination between periods preceding and following LOR, and between periods preceding and following ROR as standard spectral measures. As in our previous work, our geometric phase-space ellipse radius ratio (ERR) was significantly different between pre and post-LOR and ROR clips. Additionally, significant differences were found using 3D correlation dimension between pre and post-ROR. We also found that our ERR values were correlated with correlation dimension values during ROR. Since we observed that ERR is correlated with multiscale entropy in a previous study [24], another measure of complexity, and with correlation dimension used here, we reason that ERR may itself capture EEG complexity. In our previous study, the ERR had a greater effect size than multiscale entropy [24]. We expect similar changes across LOR would be observed using other complexity measures, such as permutation or approximate entropy, as these measures decrease with increasing anesthetic depth for GABA-ergic agents [37–39]. In fully awake patients, the attractor is a nearly perfect sphere. As patients lose consciousness, the attractor collapses into an ellipsoidal shape and explores less of the possible trajectory shape, thus demonstrating decreased complexity and information content. After LOR a clear additional compression of the attractor was evident when comparing the 20-second clips of EEG before and after patients stopped responding to verbal command. The opposite trend is observed upon ROR indicating a recovery of complexity and information upon arousal. This finding confirms a previous report of correlation dimension increase with ROR from anesthetic administration [2] and demonstrates for the first time that the dimensionality shape change is more pronounced upon recovery compared to LOR. Additionally, we also observed an abrupt change in the correlation dimension (as well as our other measures) near the ROR transition as has been previously described [2] and is distinct from the gradual decreases observed during LOR (Fig 9). We expect this dynamical difference is due to the gradual induction of anesthetic administration used in our study which through direct observation by the anesthesiologist (DC) caused a gradual loss of response in patients during LOR. Upon ROR, the transition is more abrupt. The ERR changes occur at the ROR timepoint; whereas, the correlation dimension lags the transition by about 20 s. This shows our ERR measures is more sensitive to the changes in EEG around ROR timepoints compared to CD and may weaken the positive predictive value of CD. Additionally, our observations support theories on mechanisms of hysteresis–the asymmetry of anesthetic concentration at loss and recovery of consciousness [20]. Our attractors were noticeably flatter on loss of consciousness compared to recovery of consciousness, indicating a higher degree of drug effect (higher propofol concentrations) are present at loss of consciousness, and lower concentrations are present on recovery. The differences we show in the magnitude and dynamics of complexity around LOR and ROR, support the existence of distinct neural circuits during these two conscious transition timepoints.

In this data set, attractors following ROR were observed to be more spheroidal than prior to LOR (Fig 5). Following this observed difference, the ERR, 3D and 5D CD measures consistently show lower values pre-LOR compared to post-LOR (Fig 6). This is a function of the slow propofol infusions being started significantly before (i.e., minutes) we started calculating the attractors in the LOR condition. Additionally, we can see this in the subtle differences in ERR and 3D CD between pre-LOR compared to pre-ROR. Note there is no difference in these conditions in 5D CD. This is because in both conditions the subject is on the edge of wakefulness. That is, a subject who has been exposed to a propofol infusion for ~2 minutes, before losing all responsiveness might be in exactly the same state of consciousness as when the anesthetic is slowly washing out of his/her brain, just before they regain full wakefulness. We wanted to test the limits of these measures, so we intentionally compared these close time windows to anesthetic depth transitions (loss of consciousness). The attractor shape changes that occur over LOR (Fig 9) and ROR (Fig 10), shows that the ERR measure closely follows the dynamics during these events. This is especially true in Patient G (Fig 10A), where the patient becomes less responsive with lack of stimulation after initial recovery of consciousness. This relapse to unconsciousness is not uncommon in patients who have emerged from anesthesia but are no longer being stimulated. It is as if their states of wakefulness and consciousness oscillate. Also, many drugs such as propofol and the more lipophilic volatile anesthetics have a context-sensitive half-life [40] and remain in fatty tissue. They form a depot or reservoir of anesthetic that continues to leach out even after the maintenance anesthetic is terminated. Additionally, careful observation of the dynamics of these measures around LOR and ROR timepoints show peaks in our measures shortly following the transitions. These likely reflect intubation and extubation stimulation, respectively, which can arouse patients and would cause increases in our measures. The careful tracking of our method during this period of recovery means these methods could be useful to assess level of consciousness in the recovery post-anesthesia care unit as an added measure of vigilance and safety.

Unlike previous reports [23,24] where complexity measures show enhanced performance over spectral measures, here our complexity measures had equivalent performance compared to spectral measures. We reason that this is due to the use of propofol with adjuvant agents (e.g. remifentanil, fentanyl, midazolam) which enhances the spectral signal change from high frequency, low amplitude awake activity to low frequency, high amplitude activity upon anesthetic induction [8,41–44]. For instance, we have previously shown that SEF, the bispectral index (which is thought to capture a ratio of high frequency and low frequency values), and our ERR measure all show similar changes during LOR and ROR time periods [23]. In contrast, use of agents like nitrous oxide with propofol make spectral measurements less reliable [45–47]. Note, however, in our previous publication the influence of remifentanil actually enhanced the ability of SEF and BIS to detect changes during consciousness transitions at clinically relevant timepoints [23]. In the present study, the predominant anesthetic was propofol; thus, spectral measures were expected to discriminate between the two states well.

Previous reports indicate significant differences in alpha power following induction of anesthesia [32,48,49]. In our previous work, we found significant differences in geriatric patients induced with propofol and fentanyl [24]. However, here we found an increase but not significant difference in the percentage of alpha before and after LOR (Table 1), but did observe a significant reduction upon ROR (Table 2). Our time points are very close to the clinically relevant loss of response and propofol anesthesia was delivered slowly and gradually over several minutes. Whereas, in our previous study [24] propofol was delivered in a bolus, so LOR and EEG changes were likely achieved faster than those observed here. This contrast between anesthesia administration supports proposed mechanisms of propofol-induced unconsciousness, whereby disruption of corticothalamic circuits causes increases in alpha power immediately post induction, which does not continue to increase near LOR [50]. Thus, the dynamic state of the brain around LOR was not so greatly altered during the time points investigated. This highlights the non-inferiority of complexity measures to capture subtle EEG changes compared to standard spectral analyses.

As is the case with any study involving EEG signals, we cannot rule out that some unknown fraction of our signal is likely contaminated by electromyogram (EMG, muscle) activity. In the same manner as our previous studies, EEG segments were visually inspected to remove any obvious artifact events [23,24]. Unfortunately, more subtle EMG contamination can extend as low as 20Hz [51–53]. However, higher frequency EEG components correlate with cognition and are significant contributors to judging anesthetic depth [23,51,54]. The truth for all EEG studies is that EMG signals are present in real-world scenarios. However, measures that can distinguish clinically relevant states even with contamination have practical utility even if they carry limits in their mechanistic interpretations.

Since we are retrospectively analyzing previously collected data, we have some limitations in interpretation. One issue is the unbalanced gender split that we had in this dataset (17 women and only 3 men). In our previous work [23,24], we have had much more even gender splits and did not find any differences in our measurements; thus, we do not expect our gender imbalance to influence our results here. Future work, should include more evenly distributed gender splits to make results more applicable to all patients. Our analyses also appear to distinguish between these states well regardless of age. In addition, we restricted our analysis to patients who were maintained on propofol infusions; however, adjuvant agents that could influence the EEG, including lidocaine, remifentanil, midazolam, and fentanyl, were also co-administered. In humans, slowing of EEG activity from high frequency, low amplitude signals to low frequency, high amplitude signals is observed when administering these agents [8,55–58]. This is reflected by significant observations in increased total power [59], delta power [59], and decreased spectral edge frequency [55]. When combined with propofol for sedation, remifentanil decreases the amount of anesthetic needed and can allow for lighter sedation levels while maintaining patient comfort [60]. Additionally, a similar decrease in SEF was observed after patients lost response to verbal commands (LOR); however, addition of remifentanil (compared to propofol only) attenuated this effect in a dose-dependent manner [61]. It also appears to attenuate the increase in delta activity following LOR compared to the administration of propofol alone [59]. Lidocaine decreases the median frequency when added to propofol anesthesia [62,63]. Co-administration of fentanyl with propofol anesthesia decreases the amount of propofol needed to keep patients unconscious but does not influence spectral measures upon emergence [57]. All of these agents appear to influence the EEG dynamics in the same manner as propofol relatively, driving the EEG to slower, higher amplitude patterns. Additionally, the anesthesia used in patients in this study represents clinically relevant administration protocols. The ability of these complexity measures to follow the consciousness transitions as well as standard spectral measures during routine anesthesia care highlights their utility.

Conclusions

Overall, we have demonstrated the utility of complexity measures to distinguish between the subtle EEG changes that occur with loss and recovery of response. The present study demonstrates their non-inferiority to standard spectral measures; thus, future work in hard-to-monitor patient populations is warranted to explore the full advantages of these measures. As the brain naturally transitions from wakefulness to sleep or has unconsciousness imposed with small molecules as in anesthesia, neural activity becomes more synchronous. Decreases in functional connectivity with anesthesia [64–67] impose regularity, periodicity, and consistency. Complexity measures of EEG signals are related to and capture these changes in functional connectivity well [6], especially in patients who exhibit behavior associated with unconsciousness. These observations in our work, consistent with others, highlight that complexity measures from nonlinear dynamics not only have utility as measures of anesthetic depth, but also support existing theories on neural correlates of consciousness transitions in the brain.