Assessment of Nociceptive Responsiveness Levels during Sedation-Analgesia by Entropy Analysis of EEG ^†

Abstract

The level of sedation in patients undergoing medical procedures is decided to assure unconsciousness and prevent pain. The monitors of depth of anesthesia, based on the analysis of the electroencephalogram (EEG), have been progressively introduced into the daily practice to provide additional information about the state of the patient. However, the quantification of analgesia still remains an open problem. The purpose of this work was to analyze the capability of prediction of nociceptive responses based on refined multiscale entropy (RMSE) and auto mutual information function (AMIF) applied to EEG signals recorded in 378 patients scheduled to undergo ultrasonographic endoscopy under sedation-analgesia. Two observed categorical responses after the application of painful stimulation were analyzed: the evaluation of the Ramsay Sedation Scale (RSS) after nail bed compression and the presence of gag reflex (GAG) during endoscopy tube insertion. In addition, bispectrum (BIS), heart rate (HR), predicted concentrations of propofol (CeProp) and remifentanil (CeRemi) were annotated with a resolution of 1 s. Results showed that functions based on RMSE, AMIF, HR and CeRemi permitted predicting different stimulation responses during sedation better than BIS.

1. Introduction

Pain in patients should be controlled for physiological reasons as well as for moral, humanitarian and ethical reasons. Therefore, obtaining information about the responsiveness level in patients under minimally invasive procedures, such as endoscopies, is one of the most important requirements for the anesthesiologist in order to guarantee a comfortable state for the patient during the medical procedure. Sedation-analgesia is a continuum that comprises four sedation levels from minimal sedation (anxiolysis) to general anesthesia [1]. The first level corresponds to a drug-induced state during which ventilatory and cardiovascular functions are unaffected and patients can respond normally to verbal commands, although cognitive function and coordination may be impaired. Moderate sedation/analgesia is the next level that relates to a drug-induced depression of consciousness during which ventilation is adequate and cardiovascular function is usually maintained, avoiding any interventions to maintain an open airway. Patients can respond purposefully to verbal commands, either alone or accompanied by light tactile stimulation. The following level is deep sedation/analgesia, a drug-induced depression of consciousness during which cardiovascular function is usually maintained but the ability to independently maintain ventilatory function may be impaired, requiring the assistance in maintaining a patent airway. In this level, patients cannot be easily aroused but respond purposefully following repeated or painful stimulation. General Anesthesia corresponds to the last level, a drug-induced loss of consciousness during which patients are not arousable, even by painful stimulation. This level is characterized by a depressed spontaneous ventilation, requiring assistance in maintaining a clear airway and positive pressure ventilation in the patients, and also the cardiovascular function might be depressed as well as the cardiovascular control.

Adverse drug responses can be detected and treated in a timely manner by monitoring the level of consciousness, thus avoiding many of the complications associated with sedation and analgesia [1]. Several methods have been developed for the noninvasive assessment of the level of consciousness during general anesthesia, which include pharmacodynamic assessment of the anesthetic and analgesic agents, study of hemodynamic evolution during anesthesia and evaluation of the modification of the electroencephalographic (EEG) activity [2,3,4,5,6]. Since the main action of anesthetic agents occurs in the brain, the analysis of the EEG signals is one of the most applied techniques to monitor the sedation level. A central topic would be to assess hypnotic effect and pain/nociception-related activation from EEG. Most of the commercial monitors based on EEG have been developed for measuring the hypnotic effect with methods based on bispectrum (BIS, Aspect Medical System, Newton, MA, USA) [7,8], entropy (M-entropy module, GE Healthcare, Helsinki, Finland) [9], auditory evoked potentials (AEP, Monitor/2, Danmeter, Odense, Denmark) [10,11], or power spectral analysis on different frequency bands (qCON, Quantium Medical, Barcelona, Spain) [12].

Although different methods have been proposed over the last decade to monitoring nociception, this is a topic that has not been completely solved. Most of these methods include the analysis of signals from the brain activity, such as EEG or auditory-evoked potentials (AEP) [10,11], or signals from the cardiovascular activity, such as heart rate variability (HRV) [13], and also consider skin conductance (SC) [14] and combinations of these signals [15]. HRV and SC signals are influenced by the sympathetic activity, being the systems sensitive to other disturbances in addition to pain/nociception like changes in blood pressure or heart rate due to sympathomimetic drug delivery, perioperative bleeding or patient’s baseline condition (hypertension and arrhythmias of diverse etiology). There is not consensus in the relationship between different levels of analgesia and changes in EEG features, and therefore methods based on EEG signals tend to be unspecific. The AEP signals are weakly coupled to the levels of analgesia and its amplitude is very small, being difficult to record in clinical practice without significant noise levels [16]. Recently, studies based on time-frequency representation [17] and auto-mutual information function (AMIF) [18] were applied in order to determine associated changes between EEG spectrum and EEG complexity with the prediction of the levels of unconsciousness as measured via the Ramsay Sedation Scale (RSS). In addition, EEG single-scale and multi-scale entropy measures have been effective in tracking changes of drug concentration during general anesthesia [19,20]. Similarly, a recent study [21] was proposed to measure the information coupling of EEG signals in sevoflurane anesthesia by means of a modification in the AMIF, which included permutation entropy, obtaining a measure that can clearly distinguish the awake and anesthesia states. However the discrimination between sedation-analgesia levels with no nociceptive response or sluggish response to painful stimulation still remains an open problem.

