Mapping general anesthesia states based on electro-encephalogram transition phases

Cortical electro-encephalography (EEG) served as the clinical reference for monitoring unconsciousness during general anesthesia. The existing EEG-based monitors classified general anesthesia states as underdosed, adequate, or overdosed, lacking predictive power due to the absence of transition phases among these states. In response to this limitation, we undertook an analysis of the EEG signal during isoflurane-induced general anesthesia in mice. Adopting a data-driven approach, we applied signal processing techniques to track θ- and δ-band dynamics, along with iso-electric suppressions. Combining this approach with machine learning, we successfully developed an automated algorithm. The findings of our study revealed that the dampening of the δ-band occurred several minutes before the onset of significant iso-electric suppression episodes. Furthermore, a distinct γ-frequency oscillation was observed, persisting for several minutes during the recovery phase subsequent to isoflurane-induced overdose. As a result of our research, we generated a map summarizing multiple brain states and their transitions, offering a tool for predicting and preventing overdose during general anesthesia. The transition phases identified, along with the developed algorithm, have the potential to be generalized, enabling clinicians to prevent inadequate anesthesia and, consequently, tailor anesthetic regimens to individual patients.

1 Extended results

Band behavior during emergence
Although we identified a deterministic and ordered chain of EEG events during induction and maintenance of anesthesia, this chain did not extend to recovery.In particular, we found that the θ and δ activity did not have a characteristic and systematic evolution during emergence.This is illustrated in the average plots of the δ-and θ-relative powers in Figure S7A, where emergence starts at time 25 minutes.For the incremented and 1.5% protocols (left and center columns), during emergence, the relative delta power (green plot) increases to a value higher than the baseline value, while the relative theta power increases very softly (yellow plot).However, the constant 1% plot (right column) shows the opposite behavior: the relative delta power stays very low, while the relative theta power increases drastically.

1/f-parameters have characteristic behavior
We separated the EEG power spectral density into two components: the oscillatory component and the 1/f-component, with the form at ct+f p t where a t , p t , and c t are estimated on sliding time windows (Methods section 4.6, Fig. S7).The 1/f-component exponent p t had the same behavior during baseline and induction across protocols (Fig. S7B).It started at a relatively stable value between 1 and 1.4, before increasing rapidly after GA began.p t eventually stopped increasing, and immediately decreased.The maximum value reached by p t was very close across protocols: 2.07, 2.09, and 2.02 for the incremented, 1.5%, and 1% protocol respectively.However, the duration between start of GA and maximum p t varied across protocols: 12 minutes, 6 minutes 30 seconds, and 7 minutes for the incremented, 1.5%, and 1% protocol, respectively.The decrease slope of p t varied across protocols.a t and c t behaved globally the same as p t (Fig. S7C-D).
Moreover, the dynamics of p t and the relative δ-power P δ were weakly correlated (0.61 ± 0.22 Pearson correlation, see Supplementary section 1.2).We found that the 1/f exponent p t was weakly linked to the δ relative power P δ (0.61 ± 0.22 Pearson correlation, 0.58 ± 0.21 Spearman correlation, see Fig. S7F).On average, the correlation between p t and P δ was higher for the incremented protocol and the 1.5% protocol than for the 1% protocol.

Correlation of 1/-f and oscillatory EEG components
To study the independence of the two components of the spectral decomposition, we computed the correlation of the energy estimated using the area under the 1/f and the oscillatory components.For the latter, we summed the area under the θ and δ Gaussians.We found an insignificant correlation (0.10 ± 0.2 for the incremented Inj Iso protocol, 0.11 ± 0.17 for the Inj Iso = 1.5% protocol, and -0.24 ± 0.2 for the Inj Iso = 1% protocol), with an overall -0.002 ± 0.25 (Fig. S5F).The present analysis suggests that the two components should be studied separately and bring complementary information.

