EEG entropy measures in anesthesia

Highlights: ► Twelve entropy indices were systematically compared in monitoring depth of anesthesia and detecting burst suppression. ► Renyi permutation entropy performed best in tracking EEG changes associated with different anesthesia states. ► Approximate Entropy and Sample Entropy performed best in detecting burst suppression. Objective: Entropy algorithms have been widely used in analyzing EEG signals during anesthesia. However, a systematic comparison of these entropy algorithms in assessing anesthesia drugs' effect is lacking. In this study, we compare the capability of 12 entropy indices for monitoring depth of anesthesia (DoA) and detecting the burst suppression pattern (BSP), in anesthesia induced by GABAergic agents. Methods: Twelve indices were investigated, namely Response Entropy (RE) and State entropy (SE), three wavelet entropy (WE) measures [Shannon WE (SWE), Tsallis WE (TWE), and Renyi WE (RWE)], Hilbert-Huang spectral entropy (HHSE), approximate entropy (ApEn), sample entropy (SampEn), Fuzzy entropy, and three permutation entropy (PE) measures [Shannon PE (SPE), Tsallis PE (TPE) and Renyi PE (RPE)]. Two EEG data sets from sevoflurane-induced and isoflurane-induced anesthesia respectively were selected to assess the capability of each entropy index in DoA monitoring and BSP detection. To validate the effectiveness of these entropy algorithms, pharmacokinetic/pharmacodynamic (PK/PD) modeling and prediction probability (Pk) analysis were applied. The multifractal detrended fluctuation analysis (MDFA) as a non-entropy measure was compared. Results: All the entropy and MDFA indices could track the changes in EEG pattern during different anesthesia states. Three PE measures outperformed the other entropy indices, with less baseline variability, higher coefficient of determination (R2) and prediction probability, and RPE performed best; ApEn and SampEn discriminated BSP best. Additionally, these entropy measures showed an advantage in computation efficiency compared with MDFA. Conclusion: Each entropy index has its advantages and disadvantages in estimating DoA. Overall, it is suggested that the RPE index was a superior measure. Investigating the advantages and disadvantages of these entropy indices could help improve current clinical indices for monitoring DoA.