The aim of the study is to evaluate complexity of EEG signals during sedation-analgesia using a multiscale approach in order to improve the prediction of pain responses. We exploited the refined multiscale entropy (RMSE) [22], a refined version of the multiscale entropy (MSE) originally proposed in [23]. RMSE was selected, instead of a more traditional entropy-based approach, because it assesses complexity of times series. This feature is extremely helpful because it is well-known that EEG dynamics is regulated by several mechanisms operating along different frequency bands. Furthermore, AMIF measure is also considered to give information about the complexity in EEG signals based on the information coupling [24]. To evaluate the prediction of response to painful stimulus, a multivariate analysis is performed on RMSE variables, AMIF measures, power spectral measures and other parameters such as drug concentration and heart rate, in order to maximize specificity and sensitivity.

2. Experimental Section

2.1. Database

The data were recorded from 378 patients (mean age 62.6 ± 12.5 years), of which 247 were men (62.4 ± 12.2 years) and 131 were women (63.0 ± 13.2 years), scheduled to undergo ultrasonographic endoscopy (USE) of the upper gastrointestinal tract under sedation-analgesia. Most of the cases were for discarding neoplastic processes in pancreatic level, gastric or duodenal, and there were no mucosectomy. Two diagnostic procedures were applied: simple endoscopy and endoscopy puncture. This last was started with a simple endoscopy for exploring the sonographic anatomy of the lesion followed by transluminal aspiration through the gastrointestinal wall, being the main purposes: (a) staging neoplastic lesion (single ultrasound) since infiltration of certain layers of the intestinal wall and the presence of surrounding adenopathies modifies the stage of tumor process; and (b) obtaining information of the anatomopathologic lesion in order to typify the neoplasia and thus guide therapy, either surgical or chemotherapeutic. USE was chosen because it is a relatively long procedure, approximately 1 h, with periods of stability of effect and others when the intensity of endoscopic stimulus is evident, allowing study of the repercussion of painful stimulus on the level of sedation. All patients belong to 1–3 ASA classification. Patients with altered central nervous system, medicated with analgesics or drugs with central effects on the perception of pain, from moderate to severe cardiomyopathy, neuropathy or hepatopathy that needed control during the anesthetic process were not included in the database. This study was done after receiving approval from the Ethic Committee of Clinical Research of Hospital Clinic de Barcelona and written informed consent from all involved patients.

After arrival in the USE room, every patient was routinely monitored, including electrocardiogram and BIS of the EEG, which was continuously measured using an A2000 monitor (software version 3.31) (Aspect Medical Systems, Newton, MA, USA). A 4-electrode sensor was placed on the forehead of the patients according to the manufacturer’s instructions. The raw EEG signal was recorded and stored into a flash memory card with the AEP monitor/2 (Danmeter, Odense, Denmark), using a montage of 3 electrodes: middle forehead (+), malar bone (−), and left forehead electrode used as reference. The EEG was recorded with a sampling frequency of 900 Hz, with a resolution of 16 bits and a recording time of about 60 min. Propofol and remifentanil were infused using a target controlled infusion (TCI) system (FreseniusVial, Chemin de Fer, Béziers, France). The TCI system administered propofol according to the predictions of the PK-PD model described by [25,26] and remifentanil according to the predictions of the PK-PD model described by [27]. In both cases, the TCI was targeting the effect site.

The RSS score (see

Table 1. Ramsay sedation scale (RSS) [28].

Score	Description
1	Patient awake, anxious, agitated or restless
2	Patient awake, cooperative, orientated and tranquil
3	Patient drowsy with response to commands
4	Patient asleep, brisk response to glabella tap or loud auditory stimulus
5	Patient asleep, sluggish response to stimulus
6	No response to firm nail-bed pressure or other noxious stimuli.

) [28] was evaluated at random times during the procedure in order to avoid the factors correlated with time, which could confound the results of the RSS measurements. The whole database contains annotated RSS scores from 2 to 6. Therefore, two observed categorical responses after the application of painful stimulation were available in the database: the evaluation of the RSS after nail bed compression and the presence of gag reflex (GAG) during endoscopy tube insertion. A GAG value of 1 corresponds to a positive nausea reflex during endoscopy tube insertion, while a GAG value of 0 corresponds to no response after tube insertion.

Therefore, the following information was available for each patient: (a) continuous raw EEG signal; (b) BIS, heart rate (HR), predicted concentrations of propofol (CeProp) and remifentanil (CeRemi), annotated with a resolution of 1 s; and (c) RSS and GAG annotated at random times.

2.2. EEG Signal Preprocessing

EEG signals were sampled at 128 Hz after applying a band-pass finite impulse response (FIR) filter of 100th order with cut-off frequencies of 0.1–45 Hz. Then, the EEG signals were segmented into windows of length of 1 min between 90 s and 30 s before the response annotation of RSS or GAG. Assuming that, in an ideal situation, the sedation level should be constant if plasma concentration of the anesthetic and analgesic agents remains without changes, the annotated RSS was assigned to the previous 1-min length window if the differences in the predicted concentrations of remifentanil (ΔCeRemi) and propofol (ΔCeProp), calculated between the first and the last second of the window, were ΔCeRemi < 0.1 ng/mL and ΔCeProp < 0.1 μg/mL. Otherwise, the window was cut at the sample where the conditions were satisfied. Windows of EEG containing high amplitude peak noise were filtered with a filter based on the analytic signal envelope [29]. Furthermore, if the difference between adjacent samples were higher than 10% of the mean of the differences of the previous samples, the window length was shortened to avoid the artifact. In this way, the smallest window resulted to be of 50 s.