Correlation between IES and γ-rebound parameters
In Results section 2.4, we uncovered a phenomenon happening in the EEG γ frequency range.We found that the area under the γ−power was correlated to the total time spent in IES (∆ IES ), which we called cumulative IES duration.We now detail additional γ−rebound parameters and the possible correlation with IES duration occurring during GA.We define the γ-rebound duration (∆ γ ) as the time difference between the end of the last γ spindle and the beginning of the first γ spindle during recovery.When we computed ∆ γ in recordings, we did not observe a statistical difference between the anesthesia protocols (Fig. S8A).We next explored whether the γ-rebound durations could be correlated with ∆ IES .To address this question, we plotted the points (∆ IES , ∆ γ ) for each recording containing a γrebound.A linear regression y = ax + b (Fig. S8B) led to a = 0.58 and b = 3.69 with a determination coefficient R 2 = 0.2, revealing no obvious correlation between the γ-rebound duration and IES duration.We note that in 8/23 recordings, the γ-rebound extended beyond the end of the recording, preventing an accurate estimation of ∆ γ .To further study the possible correlation between the γ-rebound and ∆ IES , we now introduce how the γ-power p γ (k, t 1 , t 2 ) is computed for recording k on the time interval [t 1 , t 2 ]: where S k,γ is the EEG signal of recording k bandpass filtered in 50−70Hz, which is the frequency domain activated during γ rebounds.We recall the discretized version: where f S is the sampling frequency.We now introduce the recovery power p recov γ during the recovery phase for the gamma band as the power computed for S γ (t) between the end (time τ stop Iso (k)) of anesthesia and the end (time τ stop rec (k))) of recording k as To further investigate the possible correlation between the gamma recovery power p recov γ (k) and ∆ IES (k), (Fig. S8C), we performed a linear regression y = ax + b and found a = 32.4,b = 84.5 and R 2 = 0.5, confirming a weak linear dependency.Finally, we define the normalized γ-recovery power between the end of anesthesia and the end of recording as where τ 0 (k) is the beginning time of the EEG recording k.We performed a linear regression y = ax + b (Fig. S8D) leading to a = 0.17, b = 0.81, and R 2 = 0.2, suggesting no correlation.
To conclude, we found no statistical correlation between the IES duration and γ−rebound durations, or IES duration and the γ−rebound energy, normalized or not.

EEG γ-rebound is not EMG contamination
During the γ-rebound described in Results section 2.4, EEG γ bursts were synchronized with EMG bursts (Fig. S10A-B).But spectral decomposition of these signals showed key differences (Fig. S10C-D).The EEG power spectrum was much more powerful, with power going up to 200 µV as opposed to 20 µV for the EMG.Moreover, during EEG bursts, the γ activity was in a narrow frequency domain, around 60 Hz.On the contrary, in EMG bursts, the entire 50−100 Hz range was active.Therefore, EMG contamination was not the cause of the high-power activity of the γ-rebound EEG.However, because of their synchronicity, we hypothesize that EEG and EMG bursts are linked, possibly caused by a common phenomenon that has yet to be unveiled.

