All experiment protocols involving humans were conducted in accordance with the principles of the Declaration of Helsinki. The current study was approved (No. ERB-C-1074-2) by the Institutional Review Board for Human Experiments at the Kyoto Prefectural University of Medicine (IRB of KPUM), and for this non-interventional and noninvasive retrospective observational study, informed patient consent was waived by the IRB of KPUM; patients were provided an opt-out option, of which they were notified in the preoperative anesthesia clinic. In our facility, the use of a BIS monitor is routine for adult and pediatric patients who undergo surgery involving GA. The anesthesiologists in charge of management did not receive notice of the study and planned the anesthesia methods for scheduled surgeries in accordance with our facility’s standard care protocol, without any feedback regarding the on-line analysis of processed EEG signals. Patients were not premedicated before anesthesia induction, in accordance with our facility’s standard protocol. An anesthetic gas monitor (IntelliVue G5 Anesthesia Gas Module, Philips, Amsterdam, Netherlands) was used for the measurement of end-tidal sevoflurane (et_SEV) and desflurane (et_DES) concentrations. The end-tidal anesthetic gas concentration (et_AG) was automatically recorded at 1-min intervals on a data server, then retrieved after the completion of anesthesia management. Based on the et_AG data for each minute, et_AG values at 20-s intervals were obtained by spline interpolation. In pediatric patients, GA was induced with an inhalational mixed gas, whereas in adult patients, rapid induction with propofol was employed (because of the irritating effects of desflurane on the respiratory tract during slow induction, all anesthesiologists in our facility refrain from the use of desflurane in children). After anesthesia induction, rocuronium (0.8–1.0 mg × kg^–1) was intravenously administered, and tracheal intubation was conducted. The anesthesia was maintained with approximately 2.5% sevoflurane or 6% desflurane, small doses of fentanyl (1 μg × kg^–1 per dose), and additional maintenance doses rocuronium (0.2 mg × kg^–1 at intervals of 20–30 min). The timing of sugammadex administration was significantly earlier in the adult_{_SEV} group than in the adult_{_DES} group (p < 0.05) (Table 1). The timing of the end of surgery was significantly earlier in the adult_{_SEV} group than in either adult_{_DES} (p < 0.05) or ped_{_SEV} groups (p < 0.01) (Table 1). In this EEG analysis, we focused on the 10-min period beginning around the end of surgery and ending with emergence from anesthesia; the study of EEG during the emergence process is important for preventing accidental intraoperative awakening. We analyzed the changes in correlations over time of various Poincaré plot parameters with et_AG, EEG parameters such as observed BIS, SEF₉₅, and total power (TP), and electromyography (EMG) parameter EMGlow.

A BIS Quatro sensor was mounted on the frontal region, in accordance with the manufacturer’s recommendations. The digitized EEG packets with a sampling frequency of 128 Hz were obtained through the serial output of the BIS monitor that sent a packet of sixteen sets of EEG μV data (32 bits) and eight packets per second (128 Hz); 8 s of EEG yielded 1,024 data points. Besides, a packet containing processed EEG valuables such as EMGlow (absolute power in the 70–110 Hz range, and, and values in decibel [dB] with respect to 0.0001 μV²) was obtained through the serial output of the BIS monitor once every second. For each 8-s data set, raw EEG signals were reconstructed from each digitized EEG packet and processed through six settings of FIR bandpass filters, f0: 0.5–47 Hz, f1: 0.5–8 Hz, f2: 8–13 Hz, f3: 13–20 Hz, f4: 20–30 Hz, and f5: 30–47 Hz (Figure 1A). For each FIR bandpass-filtered EEG dataset, the SEF₉₅, TP, and power spectrum were calculated. The Poincaré plot was constructed from the FIR-filtered EEG using 8-s epochs of the EEG signal (1,024 data points). The calculation of Poincaré plot parameters using FIR-filtered EEG A pair of EEG voltages with 1/128-s time lag (the shortest time lag under 128 Hz sampling rate of BIS monitor’s EEG) was plotted in the XY plane, first at a specific time (X: x_(k)) and then after a time delay (Y: x_(k+1)). To analyze the distribution of EEG patterns in the Poincaré plot, the standard deviation (SD) of the EEG voltage dispersion was measured along and perpendicular to the diagonal line of identity (SD1 and SD2, respectively, Figure 2). SD1 and SD2 represent the minor and the major semi-axes of this fitted ellipse. SD1 is the standard deviation of the distances of points from axis 1 and determines the width of the ellipse (short-term variability), while SD2 equals the standard deviations from axis 2 and length of the ellipse (long-term variability). The equations for the Poincaré plot parameters are as follows (Brennan et al., 2001; Golińska, 2013; Khandoker et al., 2013):