The power spectral density (PSD) for each EEG window was calculated using the Welch method of averaged modified periodograms. In this analysis, EEG segments were divided into one-second segments with 25% of overlap and a Hamming window was applied to each segment. The final spectral density was achieved as the average of spectral densities of all the segments. Then, spectral power in each band (P_δ, 1–4 Hz; P_θ, 4–8 Hz; P_α, 8–12 Hz; P_β, 12–30 Hz) was computed as the area under the normalized PSD curve for the given frequency range.

2.3. Entropy Based Measures

2.3.1. RMSE

Refined multiscale entropy (RMSE) [22,30,31] was applied as a technique of multiscale analysis in which sample entropy (SampEn) [32] is computed at different time scales (ts) of the EEG windows, given information of the regularity of the signal in each one of the time scales. In this way, given a time series x = {x(i), i = 1, …, N}, where i is the sample counter and N is the series length, SampEn estimates the probability that two patterns of length L, x_L(i) = (x(i), x(i − 1), …, x(i − L + 1)) and x_L(j) = (x(j), x(j − 1), …, x(j − L + 1)), that are similar in the L-dimensional phase space remain similar after adding new samples x(i + 1) and x(j + 1), i.e., x_{L + 1}(i + 1) and x_{L + 1}(j + 1) are still similar in the (L + 1)-dimensional phase space. Two patterns are similar if the distance between them is less than a tolerance parameter “r”. In this work, RMSEts represents the sample entropy computed at the time scale ts. RMSE provides two refinements of MSE: (i) it offers a way to improve the procedure devised to eliminate fast time scales, avoiding the aliasing effects; and (ii) it modifies the coarse graining procedure in such a way that reduction of the variance in the EEG, because the elimination of fast time scales, does not tend to artificially decrease entropy rate as a function of the time scale. The result is a more reliable method for the assessment of entropy-based complexity as a function of the temporal scale. The process was iterated with scale factors ranging from 1 to 20 and SampEn was plotted as a function of ts. The procedure to calculate RMSE is similar to that to compute MSE except for two substantial variations: (i) the suboptimal FIR filter that is used in MSE was substituted with an sixth order low-pass Butterworth filter, limiting as much as possible aliasing at any scale factor during undersampling; (ii) the parameter r setting the tolerance for the calculation of the SampEn was continuously adapted with the scale factor ts, thus preventing the dependence of RMSE on the change of variance induced by the procedure of elimination of the temporal scales. In this work, SampEn was calculated with r = 0.15 × SD(ts) and L = 2, where SD(ts) is the standard deviation of the EEG signal at each ts and L is the length of the patterns used to find similitudes in the signal. We make reference to [22] for specific details on the method.

A strategy similar to that set [31] and adopted to EEG analysis was followed. RMSE at short time scale (ts = 1) is equal to assess complexity of the original series, thus taking into account all the time scales present in the EEG signal. The application of the low-pass Butterworth filter at medium time scale in RMSE (ts = 2–4) reduces the superior limit of the pass-band from 32 Hz (ts = 2) to 16 Hz (ts = 4) with a sample frequency of 128 Hz, gradually filtering contributions in the β band, in addition to those at the highest frequencies, and leaving intact fluctuations in the δ, θ and α band. At long time scale, low-pass filtering approach reduces the superior limit of the pass-band from 12.8 Hz (ts = 5) to 3.2 Hz (ts = 20), progressively canceling contributions in α and θ bands, in addition to those in the β band and at highest frequencies, and leaving intact fluctuations slower than θ band. In addition, the slope of the RMSE course at long time scales (5 ≤ ts ≤ tsn, with tsn = 6 to 20) was calculated. Special emphasis was focused in the slope computed between 5 ≤ ts ≤ 8 (RMSE_α, 8–12 Hz), 8 ≤ ts ≤ 16 (RMSE_θ, 4–8 Hz), and 16 ≤ ts ≤ 20 (RMSE_δ, 3.2–4 Hz). These slopes provide information about the increase or decrease of RMSE profile at long time scales.

2.3.2. AMIF

Auto mutual information function (AMIF) is calculated as the mutual information between two measurements x_i and x_{i + τ} taken from a single time series. This function describes how the information of a signal (AMIF value at τ = 0) decreases over a prediction time interval (AMIF values τ > 0). In the case of a completely regular and deterministic signal, the AMIF would repeat a specific pattern for all τ. In the case of an uncorrelated random signal, the AMIF would become zero for all τ apart τ = 0. Therefore, an information loss increasing with τ is related to a gradual decrease of signal predictability and a gradual increase of signal complexity [33]. AMIF can be computed from Rényi entropy as:

