Using Permutation Entropy to Measure the Electroencephalographic Effects of Sevoflurane

Test subjects included 19 patients (aged 18–63 yr) with American Society of Anesthesiologists physical status I or II scheduled to undergo elective gynecologic, general, or orthopedic surgery.19Patient exclusion criteria were preoperative use of medication acting on the central nervous system; excessive weight or a history of gastroesophageal reflux that would not permit gaseous induction with sevoflurane; a history of cardiac, pulmonary, hepatic, or renal disease; and use of any premedication. All subjects fasted for at least 6 h before anesthesia and received no premedication. All subjects were recruited from Waikato Hospital, Hamilton, New Zealand. Informed consent was obtained from all subjects following Waikato Hospital ethics committee approval.

Before electrode application, the skin was carefully cleaned with an alcohol swab and allowed to dry; the electrode skin impedance was checked to be less than 7.5 kΩ. A composite electrode composed of a self-adhering flexible band holding three electrodes was used to record the electroencephalography between the forehead and temple. The spectral entropy was measured with a plug-in M-Entropy S/5 Module (Datex-Ohmeda). Response entropy and state entropy were sampled at 0.2 Hz. The sampling rate of electroencephalographic data was 100/s. Inspired and expired sevoflurane concentrations were measured at the mouth and sampled at 100/s. The data were recorded on a laptop computer and stored for later off-line analysis using MatLab (version 7; MathWorks, Natick, MA) computational and data analysis software.

All of the patients were connected to a closed anesthesia breathing circuit via  a facemask, and fresh gas flow was set at 4 l/min. The patients were preoxygenated to the satisfaction of the anesthetist in charge. Then, sevoflurane was delivered by vaporizer at 3% inspired for 2 min, followed immediately by 7% inspired concentration. State entropy and response entropy were recorded for comparison with awake values. The time at which response entropy decreased to 20 or less was noted, and 7% sevoflurane was continued for a further 2 min. The sevoflurane was then turned off until the response entropy returned to a value of 70 (lightening). All data in the current study are based on this single deepening and lightening protocol. No attempt was made to rouse the subjects, and no supplemental medications were administered. Loss of consciousness (LOC) was assessed by loss of response to verbal command. This data set has been used previously to model the pharmacokinetic–pharmacodynamic effects of sevoflurane by using the spectral entropy of the electroencephalographic signals.19 

Several methods based on frequency band, modeling, and time-scale transform have been proposed to remove or reduce the artifacts in scalp electroencephalographic recordings. In this study, a filter based on wavelet transform was designed to preprocess artifacts embedded in the electroencephalographic recordings.20The wavelet-based de-noising method is composed of three steps21: (1) the wavelet transform of a signal x  (t ), (2) thresholding the wavelet coefficients, and (3) the inverse wavelet transform to obtain the de-noised signal. Each step is described briefly below.

1. In this study, we used a six-level discrete wavelet transform, using a Sym8 wavelet, which was applied to electroencephalographic data epochs 10 s in length. Previous work has shown that the Sym8 better approximates the electroencephalographic waveform, resulting in preservation of the electroencephalographic waveform after application of the filter.22 

2. Optimum threshold selection is core to the success of the wavelet-based de-noising method. Several algorithms have been proposed to estimate the threshold for wavelet coefficient-based noise removal, such as data-adaptive wavelet thresholding, block thresholding, and bayesian methods.21In this study, we used bayesian block thresholding because previous work has shown that this approach outperforms the classic data-adaptive method.23In the bayesian approach, a distribution is calculated based on the wavelet coefficients, which captures the sparseness of wavelet expansions. A suitable bayesian rule is then applied to estimate posterior distribution of the wavelet coefficients.23The criteria of the thresholds depend on the bayesian formalism and assumed prior distribution, the details of which have been described previously.22,24–26 