δ disappearance predicts the emergence of suppressions several minutes in advance
We investigated the dynamics of the δ-and θ-rhythms and their roles in predicting BS appearance.Interestingly, we found that the δ-rhythm disappearance could anticipate the appearance of BS (Fig. S9).
To obtain a predictive estimate of the arrival of BS, we computed the normalized power P δ of a sub-delta band in the range 2−4 Hz (see eq. 7).We found that when P δ leaves for the first time a region where its value exceeds a threshold T δ = 0.15, BS appear in average with a delay of few minutes (Fig. S9A), in 81% of the cases.Interestingly, this delay was, on average, 4.5 minutes for the 1.5% protocol, and 1.9 min for the incremented protocol (Fig. S9B).
1.7 Repeated isoflurane exposure does not induce group differences in response to anesthesia We noticed that one individual mouse could respond differently to anesthesia when anesthetized several times.We however found no evidence that this variability was a consequence of repeated exposure to isoflurane.For instance, two mice received the constant 1% protocol three times in a row.Interestingly, the intra-individual response varied greatly.While the first mouse had a cumulative IES duration of 0 min, 0 min and then 1.32 min (suggesting that exposure to isoflurane increases isoflurane sensitivity), the second mouse had a cumulative IES duration of 2.4 min, 0.85 min, and then 1.5 min (suggesting that isoflurane exposure does not have a strong impact on isoflurane sensitivity).We therefore found no general trend as to whether isoflurane exposure history increases or decreases isoflurane sensitivity over time.We further evaluated this question by comparing the cumulative IES duration of two groups of recordings undergoing the incremented protocol.Indeed, the cumulative IES duration can be used as a quantification of sensitivity to anesthesia.In the first group, the mice were naive, they had not been anesthetized before (except during surgical implantation, weeks prior to the experiment) (n=4, mean cumulative IES duration of 4.7 ± 2.3 min).In the second group, the mice had recently been anesthetized using the constant 1.5% protocol (n=4, mean cumulative IES duration of 2.8 ± 2.1 min).The Mann-Whitney U rank two-sided test returned a p-value of 0.47, meaning we cannot reject its null hypothesis, which is that the two groups have the same underlying distribution.We conclude that in our recordings, previous isoflurane exposure does not predict the response to isoflurane.
2 Extended methods

Choice of the threshold value to segment IES
IES detection relies on the choice of a threshold value T IES .Choosing this threshold depends on the amplitude of the EEG signal, which can vary with the type and position of the electrodes, and with inter-individual variability [1].We describe here how we selected this threshold and adapted the value for each recording.We express the threshold for IES as where RMS EEG is the Root Mean Square (RMS) of the EEG signal computed over the entire duration of each recording and the parameter r IES has to be determined.If r IES is too small, then IES epochs are not correctly detected (Fig. S1A).If r IES is a relevant value, then IES epochs are correctly detected (Fig. S1B).And if r IES is too big, then non-IES epochs are wrongly detected as IES (Fig. S1C).To determine r IES , we first varied r IES between 0 and 1.5 , and for each r IES value, we ran the IES detection algorithm (Methods) during the entire duration of each recording.For each r IES and each recording k, we computed the total time ∆ k
The division in G 1 and G 2 was done by visually inspecting the data.The average curves for each subgroup are defined by where n 1 (resp.n 2 ) is the number of recordings in G 1 (resp.G 2 ), and the relative threshold values r IES were between 0 and 1.5, with an increment of 0.05.The result is plotted with the standard deviation in Fig. S2B1-B2.For mice in G 1 , the average curve ∆ 1 (r IES ) follows a fast increase, followed by a region with a much smaller slope value that we refer to as plateau P l G1 (red in Fig. S2B1).The curve ends with another fast increase.In contrast, the average curve ∆ 2 (r IES ) stays equal to zero, before rapidly increasing (Fig. S2B2).
We then approximated the first derivative of the average curve ∆ 1 : with r ′ n = rn+r n+1 2 .The plateau P l G1 is visible as a well in the first derivative of ∆ 1 (Fig. S2C) and provides an admissible range for the relative threshold r IES , i.e. a set of values with which IES detection is satisfactory.To minimize false positives (IES wrongly detected), we chose the relative threshold at the left bound of the plateau, r * IES = 0.7.We note that with this value, no IES are wrongfully detected in group G2 (Fig. S2B2).

