Sample entropy analysis for the estimating depth of anaesthesia through human EEG signal at different levels of unconsciousness during surgeries

Filtering by MEMD

EEG signal can be easily contaminated by artefacts such as noise from electromyography (EMG), electrooculography (EOG), and electric circuits (Huang et al., 2013). Thus, it is necessary to filter the signals to reduce the computation error. Regarding the nonlinear and nonstationary characteristics of EEG, a method known as the empirical mode decomposition (EMD) was proposed by Huang et al. (1998). The EMD method decomposes signal into a series of IMFs. The extraction procedure of an IMF from a signal is known as sifting. In this study, all the local extrema were found, and then, the upper envelope was generated by drawing a cubic spline line through all the local maxima. Similar steps were conducted for the local minima to produce the lower envelope.

All the data should fall in the range between the upper and lower envelopes. The first component h₁ is obtained by computing the difference between the data and mean: (1)${\displaystyle X\left(t\right)-m_1=h_1}$ where m₁ denotes the mean, and h₁ can only be considered as a proto-IMF. In the next step, h₁ is treated as data: (2)${\displaystyle h_1-m_{11}=h_{11}.}$ After sifting up to k times, h₁ becomes an IMF as follows: (3)${\displaystyle h_{1\left(k-1\right)}-m_{1k}=h_{1k}.}$ Then, the first IMF is acquired as h_1k: (4)${\displaystyle C_1=h_{1k}.}$ One can repeat these steps for the residual generated by previous loops to obtain successive IMFs.

Finally, the signal can be decomposed as: (5)${\displaystyle X\left(t\right)=\sum _{i=1}^N{c}_i\left(t\right)+r_N\left(t\right)}$ where N is the number of IMFs, c_i(t) is the ith IMF and r_N(t) is the residual. A combination of different IMFs can reconstruct a signal to eliminate the artefacts. However, EMD induces the mode mixing problem. An improved EMD known as MEMD was presented by Rehman & Mandic (2009) and noise assisted MEMD (N-A-MEMD) was also introduced by (Ur Rehman & Mandic, 2011). The new N-A-MEMD solves the mode mixing problems, and also can be applied to single channel data by adding Gaussian noise together to constitute multichannel data. The definition is as in Eq. (6). (6)${\displaystyle m\left(t\right)=\frac{1}{K}\sum _{i=1}^K{e}^{\theta i}\left(t\right)}$ where e^θi(t) is the multivariate envelope curves of the whole set of direction vector, K is the length of the vectors, and m(t) is calculated by means of the multivariate envelope curves. Unlike the traditional method that has constant amplitude and frequency in the harmonic component, an IMF provides a simple oscillatory mode as a counterpart and has a variable amplitude and frequency over time. Because decomposition is based on the local characteristics of a data series, it is applicable to nonlinear and nonstationary processes, thus overcoming both pseudo-linearity and stationarity. In our previous study (Huang et al., 2013), the EEG signals were reconstructed by summing IMF2 and IMF3 after decomposition to obtain the filtered signals due to their superior discrimination ability of different states of anaesthesia.

SampEn analysis

Due to the non-stationary characteristic of EEG data, SampEn is appropriate to quantify the degree of irregularity of EEG data. A highly complicated series generates high SampEn values. SampEn, proposed by Richman & Moorman (2000), improves ApEn to more precisely measure the complexity of a physical time series by following steps:

Step 1: Obtain a time series with N points $\left\{X\left(i\right),1\le i\le N\right\}$. Set parameters r and m, which represent tolerance and embedding dimension, respectively.

Step 2: The m-dimension template vector $X_m\left(i\right)$ is defined as: (7)${\displaystyle X_m\left(i\right)=\left\{X\left(i\right),X\left(i+1\right),\dots ,X\left(i+m-1\right)\right\},1\le i\le N-m+1.}$ Step 3: Calculate the distance between two vectors: (8)${\displaystyle d\left[X_m\left(i\right),X_m\left(j\right)\right]=max\left(\left|X_m\left(i+k\right)-X_m\left(j+k\right)\right|\right)0\le k\le m-1,i\ne j.}$ Step 4: Let B_i be the number of vectors $X_m\left(j\right)$ within r of $X_m\left(i\right)$, then: (9)${\displaystyle B_i^m\left(r\right)=\frac{1}{N-m+1}{B}_i}$ (10)${\displaystyle B^m\left(r\right)=\frac{1}{N-m}\sum _{i=1}^{N-m}{B}_i^m\left(r\right).}$ Step 5: If we set m = m + 1 and repeat Step 1 to Step 4, we obtain the following: (11)${\displaystyle A^m\left(r\right)=\frac{1}{N-m}\sum _{i=1}^{N-m}{A}_i^m\left(r\right).}$ Finally, SampEn is obtained as: (12)${\displaystyle \text{SampEn}=-ln\frac{A^m\left(r\right)}{B^m\left(r\right)}.}$ Before computing the SampEn, three crucial parameters must be set: the length of the time series N, tolerance r, and dimension m. N was chosen to be 3,750 points (30 s). Various theoretical and clinical applications have described the validity of m = 1 or 2 and r in the range of 0.1 × std − 0.2 × std, where std is the standard deviation. The algorithm was initialized with the following parameters: m = 1 and r = 0.1 × std.