3. After filtering using the bayesian estimation method described in step 2, the retained information, above thresholds, was used to reconstruct a new electroencephalographic series via  an inverse discrete wavelet transform. The electroencephalographic data in the frequency bands of approximately 25–50, 12–25, 6–12, 3–6, 1.5–3, and 0.5–1.5 Hz were obtained. The slow frequency band of 0–3 Hz was removed in this study because artifacts derived from electrical activity of the heart, muscles, and eyes often lie at this frequency band (BIS VISTA27 ^(p458)). The electroencephalographic data in the frequency band over 25 Hz often includes a large amount of electromyographic artifact28and was also discarded before further analysis. Therefore, the electroencephalographic data within the 3- to 25-Hz band was used to calculate the AE and PE values in the following sections.

Considering the chaotic and nonstationary nature of electroencephalographic data, AE29has been applied instead of spectral entropy30to measure the predictability of electroencephalographic series. AE is based on the similarity of events in phase space and is a favorable method to address randomness of a dynamical system.14The predictability of subsequent amplitude values of the electroencephalographic signals can be quantified based on the knowledge of the previous amplitude values. If the electroencephalographic series is regular, the position of a particular point can be predicted by using its previous points, whereas in an irregular electroencephalographic series, the position of the point cannot be easily predicted. It was demonstrated that AE of electroencephalographic recordings could be used to quantify the effect of anesthesia drugs on the brain activity better than other descriptors.14The number of previous (lagged) points used in making the prediction is the embedding dimension (m ). The AE looks at sequences of length m  and then establishes the negative logarithmic probability that these sequences predict a new sequence of m + 1 points to within an error range of r , typically set at 0.2 × SD. For a regular signal, most sequences can successfully predict the next data points; therefore, the AE is low. In turn, for an irregular signal, there are few successful predictions, so the AE is correspondingly high. The exact value of the AE (m , r , N ) depends on the values chosen for the three parameters of the statistic: N  (number of samples), m  (embedding dimension), and r  (noise threshold).31In this study, N = 1,000 points, r = 0.2 × SD, and m = 2 were selected in the light of previous work.14 

Ordinal time series analysis is a new method to describe the characteristics embedded in a complex time series. This method transforms a given time series into a series of ordinal patterns each describing the order relations between the present and a fixed number of equidistant past values at a given time.32Given a scalar time series x^t(t = 1,2, …), an embedding procedure forms vectors X^t[x¹,x^t+τ, …, x^t+mτ ] with the embedding dimension, m , and the lag, τ. Then, X^t  can be arranged in an increasing order. For m  different numbers, there will be m ! possible order patterns, which are also called permutations.15To quantify and visualize the changes in time series, Bandt and Pompe15have proposed PE based on the Shannon entropy. Considering each permutation as a symbol, the vectors X^t  can be represented by a symbol sequence; the distinct number of symbols (J ) should be less than or equal to m !, namely J ≤m ! For the time series x^t , the probability distributions of the distinct symbols are defined as p¹…, p^j ; the PE of this time series is defined by

The corresponding normalized entropy can be defined as follows:

The largest value of H^p  is one, meaning all of permutations have equal probability. The smallest value of H^p  is zero, meaning the time series is very regular. That is, the smaller the value of H^p  is, the more regular the time series is. PE refers to the local order structure of the time series, which can give a quantitative complexity measure for a dynamical time series. Mathematical details of the PE can be found in the references.15,32–34 

Permutation entropy calculation depends on the selection of time interval N  and embedding dimension m . Similarly to AE, interval N  was selected as 1,000 points (data length of 10 s in this study). When m  is too small (less than 2), the scheme will not work because there are too few distinct states. Often, for a long electroencephalographic recording, a large value of m  is better. However, this study concentrates on the detection of dynamical changes in the electroencephalographic recording, so too large a value of m  (greater than 10) is inappropriate.16In this study, we found that m = 6 was appropriate. The parameter selection of PE is addressed in appendix 1.