$$AMIF\left(\mathsf{\tau }\right)=\frac{1}{q-1}lo{g}_2{\displaystyle \sum }_{x_iϵX}^{}{\displaystyle \sum }_{x_{i+\mathsf{\tau }}ϵX}^{}\frac{{P_{xx}}^q\left(x_i,x_{i+\mathsf{\tau }}\right)}{{P_x}^{q-1}\left(x_i\right)P_x^{q-1}\left(x_{i+\mathsf{\tau }}\right)}$$

(1)

where q is the control parameter that defines Rényi information, which was selected from the set of values {0.1, 0.2, 0.5, 2, 3, 5, 10, 30, 50, 100}. Indeed, the largest probabilities most influence the AMIF when q > 1 and the smallest probabilities most influence the values of AMIF when 0 < q < 1. The AMIF converges to the traditional definition of the AMIF when q tends to 1. The probabilities Pxx and Px were constructed on the series x_i and their delayed series x_{i + τ}, for τ = {1, 2, …, 128} samples. Finally, AMIF was normalized by dividing each value for the AMIF (τ = 0).

In this study, two indexes were defined on the AMIF along the delay τ (Figure 1): the first relative maximum (max) and the decay for τ = 1 (FD). These indexes were calculated from the EEG signal filtered in each one of the characteristic frequency bands (δ, θ, α, and β).

2.3.3. Statistical Analysis

A non-parametric test, U of Mann–Whitney test, was applied and a significance level p-value < 0.05 was set. Variables that satisfy this condition were considered for building univariate and multivariable discriminant functions, in order to predict presence or absence of response to noxious stimuli. In this work, the analysis was focused on the observed categorical responses after nail bed compression (RSS evaluation) and the best univariate and multivariate functions were considered to study GAG reflex during endoscopy tube insertion in the same database.

For the RSS evaluation, the leave-one-out method was performed as validation in two different tests: Trial 1 (2 ≤ RSS ≤ 5 vs. RSS = 6) and Trial 2 (RSS = 5 vs. RSS = 6). Indeed, 2 ≤ RSS ≤ 5 and RSS = 5 correspond to groups with presence of response to noxious stimuli and RSS = 6 to group with absence of response to noxious stimuli. In both trials, it is hypothesized that the unresponsive state is caused by an appropriate interaction between hypnotic and analgesic effects, while the responsive states are characterized by an insufficient presence of analgesic effect, which would be happening especially in Trial 1 where RSS from 2 to 5 are included in group 2 ≤ RSS ≤ 5. The validation analysis was performed taking into account the presence or the absence of gag reflex, GAG = 1 and GAG = 0, respectively. From this statistical analysis, sensitivity (Sen) and specificity (Spe) were computed, where Sen represents responsive states (RSS < 6 and GAG = 1) and Spe represents unresponsive states (RSS = 6 and GAG = 0). In order to increase the values of Sen and Spe, RMSE indexes were combined with AMIF indexes, temporal and spectral HR measures, predicted concentrations of CeRemi and CeProp. In all combinations, a maximum of four uncorrelated indexes were considered.

Additionally, the ability of the indexes to describe pain responses was also evaluated using prediction probability (Pk), which compares the performance of indicators [34]. The Pk coefficient is a statistic commonly used to measure how well an index predicts the state of the patient. A Pk of 1 represents a perfect prediction and 0.5 is not better than tossing a fair coin. The Pk avoids the shortcomings of other measures being independent of scale units and it does not require knowledge of underlying distributions.

3. Results

3.1. CeProp, CeRemi, BIS, Time and Spectral HR Indexes

Table 2. CeProp, CeRemi, BIS, HR and spectral power bands indexes.

Index	Groups (mean ± std)			Trial 1			Trial 2
				(2 ≤ RSS ≤ 5 vs. RSS = 6)			(RSS = 5 vs. RSS = 6)
	2 ≤ RSS ≤ 5	RSS = 5	RSS = 6	Pk	Sen	Spe	Pk	Sen	Spe
CeProp	1.829 ± 0.907 †	2.363 ± 0.705	2.382 ± 0.669	0.693	68.2	58.1	0.502	45.0	53.8
CeRemi	1.106 ± 0.815 †	1.034 ± 0.820 †	1.386 ± 0.598	0.622	56.9	56.6	0.642	58.3	56.8
BIS	76.1 ± 13.6 †	65.2 ± 13.7 †	59.4 ± 14.3	0.799	75.7	68.6	0.622	64.3	56.6
mHR	72.2 ± 13.3 †	69.6 ± 12.4 †	66.7 ± 11.9	0.616	54.0	61.2	0.559	51.0	55.5
sdHR	2.984 ± 3.699 †	2.567 ± 3.490 †	1.995 ± 2.862	0.627	37.2	80.3	0.554	30.8	78.2
Pβ	0.252 ± 0.164 †	0.204 ± 0.125 †	0.144 ± 0.094	0.721	52.5	82.2	0.663	48.2	77.2
Pα	0.300 ± 0.159 †	0.400 ± 0.128 †	0.422 ± 0.147	0.712	63.1	67.4	0.555	46.7	57.6
Pθ	0.165 ± 0.061 †	0.179 ± 0.057 †	0.201 ± 0.062	0.657	63.3	58.1	0.597	57.2	52.6
Pδ	0.313 ± 0.194 †	0.250 ± 0.150 †	0.266 ± 0.142	0.552	43.7	66.0	0.544	66.5	41