where SD(x_k) is a standard deviation of the time series x_k.

The SD1/SD2 ratio, which characterizes the sharpness of the scattered pattern, has been reportedly used for estimation of anesthesia depth (Hayashi et al., 2015a,b). However, in the present study, the Poincaré plot area (PP_A) was calculated with the equation SD1 × SD2 × π (Data Sheet 1, Supplementary Method, the example codes of Poincaré plot parameters–calculations with sample data using Python in Jupyter Notebook, and Processing). Our novel approach comprised dividing the PP_A of each frequency range by the PP_A of the f0 range (PP_{A_f0}) and defining the result as the PP_AR (PP_AR = PP_{A_fx}/PP_{A_f0}, fx = f1∼f5). The value of BIS was simultaneously collected from digital packets received from the BIS monitor, then recorded at 3-s intervals into the output data file.

Prior to this study, we created a software application named EEG Analyzer (downloadable for free use from our blog site, Science to Medicine, EEG Analyzer ver. 54_GP,¹) (Figure 1B), through which raw EEG signals are obtained from a VISTA A-3000 BIS monitor (Application revision 3.22, Medtronic, Minneapolis, MN, United States) in the legacy mode using an RS-232 interface to a personal computer (Surface Pro 4, Microsoft Co., Redmond, WA, United States). The Processing 3.5 software package (ver. 3.5.3,², MIT Media Lab, Massachusetts Institute of Technology, Cambridge, MA, United States) was used with the Apache Commons Mathematics Library (version 3.6.1,³, Apache Software Foundation, Forrest Hill, MD, United States) and a Java Virtual Machine (Oracle, Redwood Shores, CA, United States) with Java class libraries (javax.swing, java.AWT, and java.io packages). Through Processing’s integrated development environment, the program code for the EEG analysis was compiled using a Java Virtual Machine to build an execution file as a standalone application that functions in both Microsoft Windows 10 (Microsoft Co.) and Mac OS X (Apple, Inc. Cupertino, CA, United States).

For parametric regression with curve fitting function and creation of graphs, Microsoft Excel for Mac (ver. 16.16.5, Microsoft Corp., Redmond, WA, United States) and RINEARN Graph 3D (ver. 5.6,⁴, RINEARN, Kyoto, Japan) were used, respectively. Changes in various EEG parameters between two-time points, 10 min before emergence (EM_–10) and at the time of emergence (EM₀), were statistically compared using paired t-tests with InStat 3 (GraphPad Software, San Diego, CA, United States). Root mean square error (RMSE) and coefficients of determination (R²) were calculated as estimators for the regression analyses. For group comparisons of patients’ background characteristics, one-way analysis of variance (ANOVA) was first used to evaluate variation among group means. Only when the variation among group means was greater than expected by chance (p < 0.05), Student-Newman-Keuls multiple comparisons test was performed to evaluate the mean variation between the two groups. p < 0.05 was considered significant. For time-course comparisons, Student t-tests (paired, two tail p-value) were used. The p-values < 0.05 were considered significant.

Data were analyzed from 20 adults anesthetized with sevoflurane (adult_{_SEV} group), 20 adults anesthetized with desflurane (adult_{_DES} group), and 20 pediatric patients (aged 1–10 years) anesthetized with sevoflurane (ped_{_SEV} group). All 60 patients from the three groups were also included in a combined group (all_{_combined} group) (Table 1). All of the patients had the American Society of Anesthesiologists’ physical status values of 1–2 without any neurological diseases (Supplementary Table 1).