The correlation between anesthetic drug effect index and anesthetic drug concentration provides construct validity for anesthetic drug effect monitoring.19,35In this study, we used a standard pharmacokinetic–pharmacodynamic model to describe the relation between sevoflurane concentrations and the electroencephalographic response (measured by the PE and AE). This was done by modeling the movement of sevoflurane from the arterial blood (end-tidal) using the first-order rate constant, K^e0. Briefly, the effect site partial pressure is estimated by a first-order effect site model19:

where C^etis the end-tidal concentration of the drug, C^effis the sevoflurane concentration at the effect site, and K^e0is the first-order rate constant for efflux from the effect compartment. The C^effis estimated by iteratively running this above model with a series of K^e0steps. For each iteration, a nonlinear inhibitory sigmoid E^maxcurve is fitted to the data by the following equation19:

where Effect is the processed electroencephalographic measure, the E^maxand E^minare the maximum and minimum Effect for each individual patient, and EC⁵⁰is the sevoflurane concentration at which Effect is midway between this maximum and minimum. γ describes the slope of the concentration–response relation. K^e0is determined from the iteration yielding the greatest coefficient of determination (R  ²) for measured and modeled electroencephalographic Effect for each patient.35A nonlinear inhibitory sigmoid E^maxcurve was fitted to the brain (effect site) concentration–PE or –AE relation. From the fitted curve, values of pharmacodynamic parameters describing this relation were derived, including γ and EC⁵⁰. The coefficient R  ²is calculated by

where yⁱ  and ŷⁱare PE or AE for a given time and the corresponding model prediction, respectively. The ȳ is the average of all the measurements (PE and AE).

The raw electroencephalographic signals were filtered between 3 and 25 Hz and divided into epochs of 10 s duration. For each patient, the PE and AE series were computed and effect site concentrations (C^eff) were estimated with pharmacokinetic–pharmacodynamic modeling. The correlations between effect site concentrations and PE and AE were investigated with the model-independent prediction probability (P  ^K) that was described by Smith et al.  36The P  ^Kvalue describes the probability that the PE or AE predicts effect site concentration. Considering the decrease of PE and AE with the increase of sevoflurane concentration, the prediction probability (P  ^K) can be calculated by9 

where P  ^c, P  ^d, and P  ^txare the respective probabilities that PE or AE values and effect site concentrations are a concordance, a discordance, or an x-only tie. In this study, we randomly selected a pair of PE values from before and after LOC and then determined whether this pair of PE values correctly predicted the trend of effect site concentration (up or down). By repeating the above procedure 300 times, the P  ^c, P  ^d, and P  ^txvalues were estimated, so a P  ^Kvalue for each patient could be obtained. The same procedure was used to deal with AE. The details of this procedure can be found in reference 9. A P  ^Kof 1 for the electroencephalographic parameter (PE or AE) means that the indictors can correctly predict the estimated sevoflurane effect site concentration. Alternatively, a P  ^Kvalue of 0.5 means that the electroencephalographic parameters (PE or AE) are no better than chance at predicting the estimated sevoflurane effect site concentration.

The Kolmogorov-Smirnov test was used to determine whether data sets were normally distributed. PE and AE values during the awake state (averaged over a 10-s epoch, 50 s before LOC for all 19 subjects) were compared with values during anesthesia (averaged over a 10-s epoch, 50 s after LOC for all 19 subjects). The PE and AE values at the awake and anesthetized states were analyzed using the nonparametric Wilcoxon test. To determine the changes in PE and AE as anesthesia was deepened, mean PE and AE mean values over 10-s epochs were compared at four time points: awake (50 s before sevoflurane delivery), anesthetic LOC (sevoflurane delivery to LOC), LOC, and LOC + 50 s. The correlation (R  ²) between the entropy values and the sevoflurane effect site concentration after pharmacokinetic–pharmacodynamic modeling was calculated for each patient, and the average of the correlation was calculated for all 19 patients. This analysis was performed to compare the ability of each algorithm track sevoflurane effect site (brain) concentration. In addition, a relative coefficient of variation (CV) (the ratio of the SD to the mean) was used to determine the predictive performance of the derived entropy estimators PE and AE for sevoflurane effect compartment concentrations.