^†: p-value < 0.05 vs. RSS = 6 group. Pk: prediction probability; Sen: sensitivity (%); Spe: specificity (%). Values are expressed as mean ± std.

shows the mean and standard deviation (mean ± std) of the predicted concentrations of propofol (CeProp) and remifentanil (CeRemi), BIS parameter, mean heart rate (mHR), variability of the heart rate (sdHR) and spectral power in bands δ, θ, α, and β. Each index was obtained from EEG windows of length of 1 min between 30 s and 90 s before the response annotation of RSS. The values of Pk, Sen and Spe calculated with a univariate linear discriminant function are also indicated in this table for Trials 1 and 2.

It can be noted from Table 2 that the mean value of indexes BIS, mHR, sdHR, P_β was higher (p-value < 0.05) in the responsive groups (2 ≤ RSS ≤ 5 and RSS = 5) than in the unresponsive group (RSS = 6) in both trials. Indexes as CeRemi, P_α and P_θ had a contrary behavior, showing lower mean value (p-value < 0.05) in the responsive groups (2 ≤ RSS ≤ 5 and RSS = 5) than in the unresponsive group (RSS = 6). The index CePropo showed only statistically significant differences in Trial 1, being the mean value lower in the responsive group (2 ≤ RSS ≤ 5) compared with the unresponsive group (RSS = 6). In Trial 1, the P_δ showed higher values in the responsive group (2 ≤ RSS ≤ 5) than responsive group (RSS = 6), while it was lower in the responsive group (RSS = 5) than the unresponsive group (RSS = 6) for Trial 2.

In Trial 1, the BIS was the index with the best Pk (0.799), also presenting the highest Sen (75.7%). The best Spe (82.2%) was obtained with the P_β index, although the sdHR index showed also high Spe (80.3%) but it had the worst Sen (37.2%). In Trial 2, the spectral power P_β was the index with the best Pk (0.663), although CeRemi showed also an important Pk value (0.642). Again, the sdHR index had the best Spe value (78.2%).

3.2. Multiscale Entropy Analysis of EEG: RMSE and AMIF

Mean ± standard error values of the SampEn computed along the time scale factor ts (RMSEts) are shown in Figure 2. These values were obtained from the EEG segments in responsive states (2 ≤ RSS ≤ 5 and RSS = 5 for Trials 1 and 2, respectively) and unresponsive states (RSS = 6 for both Trials). All the RMSE curves exhibited an initial fast increase at short and medium time scale from scale ts = 1 to ts = 4, followed by a slow decrease at long time scales, being the entropy values higher at long time scales (ts = 4 to ts = 20) than at short time scales (ts = 1), in both responsive and unresponsive groups. The comparison between responsive and unresponsive states showed a larger difference for Trial 1 than for Trial 2. In general, entropy values were lower in unresponsive state compared with responsive state at short time scales ts, whereas a contrary behavior was observed at long time scales, where entropy values were higher in the unresponsive group.

RMSE indexes with p-value < 0.001 and Pk > 0.60 in at least one Trial were included in

Table 3. Indexes obtained from multiscale entropy analysis: RMSE and AMIF.

Index	Groups (mean ± std)			Trial 1			Trial 2
				(2 ≤ RSS ≤ 5 vs. RSS = 6)			(RSS = 5 vs. RSS = 6)
	2 ≤ RSS ≤ 5	RSS = 5	RSS = 6	Pk	Sen	Spe	Pk	Sen	Spe
RMSE1	1.489 ± 0.302	1.342 ± 0.228	1.271 ± 0.235	0.725	60.1	74.6	0.599	53.3	62.4
RMSE2	2.052 ± 0.253	2.029 ± 0.164	1.957 ± 0.192	0.668	64.9	62.1	0.617	58.9	58.6
RMSE3	2.185 ± 0.260	2.262 ± 0.123	2.207 ± 0.175	0.553	33.5	61.6	0.625	66.4	49.7
RMSE10	2.100 ± 0.301	2.265 ± 0.154	2.296 ± 0.174	0.741	50.7	86.4	0.577	47.2	63.7
RMSE11	2.091 ± 0.308	2.261 ± 0.166	2.292 ± 0.178	0.738	51.6	85.5	0.563	45.8	64.7
RMSE16	2.036 ± 0.315	2.213 ± 0.209	2.254 ± 0.202	0.741	56.2	80.1	0.563	46.9	61.4
RMSE17	2.021 ± 0.316	2.202 ± 0.205	2.249 ± 0.207	0.754	58.0	80.7	0.572	48.3	62.9
RMSEα	−0.016 ± 0.033	−0.009 ± 0.028	0.002 ± 0.029	0.663	59.9	63.8	0.604	55.0	58.3
max(Req = 2)δ	0.232 ± 0.055	0.211 ± 0.054	0.191 ± 0.040	0.734	60.9	76.2	0.606	49.0	67.0
max(Req = 2)α	0.524 ± 0.027	0.537 ± 0.022	0.542 ± 0.024	0.686	62.7	62.6	0.567	58.4	52.6
FD(Req = 0.5)TB	0.862 ± 0.054	0.850 ± 0.039	0.832 ± 0.036	0.700	63.1	71.7	0.642	57.6	64.0