Influence of the correcting factor c t on the 1/f-decomposition
In Methods sections 4.5.1 and 4.5.2,we decomposed the power spectrum into 1/f and oscillatory components on sliding time windows.In this section, we discuss the influence of the parameter c t on the 1/f -component estimation.For that purpose, we compared the fits y 1 (f ) = a c+f p and y 2 (f ) = a f p in their capacity to approximate the 1/f -component in the frequency domain.The parameter estimation for a, c and p for y 1 and a, p for y 2 is described in the Methods section 4.5.1.These fits are performed for all recordings in our dataset and computed over successive time windows with a width of 60 seconds and an overlap of 30 seconds (Fig. S3).We now describe how we evaluate the local error made by each fitted function on a sliding time window: The 1/f-component of the EEG signal ỹk t,w is estimated from the time window W w (t) of recording k by parameterizing y 1 (resp.y 2 ) ỹk t,w , resulting in y k 1,t,w and y k 2,t,w respectively.The local normalized error made by y 1 on time window W w (t) of recording k is defined by and the local normalized error made by y 2 on W w (t) of recording k is: The discretized version of the local normalized error is given by: Err(ỹ, y j , t, w, k) =     3.6 ± 5.3 3.8 ± 3.1 8.2 × 10 −6 τ S − τ disp δ 1.9 ± 1.8 4.5 ± 3.1 6.0 ± 0.3 3.2 ± 2.7 9.9 × 10 −5 Table 1: Time difference between key EEG events (mean ± std in minutes), and p-value of the one-sided Wilcoxon signed-rank test.
Figure S2: Choice of relative threshold r IES for IES detections.Detected IE duration ∆k IES versus relative threshold r IES for recordings with IES (A1), and recordings without IES (A2).(B1) Average curve (full) with the standard deviation (blue area) of the curves presented in (A1).(B2) average with the standard deviation of curves in (A2).(C) Derivative of the mean detected IE duration ∆1IES .The admissible range for the relative threshold r IES is the slowest increase region in (B1), associated with a well (C).We chose the value r * IES = 0.7 (red), located at the left of the well.

Figure S3: 1
Figure S3: 1/f-component in the frequency space.(A) The power spectral density (dark blue) is fitted by an estimate (y) computed over one minute (IRASA 1/f component extraction in light blue), y 1 = a c+f p (red) and y 2 = a f p .(y) (green).The parameters a, c, and p are estimated in the Methods section.(B) linear scale.

Figure S4 :
Figure S4: Local errors associated with y 1 = a c+f p and y 2 = a f p decomposition, computed over the EEG recording.(A) log and (B) linear scale.

Figure S5 :
Figure S5: Estimation of oscillatory and fractal parameters on a time window.(A) Input window of the signal.(B) Power Spectral Density of the input signal.(C) Estimation of the fractal component and the fractal parameters a t , p t , c t .(D) Estimation of the oscillatory component.(E) Decomposition of the oscillatory component in Gaussians and parameter extraction of the θ-and δ-rhythms.(F) Correlation between the areas under the fractal and oscillatory components.

Figure S6 :
Figure S6: EMG and associated spectrogram over an entire recording.Loss of Movement (LOM), indicated by the first vertical black arrow, and Return of Movement (ROM), indicated by the second black arrow.

Figure S7 :
Figure S7: Trajectories of spectral parameters for anesthesia protocol: relative powers and 1/f-component at ct+f p t .(A) Relative δ-power (dark green) and relative θ-power (yellow).(B) exponent p t , (C) intercept a t and (D) correction term c t .(E) Distribution of Pearson and Spearman correlations between p t and P δ .Error bands indicate the 95% confidence intervals computed using the t-distribution.13

Figure S8 :A 1 Figure
Figure S8: Statistical parameters associated with the γ−rebound.(A) The γ-rebound duration is not statistically different among the three anesthesia protocols (two-sided Wilcoxonrank U test).(B) γ-rebound duration is not linearly correlated to IES duration.Recovery γ power (C) and normalized recovery γ-power (D) are not linearly correlated to IES duration.

Figure S10 :
Figure S10: Magnification of EEG and EMG during a γ-rebound.(A) EEG signal.(B) EMG signal.(C) EEG spectrogram.(D) EMG spectrogram.The EMG spectrogram differs greatly from the EEG spectrogram, showing that the EEG γ-rebound is not EMG contamination.