Figure 1shows the time course of measurement for one patient. The electroencephalographic series during deep anesthesia is more regular than at low anesthetic concentrations (see the enlarged sections of electroencephalographic data in fig. 1A). Figure 1Bshows that PE and AE of the filtered electroencephalographic epochs decrease with increasing sevoflurane concentration and then reach a plateau after approximately 3 min. It is clear that both PE and AE track the gross changes in electroencephalographic data with increasing anesthetic drug effect (fig. 1B) and end-tidal sevoflurane concentration (fig. 1C).

To evaluate the performance of wavelet filtering to remove artifacts in the electroencephalographic data, the baseline variabilities of PE and AE before anesthetic delivery were estimated by calculating the CV. The CV values of PE and AE before wavelet filtering were 0.007 ± 0.008 and 0.192 ± 0.079 (n = 19), respectively. After wavelet filtering, the CV for PE was 0.004 ± 0.006, compared with 0.059 ± 0.047 for AE (n = 19). Clearly, the difference in PE CV values before and after filtering is smaller than the corresponding difference in AE CV values. Assuming the difference of CV values before and after wavelet filtering is derived from the artifacts in the electroencephalographic recordings, these results demonstrated the robustness of the PE algorithm to noise in the electroencephalographic data.

The effect of increasing sevoflurane concentration on AE and PE was quantified by comparing awake and anesthetized values across all subjects. The median (95% confidence interval) PE values during the awake and anesthetic epochs were 0.90 (0.88–0.91) and 0.66 (0.63–0.73), respectively; the median (95% confidence interval) AE values were 1.24 (1.04–1.33) and 0.80 (0.66–0.87), respectively. The differences in parameter values between awake and anesthetized states were statistically significant for both PE (P < 0.0001) and AE (P < 0.005, Wilcoxon test). The very small 95% confidence limits for PE calculated during the awake electroencephalographic epochs further highlights the resistance of the PE algorithm to electroencephalographic artifact. To further examine the ability of PE and AE to differentiate awake from anesthetized states, the prediction probability P  ^Kwas calculated. The individual P  ^Kvalue for each patient was calculated for the PE and the AE, and the means of P  ^K(95% confidence intervals) for PE and AE were 0.86 (0.82–0.89) and 0.79 (0.74–0.83). The difference between parameters was statistically significant (P < 0.05).

Plotting PE versus  AE (1,425 pairs from all patients) revealed an approximate linear relation between PE and AE, as shown in figure 2, demonstrating the ability of both PE and AE to track gross changes in anesthetic effect. Figure 3shows the relation between the point of LOC and PE/AE. LOC typically occurred before the rapid decrease in PE and AE (0.830 ± 0.069 [mean ± SD] at LOC for PE, compared with 1.085 ± 0.23 at LOC for AE). The variability in PE at LOC was smaller than for AE, as shown in the above SD values. This suggests that PE may be a more robust indicator of LOC. To compare the changes in PE and AE as anesthesia was deepened, PE and AE values at the awake (10-s epochs), anesthetic LOC, LOC, and anesthetic (after LOC + 50 s, 10-s epochs) states were analyzed for each patient, and a box plot was constructed (fig. 4). The PE values for all patients averaged 0.935 ± 0.003, 0.890 ± 0.041, 0.830 ± 0.069,and 0.702 ± 0.113 (mean ± SD) in the four states, respectively; the corresponding AE values were 1.502 ± 0.062, 1.175 ± 0.291, 1.085 ± 0.230, and 0.785 ± 0.218. The PE and AE values decrease monotonically from the awake to the anesthetized state; however, the root mean squared error of PE linear fitting is 0.029, compared with 0.165 for AE.

The sigmoid dose–response relations between the PE and AE values and the end-tidal sevoflurane concentrations for all subjects are plotted in figure 5. Note that these data do not include any electroencephalographic burst suppression states. The main difference between these plots is the obvious accumulation of values in the region of low AE and low sevoflurane concentration in the AE plot (fig. 5B). This reflects the relative susceptibility of AE to artifacts in the awake state (e.g. , blinks and forehead muscle activity) and the resistance of PE to these same artifacts.