Indexes in both trials presents p-value < 0.001; Pk: prediction probability; Sen: sensitivity (%); Spe: specificity (%); FD: first derivative; max: first local maximum. Values are expressed as mean ± std.

, which shows the mean and standard deviation (mean ± std), Pk, Sen and Spe values for both trials, Trials 1 and 2. Only one of the indexes obtained from the slope of the RMSE course at long time scales fulfill the conditions of Table 3, and this was RMSEα. This slope showed higher absolute values in the responsive group than in the unresponsive group, indicating that RMSE in the responsive group decrease faster than in the unresponsive group as function of the long time scale. In Trial 1, the highest Pk (0.754) was obtained at long time scales (ts = 17), while in Trial 2 it was obtained at medium time scales (Pk = 0.625 at ts = 3).

Indexes computed from the AMIF analysis with p-value < 0.001 and Pk > 0.60 are also included in Table 3. In both Trials, it is observed that the values of the indexes FD(Re_{q = 0.5})_TB and max(Re_{q = 2})_δ were higher in the responsive group compared with unresponsive group, where max(Re_{q = 2})_α shows a contrary behavior. For these indexes, the best Pk (0.734) in Trial 1 was obtained with max(Re_{q = 2})_δ, while FD(Re_{q = 0.5})_TB gave the best results (Pk = 0.642) in Trial 2.

The Pearson correlation analysis between indexes provided that the spectral index P_β was highly correlated with RMSE₁, RMSE₂ and FD(Re_{q = 0.5})_TB, with the following correlations (p-value < 0.0005): ρ(P_β, RMSE₁) = 0.797; ρ(P_β, RMSE₂) = 0.869; ρ(P_β, FD(Re_{q = 0.5})_TB) = 0.959. The non-linear index FD(Re_{q = 0.5})_TB was also correlated (p-value < 0.0005) with RMSE₁ and RMSE₂: ρ(FD(Re_{q = 0.5})_TB, RMSE₁) = 0.812; ρ(FD(Re_{q = 0.5})_TB, RMSE₂) = 0.842.

3.3. Multivariate Statistical Analysis

The results that were obtained using a multivariate discriminant analysis are presented in

Table 4. Multivariate statistical analysis in Trial 1 and in Trial 2.

Index	Trial 1			Trial 2
	(2 ≤ RSS ≤ 5 vs. RSS = 6)			(RSS = 5 vs. RSS = 6)
	Pk	Sen (%)	Spe (%)	Pk	Sen (%)	Spe (%)
f1	0.713	64.8	68.1	0.718	63.3	61.6
f2	0.732	65.7	70.2	0.722	60.3	65.3
f3	0.747	62.4	74.9	0.660	59.5	64.1
f4	0.802	68.4	80.1	0.683	58.4	69.5
f5	0.776	66.2	76.7	0.701	60.3	67.3

. Only functions with Pk > 0.60, Sen > 60% and Spe > 60%, simultaneously in at least one trial, are included in this table: f₁ = f(RMSE₂, RMSE_α); f₂ = f(RMSE₂, RMSE_α, CeRemi); f₃ = f(RMSE₁, max(Re₂)_δ, mHR); f₄ = f(RMSE₁, RMSE_α, max(Re₂)_δ); and f₅ = f(RMSE₂, RMSE_α, max(Re₂)_δ). The correlation between the indexes of the linear functions was weak, being the highest 0.381 (p-value < 0.0005) which states between RMSE_α, and max(Re₂)_δ.

The function f₄, which includes RMSE and AMIF indexes, had the best performances (Pk = 0.802) in discriminating responsive from unresponsive group when 2 ≤ RSS ≤ 5 vs. RSS = 6 were considered in Trial 1, but it showed a relative low Pk (0.683) in relation with other functions when RSS = 5 vs. RSS = 6 were compared in Trial 2. On the other hand, the multivariable function f₂, which combines RMSE indexes with CeRemi, yields the highest Pk (0.722) in Trial 2.

3.4. Validation in GAG Reflex during Endoscopy Tube Insertion

To validate the robustness of the proposed univariate indexes and functions f_i for Trials 1 and 2, the same indexes and functions were also considered to study GAG reflex during endoscopy tube insertion in the same database patients.

Table 5. Presence and absence of GAG reflex: univariate analysis.

Index	Groups		–
	GAG = 1 (mean ± std)	GAG = 0 (mean ± std)	Pk	Sen (%)	Spe (%)
CeProp	2.307 ± 0.780	2.415 ± 0.718	0.549	54.4	52.3
CeRemi	1.147 ± 0.811 †	1.396 ± 0.646	0.624	58.4	54.7
BIS	80.5 ± 11.9 †	68.3 ± 15.0	0.738	72.1	63.1
mHR	79.0 ± 15.4 †	72.2 ± 13.4	0.637	57.3	59.9
sdHR	3.16 ± 4.56 †	2.31 ± 3.24	0.615	36.6	76.6
Pβ	0.321 ± 0.199 †	0.229 ± 0.148	0.642	52.5	70.3
Pα	0.215 ± 0.138 †	0.369 ± 0.161	0.766	74.3	67.3
Pθ	0.138 ± 0.063 †	0.170 ± 0.062	0.656	68.3	57.7
Pδ	0.350 ± 0.211 †	0.262 ± 0.162	0.615	50.5	71.5
RMSE1	1.611 ± 0.280 †	1.407 ± 0.297	0.702	62.4	66.4
RMSEα	−0.026 ± 0.032 †	−0.005 ± 0.031	0.687	61.6	66.2
max(Re2)δ	0.254 ± 0.061 †	0.217 ± 0.054	0.690	54.5	70.0