In the time-course from the hypnotic state to arousal during the 10-min process of emergence from GA (Figure 3), et_SEV in the adult_{_SEV} group, et_DES in the adult_{_DES} group, and et_SEV in the ped_{_SEV} group decreased from 0.97 ± 0.32, 3.71 ± 0.84, and 1.34 ± 0.75 to 0.17 ± 0.10, 0.68 ± 0.42, and 0.27 ± 0.16, respectively (p < 0.05 for all three groups). BIS values of the three groups increased from 52.9 ± 8.9, 41.9 ± 8.8, and 65.4 ± 14.9 to 80.5 ± 13.8, 84.7 ± 16.4, and 88.4 ± 9.1, respectively (p < 0.05 for all three groups). SEF₉₅ increased from 16.8 ± 2.5 Hz, 13.6 ± 1.9 Hz, and 20.8 ± 4.5 Hz to 21.7 ± 4.4 Hz, 20.3 ± 4.6 Hz, and 24.1 ± 3.3 Hz, respectively (p < 0.05 for all three groups). RelativeBetaRatio (RBR), a subparameter of BIS, tended to increase (p < 0.05). EMGlow, which corresponds to the EMG power of 70–110 Hz, also increased for all three groups (p < 0.05). In the all_{_combined} group, significant differences (p < 0.05) between EM_–10 and EM₀ were detected for et_AG, BIS, SEF₉₅, RBR, and EMGlow (Figure 3). Supplementary Tables 2, 3 show the detailed measurements in each of the four groups.

The time-course changes of TP were analyzed at each frequency band. In each frequency band, TP decreased over time, although small increases were observed in the last few minutes of TP_f5 in the adult_{_SEV} group, as well as in the last few minutes of TP_f4 and TP_f5 in the adult_{_DES} group (p < 0.05) (Supplementary Figure 1A). Next, Poincaré plot parameters (SD1 and SD2), which were obtained from the Poincaré plot-analysis of a pair of FIR-filtered EEG μV with a 1/128-s time lag, were analyzed at each frequency band. The time-course changes of SD1/SD2 depended strongly on the frequency band: gradual increases of SD1/SD2 were observed in f0, f3, f4, and f5, while gradual decreases were observed in f1 and f2 (Supplementary Figure 1B). However, the overall mean value changes of SD1/SD2 across all frequency bands were small compared to their standard deviations.

PP_A, which is calculated from SD2 × SD2 × π, changed dynamically during the emergence process (Supplementary Figure 1C). The absolute values of PP_A in the ped_{_SEV} group were approximately five times larger than those values in the adult groups. The changes in PP_A were not consistent over time: the values of PP_{A_f4} and PP_{A_f5} suddenly increased over the last 2–3 min in all three groups (p < 0.05, PP_{A_f5}). A similar increase in PP_{A_f0} was observed during the final phase of emergence. Next, the calculated absolute values of PP_A were significantly different among the three groups. Therefore, the adjustment was performed: the ratios of PP_A at each frequency range to PP_{A_f0} were calculated as PP_AR (PP_AR of each frequency-range/PP_AR of the f0: 0.5–47 Hz range) throughout the time course from the hypnotic condition to the awake state. The results showed that the PP_AR values for all three groups fit into the same plotting range. Drastic changes in PP_AR were detected among f1–f5 across all three groups (Figure 4). PP_{AR_f2} and PP_{AR_f3} significantly decreased (p < 0.05 for all three groups), while PP_{AR_f5} significantly increased (p < 0.01 for all three groups). The time-course changes of PP_{AR_f5} demonstrated notable changes 2–3 min before emergence: although the values of PP_{AR_f5} remained near 0 for the first 7–8 min, the values began to increase in the final few minutes before emergence. The mean changes in TP, SD1/SD2, PP_A, and PP_AR values in each frequency range over time are summarized in Figure 5. Changes of parameters in f5, such as PP_{A_f5} and PP_{AR_f5}, were remarkable among all five hierarchical frequency bands. In the all_{_combined} group, for et_AG, significant differences (p < 0.05) were detected in both SD1/SD2 and PP_AR at all five hierarchical frequency bands (Figure 4 and Supplementary Table 2B).