Figure 6shows the data fit for the sevoflurane (effect site) concentration and the PE and AE values for the pharmacokinetic–pharmacodynamic modeling for one patient with an exponential function. Using the pharmacokinetic–pharmacodynamic model, the effect site sevoflurane concentration can be calculated, so the relation between PE and AE values and the effect site sevoflurane concentration can be obtained. Inspection of the individual data fits showed that the inhibitory sigmoid E^maxmodel approximately described the relation between the sevoflurane (effect site) concentration and the PE and AE values. The modeled parameters are listed in table 1, including E^max, E^min, γ, C^et, EC⁵⁰, K^e0, and correlation coefficient R  ². The parameter values are similar for PE and AE, although the correlation (R  ²) between the entropy values and the sevoflurane effect site concentration was 0.89 ± 0.07 (n = 19) for PE compared with 0.60 ± 0.14 (n = 19) for AE. It is worth noting that for the case shown in figure 6, above a sevoflurane concentration of approximately 1.7%, the PE decreases more slowly than indicated by the fitted curve.

In this study, the effect of sevoflurane on brain activity was demonstrated by the PE and AE values of electroencephalographic recordings. The PE and AE values decrease with increasing sevoflurane concentration, and either algorithm could be used to detect gross changes in sevoflurane effect. The variability at baseline, the prediction probability, sensitivity at LOC, the correlation of both parameters and its related measures of anesthetic drug effect on the electroencephalogram (pharmacokinetic–pharmacodynamic) were tested; the statistical analyzes demonstrate several advantages of PE over AE.

First, PE is less sensitive to “artifacts” in the electroencephalographic signal during the awake state. This was reflected in the lower CV for PE and in the higher correlation values for the modeled relation between sevoflurane effect site concentration and PE. Second, PE provides a more robust estimate of the point of LOC. Although both PE and AE could distinguish between the awake and anesthesia states (as shown in fig. 4), the variability of PE at the moment of LOC was smaller than AE, and PE monotonically decreased from awake to LOC. The greater predictive value of PE for differentiating subject state was also reflected in the slightly higher P  ^Kvalue for PE compared with AE. Third, the PE is more correlated to the sevoflurane effect site concentration. The correlation (R  ²) in table 1shows the PE can accurately reflect the effect of the sevoflurane effect site concentration on brain activity. Finally, PE calculation is simple, fast, and robust. PE calculation simply considers the order relation between the values of a time series instead of the values themselves; the value of PE is based on the distribution of ordinal patterns. Using the method shown in figure 1Bto process an electroencephalographic trace of 10 min requires less than 258 s for PE on a P4 1.6-GHz personal computer including wavelet filtering, compared with 372 s for AE. This case analysis shows that the PE method is more readily amenable for use in a real-time electroencephalographic monitor.

Approximate entropy is based on the similarity of events in phase space and is a favorable method to address randomness of a dynamical system. Because AE specifies a noise threshold, it may be better than spectral entropy in the quantification of complexity of electroencephalographic recordings.14,37The disadvantage of AE is that it is heavily dependent on the record length and is often lower than expected for short records. Another disadvantage is that AE lacks relative consistency.31PE is a new statistical parameter that also quantifies the amount of regularity in electroencephalographic data. The first step of PE calculation is to transform an electroencephalographic series into a series of ordinal patterns, which implies that a nonstationary series is transformed to a nearly stationary ordinal series. As a result, PE is less affected by the amplitude of the electroencephalographic data. Furthermore, the PE should be less sensitive to noise embedded in electroencephalographic recordings, which was confirmed by the findings of this study (appendix 2). On the other hand, the correlation coefficient of 0.8 between the PE and AE shows that the PE can be applied to quantify the complexity of electroencephalographic signals instead of AE. These factors suggest multiple possible applications of this new tool for the analysis of electroencephalographic data,33,38as suggested by case studies successfully separating consciousness from unconsciousness during anesthesia17and indicating increases and decreases in anesthetic level.18 