^†: p-value < 0.001. P_k: prediction probability; Sen: sensitivity; Spe: specificity.

shows the mean ± std, Pk, Sen (%) and Spe (%) of the univariate indexes. In order to simplify the information, this table only contains the RMSE and AMIF indexes that were included in the multivariable functions that are shown in Table 4. It can be observed that all statistically analyzed indexes showed significant differences between groups GAG = 1 vs. GAG = 0 with the exception of CeProp.

Comparing with the analysis in Trials 1 and 2, the validation of the univariate indexes showed a similar behavior in the study of the GAG reflex. Indeed, the mean value of indexes BIS, mHR, sdHR, P_β and P_δ was higher (p-value < 0.001) in the responsive groups (GAG = 1) than in the unresponsive group (GAG = 0). Indexes such as CeRemi, P_α and P_θ had a contrary behavior, showing lower mean value (p-value < 0.001) in the responsive groups than in the unresponsive group. RMSE₁ (SampEn at ts = 1) and RMSE_α (slope of RMSE course between 5 ≤ ts ≤ 8) showed higher absolute values in the responsive group than in the unresponsive group. The index max(Re_{q = 2})_δ was higher in the responsive group than in the unresponsive group. The best Pk (0.766) and Sen (74.3%) values were reached with the index P_α, while the index sdHR showed the highest Spe value (76.6%).

The results of the Pearson correlation analysis between indexes indicated that the spectral index P_β was highly correlated with RMSE₁, RMSE₂ and max(Re_{q = 2}), with the following correlations (p-value < 0.0005): ρ(P_β, RMSE₁) = 0.782; ρ(P_β, RMSE₂) = 0.872; ρ(P_β, max(Re_{q = 2})) = 0.948. The non-linear index FD(Re_{q = 0.5})_TB was also correlated (p-value < 0.0005) with RMSE₁ and RMSE₂: ρ(FD(Re_{q = 0.5})_TB, RMSE₁) = 0.861; ρ(FD(Re_{q = 0.5})_TB, RMSE₂) = 0.842.

Finally, the multivariate functions with the best performance in the GAG reflex study were: f₃ = f(RMSE₁, max(Re₂)_δ , mHR) and f₄ = f(RMSE₁, RMSE_α, max(Re₂)_δ). The Pk, Sen and Spe values that were obtained with the function f₃ were 0.723, 64.4% and 70.6%, respectively, and with function f₄ were 0.827, 68.3% and 76.9%, respectively.

4. Discussion

The prediction of response to noxious stimulation is an issue that has not been completely solved, particularly in sedation-analgesia procedures where the patients may be under a sedation state in which the patient cannot be aroused but respond purposefully following repeated or painful stimulation, as the tube insertion during endoscopies, a minimally invasive procedure. Therefore, physiological conditions in these sedation levels are very different from the conditions detectable under general anesthesia where the responses to noxious stimulation are absent. This study demonstrates that a multivariate function of complexity measures as RMSE and AMIF can predict the responsiveness to noxious stimulation, which can be used by anesthesiologist in order to guarantee a comfortable state for the patient during the medical procedure.

Specifically, the function f₄ = f(RMSE₁, RMSE_α, max(Re₂)_δ) had the best performances (Pk > 0.8) in discriminating responsive from unresponsive groups when considering 2 ≤ RSS ≤ 5 vs. RSS = 6 and GAG = 1 vs. GAG = 0, while the function f₂ = f(RMSE₂, RMSE_α, CeRemi) yielded the highest Pk (0.722) when comparing RSS = 5 vs. RSS = 6, being these values higher than the values computed with the BIS index. It is well known that BIS is able to describe hypnotic effect as it was confirmed by results of Trial 1 (Table 3). However, as it was seen in Trial 2 where (2 ≤ RSS ≤ 5) and (RSS = 6) were compared or comparing (GAG = 1) vs. (GAG = 0), the functions based on RMSE and AMIF indexes showed a better capability than BIS to describe the analgesic effect and then to predict the response to noxious stimulation.

Analyzing the results in Table 2 and Table 5, the statistically significant differences of BIS values between responsive and unresponsive groups indicate that higher sedation levels can be associated with low probability of nociception. Since in this study the patients were not under general anesthesia but only under sedation, the expected value of the BIS fluctuates between 100 and 60. Low BIS values are associated with higher sedation levels, which are ideally reached by an adequate combination of hypnotic and analgesic concentrations, generating less probability of nociception. Thereby, Brocas et al. [35] found BIS values significantly higher during a control period than an alfentanil period when evaluated the effect of an intravenous bolus of alfentanil on the variations in BIS level. An increase in BIS values was also observed after a control period tracheal suction. However, BIS might have the same value for different concentrations of drugs and it is possible that in case of low doses of analgesia, a response to noxious stimuli might be observed even at low BIS values.