In addition, permutation entropy and recurrence analysis are utilised to analyse the nonstationary EEG data. Permutation entropy is an ordinal analysis method and is quantified by dividing time series into ordinal patterns for describing the order of relations between the present values and a fixed number of equidistant past values. Researchers claim it offers simplicity, robustness, and low computational complexity (Zanin et al., 2012); these advantages are mostly required in nonstationary data analysis, including studies related to anaesthesia (Li, Cui & Voss, 2008). In this study, we chose dimension m = 6 and lag τ = 1, as suggested in previous studies (Li, Cui & Voss, 2008; Liu et al., 2016). Recurrence analysis is a kind of statistical analysis method used to quantify the frequency of the phase space trajectory of a dynamical system visiting approximately the same area in some phase space (Huang, Wang & Singare, 2006; Marwan et al., 2007). We used the RQA tool to conduct the analysis (Marwan et al., 2007). The parameters were chosen: embedding dimension equal to 3, delay equal to 5, and threshold equal to 0.5, which were decided by the global false nearest neighbour and mutual information. In this study, both recurrence rate (RR) and determinism (DET) were used to measure the recurrence analysis. The data segment window size and step were the same as in the SampEn processing part.

Regression models

In order to symbolise DoA more accurately, three typical regression models were applied to train our datasetdash SVM, ANN and Random Forest. These models serve vital purposes in machine learning and have strong self-learning capabilities. Libsvm was employed to conduct SVM regression (Chang & Lin, 2011; Drucker et al., 1997). The exponential function was chosen as the kernel. From engineering perspective, three to four layers were mostly used for ANN (Kourentzes, Barrow & Crone, 2014; Ripley & Ripley, 2001). The backpropagation neural network (BPNN), a structure of 1–9–18–1, was chosen as the model after trials. The commonly recognised BIS values were regarded as the target whereas the entropy features were employed as input data. The input sample data was randomly divided into training (70%), validation (15%), and testing datasets (15%) to adjust the weights and bias for generating the eventual model. Random Forest is a type of ensemble of trees. We used 1,000 trees to train the SampEn extracted from filtered EEG. For all the training processes, the corresponding BIS acted as a target as well.

Discrimination of the different anaesthetic stages

The awake stage denotes the period before the anaesthetic injection. Patients are notably aware during this stage. Then LoC is slowly induced by drugs, eventually followed by a stable unconsciousness state. By the end of the operation, the patients gradually recover consciousness the as the drug volume decreases. Therefore, it is worthwhile to distinguish the stages prior to tracking DoA for the entire surgery. A 2-min noise-free segment was extracted from the three stages for each patient. By applying a window size of 30 s and a window step of 5 s, 24 data points were calculated for each stage, and then an average of the points of all cases were taken. Figure 5 displayed four indexes percentile distribution of the three stages—RQA-RR, RQA-DET, PEn, and SampEn. All the indexes percentile distributions displayed a significant difference between the unconsciousness and conscious state (p < 0.05). Among them (Figs. 5A and 5D), the RR and SampEn were more robust and compatible to the awake and recovery period. Therefore, these distributions were more capable of processing the consciousness stage data and reducing the risk of classifying the EEG into a wrong stage. Moreover, both distributions were more sensitive to the recovery from anaesthesia. Both RQA results increased during anaesthesia, whereas the entropy ones had a low value at LoC (Figs. 5C and 5D). This finding is consistent with the stable and less complex brain activities during unconsciousness.

Perioperative performance of evaluating DoA

To find the optimal index to symbolise DoA, a further statistical analysis was conducted. By regarding the corresponding perioperative BIS value as the gold standard, the correlation coefficient, MAE, and AUC under ROC were obtained, as listed in Table 2. From the comparisons, we observed that the entropy method had a much better capability than the RQA group. Moreover, SampEn had the highest correlation with BIS (0.72 ± 0.09), and the corresponding AUC was the highest as well (0.95 ± 0.04). The AUC usually ranged from 0.5 to 1. The more effectively a method classified input, the larger the value was. This provided strong evidence for using SampEn to assess DoA in this study.