INTRODUCTION
In the operating room, general anesthesia is important to guarantee successful surgery and ensure patients' safety and comfort. For anesthesia, the reliable monitoring of anesthetic drug effects on the brain is a clinical concern for anesthesiologists (Monk et al., 2005). The central nervous system (CNS) is the main target of anesthetic drugs. Originated in CNS, the electroencephalogram (EEG) reflects the neural activities of brain, and has been widely used as a surrogate parameter to quantify the anesthetic drug effect (Rampil, 1998;Bruhn et al., 2006;Jameson and Sloan, 2006). However, only limited information can be obtained from the EEG signals purely by waveform observation. With the development of signal processing, various methods have been applied to analyze, identify or detect mental disorders and consciousness mechanisms from EEG signals (Okogbaa et al., 1994;Natarajan et al., 2004;Abásolo et al., 2006), as well as evaluating the effects of anesthesia.
In recent decades, numerous attempts have been made to develop an index for describing anesthetic drug effects on the brain, including zero crossing frequency, spectral edge, wavelet analysis, high-order spectral analysis etc. These studies laid the foundation of commercial EEG-based monitors of depth of anesthesia (DoA), such as BIS (Aspect Medical Systems, Newton, MA) (Bruhn et al., 2006;Ellerkmann et al., 2010) and M-entropy (GE Healthcare, Helsinki, Finland) (Viertiö-Oja et al., 2004;Bruhn et al., 2006). Many of these methods are derived from linear theories. However, various studies have shown that the EEG is a non-stationary signal that exhibits non-linear or chaotic behaviors (Elbert et al., 1994;Pritchard et al., 1995;Zhang et al., 2001;Natarajan et al., 2004). This prompted many researchers to adopt non-linear analysis methods in anesthesia study, for example largest Lyapunov exponent (Fell et al., 1996), Hurst exponent (Alvarez-Ramirez et al., 2008), fractal analysis (Klonowski et al., 2006;Gifani et al., 2007;Liang et al., 2012), detrended fluctuation analysis (DFA) (Jospin et al., 2007;Nguyen-Ky et al., 2010b), recurrence analysis (Huang et al., 2006), and non-linear entropies (Bruhn et al., 2001;Li et al., 2008a). In particular, non-linear entropy methods describing the complexity of EEG signals, have received considerable attention.
Spectral Entropy is the method applied in the commercial M-Entropy Module (Viertiö-Oja et al., 2004). It consists of two parameters: Response Entropy (RE) and State Entropy (SE). SE primarily includes the spectrum of the EEG signal from 0.8 to 32 Hz, and RE includes electromyogram activity from 0.8 to 47 Hz (Viertiö-Oja et al., 2004). Shannon Wavelet entropy (SWE) is the Shannon entropy in the wavelet domain, which indicates signal variation at each frequency scale (Rosso et al., 2001). And the Hilbert-Huang spectral entropy (HHSE) is the Shannon entropy based on the Hilbert-Huang transform proposed by Huang et al. (1998). HHSE has been successfully applied to the anesthetic EEG signals (Li et al., 2008b).
The above methods are based on the frequency spectrum. Whereas many entropy methods are based on the time series and phase space analysis. ApEn is an algorithm derived from the Kolmogorov-Sinai entropy (Pincus, 1991). It quantifies the predictability of subsequent amplitude values of a signal. A previous investigation showed that ApEn correlates well with the concentration of desflurane (Bruhn et al., 2000). However, ApEn lacks relative consistency and is highly dependent on data length, SampEn was proposed to overcome ApEn's limitation by removing self-matching and relieving its bias (Richman and Moorman, 2000). SampEn has been used for analyzing EEG signals (Montirosso et al., 2010;Yoo et al., 2012). FuzzyEn was proposed by Chen et al. (2007). It is based on the fuzzy membership functions to define the vectors' similarity, using the soft and continuous boundaries of fuzzy functions to ensure the continuity and the validity of FuzzyEn's definition (Chen et al., 2009). SPE was introduced by Bandt and Pompe (2002). It is a complexity measure based on symbolic dynamics (Bandt and Pompe, 2002). Because of its simple concept and fast computation, SPE has been widely used in EEG signal analysis (Cao et al., 2004;Li et al., 2007Li et al., , 2008a. Furthermore, its derivatives, multi-scale permutation entropy  and composite permutation entropy index (Olofsen et al., 2008) have been successfully applied to analyze EEG signals during anesthesia.
However, "No one knows what entropy really is, so in a debate you will always have the advantage." This statement is true for EEG analysis today (Ferenets et al., 2006). Each entropy index has its own advantages and disadvantages, but how does their performance compare when evaluating the effect of anesthesia on brain activity? To this end, some researchers have compared the performance of different entropy methods for anesthesia monitoring (Sleigh et al., 2001(Sleigh et al., , 2005Bein, 2006). Unfortunately, these articles analyzed no more than three entropies. To our knowledge, a systematic comparison of the performance of them in assessing anesthesia drug effect is lacking. In this study, we aim to compare the capability of several commonly used entropy indices for monitoring DoA.
We noticed that definitions of all the above entropies are based on Shannon information theory, which belongs to a shortrange or extensive concept. However, the physical systems especially the biomedical systems are often characterized by either long-range interactions, long-term memories, or multifractality (Zunino et al., 2008). To describe these characters, two generalized forms of entropy were proposed: Renyi entropy (Renyi, 1970) and Tsallis entropy (q-entropy) (Tsallis et al., 1998). For example Tsallis entropy has a parameter q for non-extensity. If q > 1, the entropy is more sensitive to events that occur often, whereas if 0 < q < 1 it is more sensitive to the events that occur seldom (Maszczyk and Duch, 2008). In the limit q → 1, it coincides with Shannon entropy. These generalized entropies can provide additional informational about the importance of specific events, such as outliers or rare events. The two classes of entropies and their combinations with current signal processing methods have been already applied in EEG analysis Tong et al., 2003;Inuso et al., 2007) and often been proved advantageous than the Shannon version (Zunino et al., 2008;Arefian et al., 2009). To make the research more instructive, we believe it useful to investigate these non-extensive entropy measures along with those extensive Shannon entropies in DoA monitoring. In this study, we involved the Tsallis wavelet entropy (TWE) and Renyi wavelet entropy (RWE) proposed by Rosso et al. (2003Rosso et al. ( , 2006, as well as the Tsallis permutation entropy (TPE) proposed by Zunino et al. (2008) and a new Renyi permutation entropy (RPE).
In this work, their performance for monitoring DoA were compared. Using data sets obtained during sevoflurane and isoflurane anesthesia, we quantified for each index the responsiveness to loss of consciousness, computation complexity and the ability to detect BSP. Pharmacokinetic/pharmacodynamic (PK/PD) modeling and prediction probability statistics were applied to evaluate the efficiency of each index for tracking anesthetic concentration. Additionally, in order to prove the efficiency of the entropy approaches, two non-linear dynamic methods: DFA (Jospin et al., 2007) and multifractal DFA (MDFA) (Kantelhardt et al., 2002) are compared.

ENTROPY INDICES
The computation of each entropy index is briefly described as follows.

SPECTRAL ENTROPY (RE AND SE)
Spectral Entropy quantifies the probability density function (PDF) of the signal power spectrum in the frequency domain. The detail of the Spectral Entropy algorithm can be seen in Inouye et al. (1991) and Rezek and Roberts (1998). Spectral Entropy consists of the RE and the SE. RE is computed over a frequency range from 0.8 to 47 Hz while SE is computed over the frequency range from 0.8 to 32 Hz. The normalization step for RE and SE are defined as follows: where H sp 0.8−47 and H sp 0.8−32 means the sum of spectral power between 0.8 and 47 Hz, and 0.8 to 32 Hz, respectively. And N 0.8−47 equals the total number of frequency components in the range 0.8-47 Hz. Spectral Entropy describes the degree of skewness in the frequency distribution. For example, in the normalized case, the Spectral Entropy of a pure sine wave with a single spectral peak is 0, while that of white noise is 1.

WAVELET ENTROPY (SWE, TWE, AND RWE)
WE differentiates specific brain states under spontaneous or stimulus-related conditions and recognizes the time localizations of a dynamic process. To calculate Wavelet Entropy, wavelet energy E j of a signal is determined at each scale j as follows: where k and L j are the summation index and the number of coefficients at each scale j with in a given epoch, respectively. The total energy over all scales is obtained by: Then wavelet energy is divided by total energy to obtain the relative wavelet energy at each scale j: SWE is calculated from Shannon entropy of p j distribution between scales as follows: The detail of the algorithm used in this study can be seen in Särkelä et al. (2007).
And the TWE is defined as, where q is a non-extensity parameter.
Based on the definition of Renyi entropy (Renyi, 1970), the RWE is defined as Rosso et al. (2006): For S (S) q , the normalized SWE is where N J is the number of wavelet resolution levels. And defined by Rosso et al. (2003): Further, the normalized S (R) a is defined as Maszczyk and Duch (2008): The values of three WE measures depend on the wavelet basis function, the number of decomposed layers (n) and the data Frontiers in Computational Neuroscience www.frontiersin.org February 2015 | Volume 9 | Article 16 | 3 length (N). Furthermore, the TWE and RWE are related to the parameters q and a respectively. Among these parameters, the wavelet basis function is most important. Because of the lack of a fixed criterion, it is very difficult to select an appropriate wavelet basis function in practical applications and many studies choose it based on experiments. The details of the selection process in this study can be found in Supplement Material 1.

HILBERT-HUANG SPECTRAL ENTROPY (HHSE)
HHSE is based on the Hilbert-Huang transform, which applies the Shannon entropy concept to the Hilbert-Huang spectrum. The detail of the algorithm is seen in Li et al. (2008b). For a given non-stationary signal x(t), the EMD method decomposes the signal into a series of intrinsic mode functions (IMFs), C n (1, 2, . . . , M), where M is the number of IMFs. The signal x(t) can be written by: Apply the Hilbert transform to the IMF components, where ω (t) and a(t) are the instantaneous frequency and amplitude, respectively, of the IMFs. The Hilbert-Huang marginal spectrum is defined by: To simplify the representation, the Hilbert-Huang spectrum is denoted as a function of frequency (f ) instead of angular frequency (ω). The marginal spectrum is normalized by: Next, the Shannon entropy concept is applied to the Hilbert-Huang spectrum, and Hilbert-Huang spectral entropy is obtained by: The HHSE values are mainly affected by the frequency resolution and data length (N). For accurate computation, the frequency resolution is chosen as 0.1 Hz. N directly influences the EMD. In general, the boundary effect may be induced if N is too large or too small, which can contaminate the data and distort the power spectrum. The selection of N in this study is given in Supplement Material 1.

APPROXIMATE ENTROPY (ApEn)
ApEn is derived from Kolmogorov entropy. It was introduced by Pincus (1991). It can be used to analyze a finite length signal and describe its unpredictability or randomness. Its computation involves embedding the signal into the phase space and estimating the rate of increment in the number of phase space patterns within a predefined value r, when the embedding dimension of phase space increases from m to m + 1. For a time series x (i), 1 ≤ i ≤ N of finite length N, reconstitute the N − m + 1 vectors X m (i) following the form: where m is the embedding dimension. Let C m i (r) be the probability that any vector X m (j) is within distance r of X m (i), defined as: where d is the distance between the vectors X m (i) and X m j , defined as: and is the Heaviside function.
After that, define a parameter m (r): Next, when the dimension changes to m + 1, the above process is repeated.
Finally, the approximate entropy is defined by: The detailed algorithm is seen in Bruhn et al. (2000). The ApEn index is influenced by data length (N), tolerance (r) and embedding dimension (m). According to Pincus (1991) and Bruhn et al. (2000), N is recommended to be 1000, r 0.1∼0.25 of the standard deviation of the signal and m 2∼3. The selection of these parameters is described in Supplement Material 1.

SAMPLE ENTROPY (SampEn)
The SampEn proposed by Richman and Moorman (2000) is based on ApEn but differs from it in three ways to remove bias: (1) SampEn eliminates self-matches; (2) To avoid ln 0 caused by removing self-matches, SampEn computes the additional operation of the total number of template well-matches prior to the logarithmic operation. (3) In order to have an equal number of patterns for both embedding dimension m and m + 1, the time series reconstitution in SampEn have N − m rows instead of N − m + 1 in ApEn in embedding dimension m.
The first step of calculating SampEn is the same as ApEn. When the embedding dimension is m, the total number of template matches is: Similarly, when the embedding dimension is m + 1, the total number of template matches is: Finally, the SampEn of the time series is estimated by: SampEn is based on ApEn, so its parameter selection procedure is similar to that of ApEn (see Supplement Material 1).

FUZZY ENTROPY (FuzzyEn)
Zadeh introduced the concept of "fuzzy set" (Zadeh, 1965). Fuzzy set provides a mechanism for measuring the degree to which a pattern belongs to a given class, by introducing the concept of "membership degree" having a fuzzy function u c (x). The nearer the value u c (x) is to unity, the higher the membership grade of x in the set C will be. Inspired by this, Chen et al. (2007) developed the FuzzyEn based on SampEn. FuzzyEn uses the fuzzy membership function u(d m ij , r) to obtain the similarity between X m i and X m j instead of the Heaviside function. FuzzyEn is based on SampEn, so its parameter selection is similar to that of SampEn (see Supplement Material 1).

PERMUTATION ENTROPY (SPE, TPE, AND RPE)
There are three types of PE measures involved in this study. PE is an ordinal analysis method, in which a given time series is divided into a series of ordinal patterns for describing the order relations between the present and a fixed number of equidistant past values (Bandt, 2005). The advantage of this method is its simplicity, robustness and low computational complexity (Li et al., 2007).
For an N-point normalized time series {x(i) : 1 ≤ i ≤ N}, firstly the time series is reconstructed: where τ is the time delay, m is the embedding dimension. Then, rearrange X i in an increasing order: There are m! permutations for m dimensions. Each vector X i can be mapped to one of the m! permutations. Next, the probability of the jth permutation occurring p j can be defined as: where n j is the number of times the jth permutation occurs. Based on the probability of the jth permutation p j , we define SPE, TPE and RPE as follows.
SPE is just the Shannon entropy associated with the probability distribution p j : And the normalized SPE is: Based on the definition of Tsallis entropy, Zunino et al., proposed the normalized TPE and defined it as Zunino et al. (2008): Furthermore, the normalized RPE measure based on the Renyi entropy and permutation probability distribution p j is: In Li et al. (2008aLi et al. ( , 2010, SPE was used to evaluate the effect of sevoflurane and isoflurane anesthesia on the brain. In this study, the parameters of m = 6 and τ = 1 are selected for sevoflurane anesthesia as proposed in Li et al. (2008a). The SPE's parameters for isoflurane anesthesia are the same as those proposed by . TPE and RPE are first used in DoA measure, therefore selection of the appropriate parameters of TPE and RPE should be based on the experiments. The details of the selection process is shown in Supplement Material 1.

EEG data set during sevoflurane-induced anesthesia
In this study, the first data set we used was from a previous study (McKay et al., 2006), in which 19 patients aged 18-63 years were recruited from Waikato Hospital, Hamilton, New Zealand. The subjects were scheduled for elective gynecologic, general, or orthopedic surgery. All patients fasted for at least 6 h before anesthesia and received no premedication. Patients were American Society of Anesthesiologists physical status I or II and signed written informed consent following approval by the Waikato Hospital ethics committee. Before application of Ag/AgCl electrodes, the skin was carefully cleaned with an alcohol swab to ensure electrode-skin impedance of less than 7.5 k . A composite electrode, the Entropy™ Sensor, composed of a self-adhering flexible band holding three electrodes were used to record the EEG signals between the forehead and temple (active = FpZ, earth = Fp1, and reference = F8). RE and SE were measured every 5 s with a plug-in M-Entropy S/5 Module (Datex-Ohmeda). The sevoflurane concentration was measured at the mouth at 100/s (McKay et al., 2006). All data were recorded and stored on a laptop computer. Off-line analysis was performed using the MATLAB (version 8, MathWorks Inc.) software.

EEG data set during isoflurane-induced anesthesia
The second data set contains 29 patients (9 men and 20 women, age 33-77 year) receiving elective abdominal surgery during combined isoflurane general anesthesia and epidural anesthesia (Hagihira et al., 2002). These patients had no neurologic or psychiatric disorders and didn't receive medication with any drugs known to influence anesthesia. The data recordings were approved by the Osaka Prefectural Habikino Hospital and all patients gave written informed consent.
Each patient was injected intramuscularly with 0.5 mg atropine before entering the operating room. Initially, an epidural catheter was placed at the appropriate spinal location. Then, after confirming the effect of epidural analgesia, 3 mg/kg thiopental was used to induce anesthesia. Anesthesia was subsequently maintained with isoflurane, oxygen, and nitrogen after tracheal intubation. Vecuronium was given as required. Lidocaine 1% (80-110 mg/h; initial dose, 90-100 mg) was administered epidurally. Patients received controlled ventilation to maintain adequate oxygenation and normocapnia. To keep mean blood pressure at 60 mmHg, dopamines were administered as required at a dose of 2-5 µg/(kg·min).
Before induction of anesthesia, five EEG electrodes (A1, A2, FP1, FP2, and FPz) were attached to the patients according to the International 10-20 System. FPz was used as the ground electrode. The EEG signal used was recorded from a unipolar lead (FP1-A1) through a 514 X-2 EEG telemetry system (GE Marquette, Tokyo, Japan) with sample frequency of 512 Hz (another Fp2-A2 channel was not analyzed). Isoflurane was initially increased to 1.5% and then stepped down to 0.7%. The end-tidal concentration of isoflurane was purposely maintained at set levels (1.5, 1.3, 1.1, 0.9, and 0.7%) for 30 min at each level. The EEG recordings at 0.3 and 0.5% isoflurane were collected immediately after the operation. The concentration of isoflurane was continuously monitored and recorded by Canomac (Datex, Helsinki, Finland). The BSP was evident in six of the 29 EEG recordings.
The two data sets used can be obtained by asking the authors of corresponding original papers.

EEG PREPROCESSING
All the EEG recordings were preprocessed by following the steps outlined in Li et al. (2010) before further analysis. Firstly, data points whose amplitude values exceeded a threshold determined by mean and standard deviation (SD) statistics were removed as outliers. Then, the filter function filter.m was used to remove the frequency components higher than 60 Hz. This FIR filter ensures that phase information is not distorted. Thirdly the stationary wavelet transform was used to reduce electro-oculogram (EOG) artifact. Finally, an inverse filter was used to detect and remove EMG and other high-amplitude transient artifacts.

PHARMACOKINETIC/PHARMACODYNAMIC MODELING
To derive the relationship between effect-site anesthetic drug concentration and the measured EEG index, PK/PD modeling was used. These methods have been successfully used to evaluate the proposed EEG indices (Li et al., 2008a;Olofsen et al., 2008). It describes the relationship between drug dose and its effect through two successive physiological processes (McKay et al., 2006). The pharmacokinetic (PK) side of the model describes the changes in blood concentration of the drug over time, while the pharmacodynamic (PD) aspect shows the relation between the concentration of drug at its effect site and its measured effect. The simplest effect site model is a first order model, defined as: where C eff denotes the effect-site concentration, k eo is the firstorder rate constant for efflux from the effect compartment and C et is the end-tidal concentration.
In addition, a non-linear inhibitory sigmoid E max model was used to describe the relationship between the estimated C eff and the measured EEG indices.
where Effect is the processed EEG measure, E max and E min respectively are the maximum and minimum Effect for each individual, EC γ 50 is the drug concentration that causes 50% of the maximum Effect and γ is the slope of the concentration-response relationship.
The coefficient of determination R 2 is calculated by: where y i is the measured Effect for a given time andŷ i is corresponding modeled Effect. C eff is estimated by iteratively running the above model with a series of k eo values, with the optimal k eo yielding the greatest R 2 for each patient. Kantelhardt et al., proposed the MDFA method to describe the non-stationary time series, which is based on a generalization DFA method (Kantelhardt et al., 2002 For a time series x(t) of length N, the main computation procedure of MDFA consists of three steps.

MDFA EXPONENT
Step 1. Construct the profile as the equation showed below, where x represents the average value of x(t).
Step 2. Divide the new profile y j into N s = N/s nonoverlapping segments of equal length s. Since the record length N may not be a multiple of the considered time scale s, a short part at the end of the profile will remain in most cases. In order not to disregard this part of record, the same procedure is repeated starting from the other end of the profile y j . Thus, 2N s segments are obtained altogether.
Step 3. Calculate the local trend for each segment by a leastsquare fit of the data and calculate the variance F 2 (s, v). Thus, the qth order fluctuation function is calculated as follows: If q = 0, then It is obvious that when q = 2, we have the standard DFA procedure. MFDFA characterizes the evolution of F q (s) is a function of the segment length s. Modeling fluctuations that present a powerlaw behavior between F q (s) and s, F q (s) ∝ s h(q) , where the h(q) is generalized Hurst exponent.
For the multifractal time series, the scaling behavior is sensitive with the parameter q. For positive q, h(q) describes the scaling behavior of the segments with large fluctuations. On the contrary, for negative q, h(q) is sensitive to small fluctuations. For more detail of the MDFA method, see in Kantelhardt et al. (2002).
In this study, we only considered the influence of q with the MDFA measure. The selection of parameter is described in Supplement Material 1.

STATISTICAL ANALYSIS
To further evaluate the correlation between the measured EEG index and underlying anesthetic drug effect, prediction probability (P k ) statistics were applied, as described in Smith et al. (1996). Given two random data points with different C eff , P k describes the probability that the measured EEG index correctly predicts the C eff of the two points. Its definition is: where P c , P d and P tx separate the probability that two data points drawn at random, independently and with replacement from the population are a concordance, a discordance or an x-only tie. A value of 1 means that the EEG index is perfectly concordant with C eff ; whereas a value of 0.5 means the EEG index is obtained by chance. When the monotonic relation between the drug concentration and the EEG index is negative, the resultant P k value is replaced by 1 − P k . In addition, The Kolmogorov-Smirnov test was used to determine whether the data sets were normally distributed. To assess the index stability during the awake state and the sensitivity to the induction process, the relative coefficient of variation (CV) (Li et al., 2008a) was used. Kruskal-Wallis test was used to determine the significant difference of the index values between awake, induction, anesthesia and recovery states.

RESULTS
First we used these entropy measures on EEG data from sevoflurane anesthesia. Figure 1A shows a preprocessed EEG recording and the derivative from the EEG signal during the whole sevoflurane induction process, from awake to induction, then to deep anesthesia, and finally to recovery. With deepening anesthesia, the mean amplitude of the EEG gradually increased and then the amplitude decreased in the state of recovery. The concurrent endtidal sevoflurane concentration is represented by the black line given in Figure 1B. It can be regarded as the drug concentration in blood, derived from the recorded sevoflurane concentration at the mouth (represented by gray line). The changes in RE, SE, SWE, TWE, RWE, HHSE, ApEn, SampEn, FuzzyEn, SPE, TPE, RPE, and MDFA corresponding to the EEG recording are successively given in Figures 1C-K. As can be seen, all the entropy indices generally followed the changes in EEG pattern as the drug concentration increased and decreased. And MDFA had the opposite trend with entropy indices.
Then we analyzed the EEG recording during isoflurane anesthesia using the same entropy algorithms and MDFA methods. Figures 2A,B are the EEG recording and isoflurane end-tidal concentration respectively. It can be seen that the drug concentration increased and then decreased. Loss of consciousness (LOC) is the most important clinical time point during anesthesia. We investigated the ability of these entropies in tracking LOC. Figure 3 demonstrates the changes in each index around LOC, from LOC−30 s to LOC+30 s for all subjects during sevoflurane anesthesia. For these plots, index values were normalized to between 0 and 1. It can be seen in Figures 3A-N that MDFA(−8) decreased most rapidly, followed by SWE. Thus, the MDFA with q = −8 appeared to be the most sensitive to LOC. To verify this, we calculated the absolute slope values (mean ± SD) of the linear-fitted polynomials vs. time for these indices, as shown in Figure 3O. As can be seen, the absolute slope value for MDFA(−8) (0.44 ± 0.22) is largest, followed by SWE (0.43 ± 0.23).
To further compare the ability of the indices to distinguish different anesthesia states, the sevoflurane anesthesia procedure was divided into four states, i.e., awake, induction, deep anesthesia, and recovery. For each index, a box plot is given in Figure 4. The data was not normally distributed, so the statistics of the 19 patients undergoing sevoflurane anesthesia were expressed as median (min-max), as shown in Table 1. All the entropy indices monotonically decreased as anesthesia deepened, then increased

FIGURE 4 | Box plots of RE, SE, SWE, TWE, RWE, HHEn, ApEn, SampEn, FuzzyEn, SPE, TPE, RPE, MDFA(2) and MDFA(-8) (A-N) at awake (I), induction (II), deep anesthesia (III) and recovery (IV) states.
during recovery. The MDFA indices have an opposite trend with the entropy measures. These are consistent with the results in Figure 1. The overlap of three types of PE (SPE, TPE, and RPE) values between the awake and deep anesthesia states were smaller than the other indices. This means the PE has a better ability to separate these states and a greater robustness for individual differences.
To estimate the baseline variability and the sensitivity to the induction process of each index, the CV value of all the indices for the sevoflurane data set are computed and the results are given in Table 2. During the awake state, the CV value of SampEn was 0.095, which was the highest; The CV value of TPE was 0.003, significantly lower than MDFA (2) (2)  In order to verify the performance of all the indices for monitoring DoA and detecting the burst suppression state, we analyzed the isoflurane anesthesia data set, in which some subjects entered into the burst suppression state during deep anesthesia. The results are given in histogram form and shown in Figure 5. All the indices except SE and MDFA decreased with increasing isoflurane concentration. During burst suppression, only ApEn and SampEn continued to decrease. This means that the ApEn and SampEn algorithms could be used to evaluate DoA including detection of the burst suppression state, without the need for Supplementary Methods. The tabulated results for each index at the different isoflurane concentrations and BSP are presented in Table 3. The CV of the indices show that PE (0.033) outperformed the others in awake state (0% concentration) (see Table 4). And the CV of two MDFA measures were relative higher in awake state. It indicate that MDFA algorithms were no better than some entropy measures in anti-noise performance.
To further compare the performance of the studied indices, PK/PD modeling was performed to describe the relationship between the index values and the estimated sevoflurane and isoflurane effect-site concentration. Tables 5, 6 give these parameters for isoflurane and sevoflurane anesthesia respectively, in which the maximum coefficient of determination (R 2 ) gives the correlation between the index values and the anesthetic effect site concentration. Figures 6A,B show the R 2 values of the indices for the two data sets. Figure 6A shows the R 2 values for sevoflurane. It can be seen that R 2 for TPE (0.95, 95% confidence interval 0.92-0.98) was significantly higher than the other entropy indices. Figure 6B shows R 2 values for isoflurane. Again, R 2 for SPE (0.81) was higher than the other entropy indices. Although R 2 of MDFA with q = 8 was relative higher in sevoflurane anesthesia, the value in isoflurane anesthesia was lower. The statistical analysis also shows that for the same entropy algorithm, the mean R 2 value for sevoflurane was significantly higher than for isoflurane.
To assess the performance of the indices to correctly predict drug effect-site concentrations, we evaluated the prediction probability P k of all the indices from the PK/PD modeling for all the subjects, as shown in Figures 7A,B. And the statistical results are shown in Table 7. Overall, most P k values of indices for sevoflurane were higher than for isoflurane. For sevoflurane, P k of RPE and MDFA were equal (0.87, 95% confidence interval is 0.83-0.90 and 0.83-0.92 respectively), slightly higher than RWE (0.85) and TWE 0.81 (95% confidence interval 0.79-0.84). Also, P k of    RPE was higher than that of TPE and SPE. Similarly, P k of RWE was highest in three WE methods. It means that Renyi entropy had a better performance in predicting drug effect-site concentrations comparing with Shannon entropy and Tsallis entropy. The differences between RPE and the other indices were statistically significant (all p < 0.05, paired t-test), except for MDFA (-8). And the difference between RPE and TPE, SPE were statistically significant (p = 0.03 and 0.01 respectively, paired t-test), which means that RPE had a stronger ability to track the sevoflurane effect-site concentration during anesthesia. In order to get a more intuitive comparison, the best curve fits of all indices against the effect-site concentration are demonstrated for both sevoflurane (Figure 8) and isoflurane (Figure 9). To compare the timeliness performance of each index in tracking DoA, we recorded the computing time of each index for the same subject. 20 EEG recordings from the two data sets were selected. The calculate epoch length (N) of each algorithm is equal to 10 s, and the overlap select 5.0 s. The computing time for 1 min of EEG data compared for each index is given in Table 8. The fastest index was WE (0.025 ± 0.001 s). The RE/SE and PE computation times were 0.096 ± 0.008 s and 0.545 ± 0.016 s respectively. The MDFA (16.338 ± 0.280 s) was the slowest.

Frontiers in Computational Neuroscience
www.frontiersin.org February 2015 | Volume 9 | Article 16 | 11 The desktop computer used for this test had the following configuration: Intel Core i3 CPU, 4 cores at 2.93 GHz, with 2 GB of RAM, running Windows XP professional operating system.

DISCUSSION AND CONCLUSION
In this study, we investigated the performance of 12 entropy algorithms to assess the effect of GABAergic anesthetic agents on EEG activity, including RE, SE, SWE, TWE, RWE, HHSE, ApEn, SampEn, FuzzyEn, SPE, TPE, and RPE. Two data sets including sevoflurane and isoflurane anesthesia were employed as the test samples for evaluating the entropy algorithms. We compared their performance in estimating the DoA and detecting the burst suppression pattern. PK/PD modeling and prediction probability statistics were applied to assess their effectiveness. In addition, we compared the MDFA measure with all entropy indices to test the efficiency of entropy approach. The twelve entropy measures could be divided into two classes: time-domain-based and time-frequency-domain-based analyses. On one hand, ApEn, SampEn, FuzzyEn, and PE are time domain analysis methods. All these entropy algorithms are based on nonlinear theories, and the first three are phase space analytical methods (Chen et al., 2009). PE is based on ordinal pattern analysis of the time series (Bandt, 2005). Considering that the EEG has non-linear characteristics, these four methods have their advantages. For example, FuzzyEn and PE are less sensitive to the signal quality and calculation length (Pincus, 1991;Li et al., 2008a). Relative to ApEn and SampEn, FuzzyEn can resolve more detail in the time series and has more accurate definition in theory (Chen et al., 2009). On the other hand, RE, SE, WE, and HHSE indices are based on the time-frequency domain. The start point of RE and SE is the spectral entropy, which has the particular advantage that the contributions to entropy from any particular frequency range are explicitly separated. In order to achieve optimal response time, RE and SE adopt a variable time window for each particular frequency-called time-frequency balanced spectral entropy (Viertiö-Oja et al., 2004). Compared to the variable time windows of RE and SE, the window function of WE is variable in both time and frequency domains. The HHSE algorithm is based on the EMD and Hilbert transform (Li et al., 2008b). The advantage of this method is that it can estimate the instantaneous amplitude and phase/frequency. Also it can break down a complicated signal without a basis function (such as sine or wavelet functions) into several oscillatory modes that are embedded in this complicated signal. The marginal spectrum gives a more accurate and nearly continuous distribution of EEG energy, which is completely different from the Fourier spectrum (Li et al., 2008b).    Although each entropy algorithm has theoretical advantages with respect to the characterization of EEG recordings during GABAergic anesthesia, we still need to assess the practical performance from several perspectives. In qualitative terms, all the indices are effective at tracking the changes of drug concentration through the EEG analysis. As demonstrated in the presented figures and tables, all the entropies decreased with deepening anesthesia. However, there are quantitative differences between indices for different anesthesia states. This is because the principles underlying each algorithm are entirely different. Entropies based on the time domain, ApEn for example, measure the predictability of future amplitude values of the electroencephalogram based on the knowledge of one or two previous amplitude values. With increasing GABAergic anesthetic drug concentration, the EEG signals become more regular, which leads to a reduction in the ApEn value. Entropies based on the timefrequency domain, such as RE and SE, also decrease with increasing DoA because the EEG shifts to a simpler frequency pattern as the anesthetic dose increases (Rampil, 1998). In all 12 entropy measures, the TWE, RWE, TPE, and RPE are based on the Tsallis entropy and Renyi entropy theory respectively. Tsallis entropy and Renyi entropy theory are considered generalized concept of entropy compared to Shannon entropy. Similar to Renyi entropy, the Tsallis entropy uses the nonextensive parameter q to measure the information of specific events. The results showed that TPE and RPE were better than SPE in assessing the effect of anesthesia. Similar results can also be seen in TWE, RWE, and SWE. There are no studies using TPE or RPE in DoA monitoring before. The excellent performance indicates their potential usefulness in anesthesia analysis. Furthermore, the coefficient of determination and prediction probability statistics were used to assess the correlation of each index with the anesthetic drug effect site concentration. Three PE measures had a higher P k and R 2 compared with the other indices. Also, MDFA at q = 2 had a relative higher P k and R 2 in all indices. Comparing anesthetic drugs, the R 2 values for sevoflurane anesthesia were higher than for isoflurane anesthesia, while the P k values were similar (see Figures 5, 6 and Table 3). This means that the entropy measures were better able to track sevoflurane than isoflurane effect site concentration.
Four additional measures were considered for evaluation of each entropy index. First, the CV was used to evaluate the sensitivity of each index to artifacts during the awake state (Li et al., 2008b. The results showed that PE outperformed the other indices on this level. In all entropy measures, SWE had the highest CV during anesthesia induction, indicating that this index was superior at discriminating between the awake and anesthetized states. Secondly, the performance for estimating the point of LOC was considered. Although all the entropy measures could distinguish between awake and anesthetized states (see Figure 4), the speed of transition (slope) between the two states was fastest  for SWE, while SE had the slowest transition. Thirdly, the performance for discriminating different drug concentrations was considered, especially the ability to distinguish the burst suppression state. The mean ± SD value of the indices showed that all the entropy measures can distinguish different drug concentrations, while only ApEn and SampEn have the ability to distinguish burst suppression from the other states. This means that, if using PE as a DoA index, an additional method for detecting the burst suppression pattern would need to be incorporated, such as Nonlinear Energy Operator (NLEO) (Särkelä et al., 2002). The results are in accordance with the findings during desflurane anesthesia for ApEn (Bruhn et al., 2000) and sevoflurane anesthesia for PE and HHSE (Li et al., 2008b. Finally, the computing time was used to assess algorithm complexity. The results showed that the WE index is the fastest algorithm of all the entropy indices tested. HHSE was the slowest: its computing time for the same data length was about 580 times longer that for WE. In order to improve the computational efficiency, the parallelized method based on the graphics processing unit has been proposed (Chen et al., 2010). The efficiency of these entropy measures were compared with other two non-linear dynamic measures, the MDFA with q = 2 and −8, where MDFA with q = 2 is a standard DFA measure. The results and statistics show that MDFA were better in some aspects compared to some of entropy measures, such as sharper slope in LOC, higher P k and R 2 for sevoflurane (almost equal to RPE) measure. However, there are several shortcomings in MDFA measures. First, CVs of MDFA in awake state were higher compared to those of entropy indices. Second, MDFA could not distinguish the burst suppression state from other states. Most importantly, the computing time of MDFA is the longest in all algorithms, even longer than HHSE, which means that MDFA algorithms are not suitable for real time DoA monitoring. Therefore, entropy approaches are capable for monitoring the EEG changes in anesthesia, and are often advantageous in computation efficiency.
Although this study covers a number of entropy methods and two types of anesthesia, the research has its limitations. For instance, errors caused by individual variability, e.g., age, physical wellness, intraoperative tolerance are hard to control because of the difficulty in data collection in clinical practice. Besides, Interactions between EEG activities and drug concentrations could be studied using finer-grained paradigm, for instance by increasing the drug concentration in a stepwise pattern. Additionally, optimal parameters for each entropy measure may not have been achieved and need further investigation.
This study doesn't provide an absolute measure of "depth" of clinical anesthesia, nor of consciousness for the prevention of intra-operative recall; but rather focuses on understanding the inner workings of each entropy index, and explores whether these indices correlate with GABAergic drug effect. Having a good understanding of the strengths and weaknesses of each measure is necessary before possibly applying them within a clinical context. In conclusion, each entropy measure has its advantages, and several indices show promise as a simple open-source method for quantifying the brain effects of GABAergic drugs. In particular, the PE indices perform better than other entropy indices as an EEG derivative in several aspects, especially for RPE measure. However, further work is required to accurately quantify the burst suppression pattern. Also, to be useful as a clinical measure, each algorithm still needs additional parameter and computation efficiency optimizations.