Several studies [36,37,38] relate nociception with the parameters derived from the analysis of heart rate variability (HRV). It has been demonstrated that the tone of the autonomic nervous system is strongly influenced by anesthetic drugs. Therefore, the traditional parameters of HRV in time domain and frequency domain [39] are influenced by changes in the depth of sedation [40] and can act as indicators of inadequate analgesia [35,41]. The studies [41,42] showed that anesthesia induction decreased the heart rate and that during nociception, the HRV not change when the patient receives adequate analgesia. In this study, the lower values of meanHR in unresponsive state (RSS = 6) confirm that high sedation level can be associated with low heart rate. Furthermore, the lower values of sdHR in unresponsive state denote a lower variability of the heart rate when the analgesia is adequate. An increase of meanHR and sdHR from presence to absence of GAG reflex is most likely to be due to a stimulation of central noradrenergic neurons that could realize a kind of cortical awakens, which reflects adrenergic hyperactivation.

Concerning about EEG bands, changes in δ, θ and α activity have been associated with changes in the sedation levels [43,44,45,46]. Indeed, an increase of the hypnotic and analgesic doses in patients moves EEG activity from a low-amplitude and high-frequency signal to a high-amplitude and low-frequency signal. Concretely, there is a reduction of the β activity and an increment of the α and δ activities when the levels of propofol anesthesia are incremented [46]. In the present study, this behavior was corroborated by the statistical differences observed in the spectral power in α, θ and β frequency bands (see Table 2 and Table 5) between the responsive (2 ≤ RSS ≤ 5, RSS = 5 and GAG = 1) and unresponsive groups (RSS = 6, GAG = 0). The spectral power in δ band showed different behavior in 2 ≤ RSS ≤ 5 vs. RSS = 6 and GAG = 1 vs. GAG = 0 that between RSS = 5 vs. RSS = 6. Specifically, the increase of the δ activities with high level of sedation was only observed in the comparison between RSS = 5 vs. RSS = 6, while the δ power decreased when comparing 2 ≤ RSS ≤ 5 vs. RSS = 6 and GAG = 1 vs. GAG = 0. This behavior might be explained by assuming that the strong response to nail bed compression (RSS < 5) and the gag reflex after tube insertion (GAG = 1) are associated with lower sedation levels than the sluggish response to nail bed compression (RSS = 5), hence they might contain also ocular activity. In this way, slow eye movements can affect the spectral power in δ band by causing its increase also in low sedation levels RSS < 5 and GAG = 1. It is also important to remember that patients in this study were under sedation-analgesia procedure and not under general anesthesia, and therefore it was not observed a significant increment in δ power at high sedation levels.

The results of RMSE (Table 3 and Table 5) indicate that the complexity of the EEG signal, in responsive and unresponsive groups, was higher in the low frequency bands (ts > 5), corresponding to δ, θ and α bands, than in high frequency bands (ts < 3). These changes were also observed in the index RMSE_α which is the slope of RMSE course corresponding to time scale factors ranging from ts = 5 to ts = 10, corresponding approximately to α band at a sampling frequency of 128 Hz. Besides that, complexity of the EEG signal was higher in responsive than unresponsive group at short time scales (ts < 3). This behavior can be interpreted by considering that short time scales contain frequency bands where scalp and facial muscle artifact are presented, suggesting that this muscle activity, which is more important in patients with low sedation level, is the responsible for the higher complexity in the responsive state at those time scales. In this sense, it is possible to associate an increased activity in the facial muscles with a greater probability of pain.

At long time scales, in which muscle activity has been eliminated, the entropy values indicate that the EEG signal contains more regular patterns in the responsive state than in the unresponsive state. This regularity increases as the time scale is larger for the responsive state but remains almost constant for the unresponsive state. Therefore, the results demonstrated that patients in unresponsive group show a more complex EEG activity in low frequency bands than patients in responsive group. One explication of this result is related with the fact that EEG activity becomes slower as the sedation level increases, and thus, it is expected that patients in unresponsive group present a slower EEG signal than responsive group, increasing the signal complexity in low frequency bands. A similar situation was reported in a study during sevoflurane anesthesia where MSE curves showed an ascending or flat trend for deep anesthesia and a descending trend for waking state [20].

Finally, regarding the AMIF indexes (Table 3 and Table 5), it can be denoted that EEG behavior was more complex in total and in α frequency band for the responsive group compared with the unresponsive group, while the EEG of responsive group was more regular in δ band than unresponsive group. This behavior was in agreement with the results of RMSE at long time scale, which showed lower values, and then more regular EEG, in the responsive group compared to the unresponsive group.

5. Conclusions

This study suggests that a multivariate model based on indexes extracted from RMSE, AMIF, heart rate and CeRemi enhance the predictability of different stimulation responses during sedation-analgesia. Patients in unresponsive state showed a more complex EEG activity in low frequency bands than patients in responsive state. This can be related with the fact that EEG activity becomes slower as the sedation level increases, and thus, it is expected that patients in unresponsive state present a slower EEG signal than in responsive state, thus increasing the signal complexity in low frequency bands. To conclude, the indexes analyzed in this study can be used as guide by the anesthesiologist in reducing stress responses to nociceptive stimuli. The implementation of the proposed tool based on linear, nonlinear and clinical indexes into standards for patient monitoring during sedation might be helpfully exploited for a fine tuning of the levels of anesthesia, thus speeding up patient’s recovery after surgery and avoiding side effects.