At surgical levels of anesthesia, some high-frequency artifacts such as electromyographic may interfere with the electroencephalographic PE values. The current data cannot reveal the details on this issue, which should be discussed in future work. From a practical viewpoint, this is a difficult issue because it relates to the identification of artifacts embedded in electroencephalographic series. The PE description of the burst suppression electroencephalographic pattern is another possible problem. With increasing anesthesia concentration, the electroencephalographic waveform changes into a burst suppression pattern, resulting in an increase in PE (fig. 7). Clearly, the AE of burst suppression better reflects very deep anesthesia than PE in this case, although further work is required to draw firm conclusions. The PE method may require a burst suppression ratio component, similar to that used in the BIS monitoring system. The phase relations between the component waves of different frequencies that make up the composite electroencephalogram are considered in the BIS. In fact, the phase information is related to the ordinal pattern in the electroencephalographic data, so the PE is similar to the BIS in this respect. At low concentrations of sevoflurane, such as 1–2%, the PE approximates to the plateau (see the end of figs. 1B and C), which is similar to the BIS.39Therefore, PE does not accurately track all levels of anesthesia.

In conclusion, when comparing the performance of PE and AE as measures of anesthetic drug effect on the electroencephalogram, it was found that PE performed better than AE during sevoflurane anesthesia on a number of levels: the baseline variability of PE was lower than that of AE; the prediction probability P^Kwas higher than AE; and pharmacokinetic–pharmacodynamic modeling of PE and AE versus  the sevoflurane (effect site) concentrations revealed that the correlation (R  ²) between the PE values and the sevoflurane effect site concentrations is higher than AE. These differences largely reflect the resistance of PE to artifact in the electroencephalographic recording. Furthermore, the fast computation of PE ensures that this new measure can be applied to clinical real-time on-line monitoring of anesthetic drug effect. It is suggested that PE may be an excellent candidate for designing a monitor system for depth of anesthesia when integrated with other disparate descriptors of the electroencephalogram.

The parameter selections of PE should be discussed. The calculation of PE depends on the length of the epoch (N ) and embedding dimension (m ), m ! < N . In Bandt and Pompe,15the authors recommend m = 3, …,7. It was found that m = 3 and 4 may still be too small and a value of m = 5, 6, or 7 seems to be the most suitable for an electroencephalographic series.16 

In reference 39, prediction probability (P  ^K) between the AE and effect site concentration was used to assess the parameter selection for AE calculation. A similar approach was used to the parameter selection for PE calculation. Figure 8shows how varying the embedded dimension m  and length of epoch N  affect the PE calculation at the awake state. Given the embedded dimension m = 6, the prediction probability (P  ^K) is the biggest when the length of the epoch is 1,000 (10 s) (fig. 8A), because 7! > 1,000, the Somers d statistic is set as zero, so P  ^K= 0.5 for m = 7, 8, and 9. Figure 8Bshows that the prediction probability (P  ^K) at the embedded dimension m = 6 is the biggest for N = 1,000 (N = 300, …, 1,200) samples. These results suggest that N = 1,000 and m = 6 are appropriate for PE calculation in this study; the optimal selection could be discussed in the future work.

To show the robustness of PE, a simulation study was created by following reference 40. Given a gaussian white noise series of 12,000 samples at the sample rate of 100 Hz (mean = 0, SD = 1), the following artifacts were added to generate a pseudo-electroencephalographic series: (1) a sine wave of 2 Hz (similar to eye movement) at the duration of 10–40 s, (2) a square wave at the duration of 60–90 s, and (3) an impulse at 100 s. The same parameters were selected for calculation of PE and AE. As shown in figure 9, the PE holds steady at approximately 0.93, but the AE varies from 1.45 to 1.69. It is clear that the AE is less sensitive to the artifacts.