Next, using 4,000 data points obtained from 20 patients, the relationships of SD1/SD2 with et_AG and BIS were evaluated (Supplementary Figure 2A). Differences in the distributions of SD1/SD2, compared with et_AG and BIS, are present in each frequency range plot. Analysis of the relationship between SD1/SD2 and et_AG showed that when et_AG decreased, SD1/SD2 tended to be dispersed across all plots. Overall, both linear and logarithmic correlations between SD1/SD2 of six frequency bands and BIS in all three groups were low (R² ≤ 0.5). The relationships of PP_AR with et_AG and BIS were also evaluated (Supplementary Figure 2B). Scatter plots were created to assess PP_AR at five frequency bands, compared with et_AG. In plots between PP_{AR_f5} and et_AG, hyperbolic shapes with large variation at low et_AG values were observed in all three groups. Regarding the linear correlation between PP_{AR_f1} and BIS, the R² values in all three groups were less than 0.44. PP_{AR_f2} and PP_{AR_f3} showed biphasic changes against BIS in all three groups. Logarithmic correlations between PP_{AR_f4} and BIS showed R² values of 0.36, 0.59, and 0.27 for adult_{_SEV}, adult_{_DES}, and ped_{_SEV} groups. The relationships between PP_{AR_f5} and BIS in the three groups are summarized in Figure 6. Regarding the logarithmic correlations between PP_{AR_f5} and BIS, R² values were 0.67, 0.79, and 0.71, respectively [Figure 6A(3)]. Regarding the logarithmic correlation between PP_{AR_f5} and BIS in the all_{_combined} group, R² was 0.78 [Figure 6B(3)].

During the emergence process, the all_{_combined} group exhibited a gradual reduction of the average et_AG from 0.61 ± 0.30% to 0.12 ± 0.08%; furthermore, the average BIS gradually increased from 53.4 ± 14.7 to 84.6 ± 13.6 [Figure 3(1, 2)]. Because the relationship between PP_{AR_f5} and BIS fit well to a logarithmic curve with high R² values in all three groups, the relationship between log₁₀(PP_{AR_f5}) and BIS was fitted using linear regression (Figure 6C). Thus, the depth of anesthesia captured by the logarithmic value of PP_{AR_f5}, or by BIS, changed in a nearly linear manner from the state of anesthesia to emergence. Therefore, BIS is an index that can be easily understood by anesthesiologists because it converts the depth of anesthesia into a score ranging from 0 to 100; the logarithmic value of PP_{AR_f5} can also be easily converted to a 0–100 range. From the data of our 60 cases, PP_{AR_f5} = 0.3162 (logarithmic value −0.5) was defined as Poincaré plot-integrated score (PIS) = 100, and PP_{AR_f5} = 0.0031623 (logarithmic value −2.5) was defined as PIS = 50. The linear regression conversion formula for calculating the PIS value from PP_{AR_f5} then became PIS = 25 × log₁₀ (PP_{AR_f5})+112.5. Figure 4(6, 7) show the time-course changes of log₁₀(PP_{AR_f5}) and calculated PIS values in the four groups. R² between PIS and BIS in the three groups, which was between 0.67 and 0.79 (Figure 6C), is mathematically identical to the R² between the original PP_{AR_f5} and BIS [Figure 6A(3)]. In the all_{_combined} group, the regression equation BIS = 1.12 × PIS −11.68 was obtained (R² = 0.78) (Figure 6C).

While, during the emergence process, the all_{_combined} group exhibited a gradual increase of the average BIS from 53.4 ± 14.7 to 84.6 ± 13.6 [Figure 3(2)], the average EMGlow gradually increased from 2943.5 ± 465.7 to 5063.9 ± 1082.0 (Figure 3(5), Supplementary Table 2B). R² of the logarithmic regressions between EMGlow and BIS in the three groups were between 0.24 and 0.54 [Figure 7(1)], and R² of the exponential regressions between EMGlow and PP_{AR_f5} were between 0.39 and 0.64 [Figure 7(2)]. EMGlow is an absolute power at the range of 70–110 Hz. These results implied that, although the R² of the regressions between EMGlow and PP_{AR_f5} was not so great as that between BIS and PP_{AR_f5}, a part of EMG activity at the gamma range of 30–47 Hz probably influenced the calculations of BIS and PP_{AR_f5} with a non-negligible level.

In the case-by-case plots of BIS and PIS (Supplementary Figure 3), the overall time-course change of calculated PIS fit well to the change of measured BIS in each case, despite slight differences in correlations detected by RMSE and R² among the three groups (Supplementary Table 3). One representative case from each group of the three patient groups (A: adult_{_SEV}, B: adult_{_DES}, and C: ped_{_SEV}) is shown as a movie clip (Supplementary Video Clips) with spectrograms, time trends of BIS, PIS, and et_SEV or et_DES at 10 min before emergence from GA. Corresponding spectrograms of the above three cases from each group are shown with the corresponding time trends of BIS, PIS, et_SEV, or et_DES (Figure 8). Theta-to-alpha oscillations were observed before the dissipation points, followed by the emergence of beta-to-gamma oscillations with abrupt increases in BIS and PIS.