# Heart rate variability analysis for the prediction of pre-arousal during propofol-remifentanil general anaesthesia: A feasibility study.

Anne Wojtanowski, Maxence Hureau, Camille Ternynck, Benoit Tavernier, Mathieu Jeanne, Julien de Jonckheere

<figure><img src="image_1.png" /></figure>

## To cite this version:

Anne Wojtanowski, Maxence Hureau, Camille Ternynck, Benoit Tavernier, Mathieu Jeanne, et al.  
Heart rate variability analysis for the prediction of pre-arousal during propofol-remifentanil general anaesthesia: A feasibility study.. PLoS ONE, 2024, PLoS One, 19, pp.e0310627. 10.1371/journal.pone.0310627. hal-04844925

HAL Id: hal-04844925

https://lilloa.hal.science/hal-04844925v1

Submitted on 18 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
RESEARCH ARTICLE

# Heart rate variability analysis for the prediction of pre-arousal during propofol-remifentanil general anaesthesia: A feasibility study

Anne Wojtanowski<sup>1,2*</sup>, Maxence Hureau<sup>3,4</sup>, Camille Ternynck<sup>5</sup>, Benoit Tavernier<sup>2,3</sup>, Mathieu Jeanne<sup>1,3,4</sup>, Julien de Jonckheere<sup>1,2</sup>

1 CIC IT 1403, CHU Lille, Lille, France, 2 ULR 2694—METRICS, Univ. Lille, Lille, France, 3 Anesthésie Réanimation CHU Lille, CHU Lille, Lille, France, 4 ULR 7365 - GRITA - Groupe de Recherche sur les Formes Injectables et les Technologies Associées, Univ. Lille, Lille, France, 5 Département de Biostatistique, CHU Lille, Lille, France

* anne.wojtanowski@chu-lille.fr

<figure><img src="image_3.png" /></figure>

## Abstract

Accidental awareness during general anaesthesia is a major complication. Despite the routine use of continuous electroencephalographic monitoring, accidental awareness during general anaesthesia remains relatively frequent and constitutes a significant additional cost. The prediction of patients' arousal during general anaesthesia could help preventing accidental awareness and some researchers have suggested that heart rate variability (HRV) analysis contains valuable information about the patient arousal during general anaesthesia. We conducted pilot study to investigate HRV ability to detect patient arousal. RR series and the Bispectral Index™ (BIS™) were recorded during general anaesthesia. The pre-arousal period T0 was defined as the time at which the BIS™ exceeded 60 at the end of surgery. HRV parameters were computed over several time periods before and after T0 and classified as "BIS™<60" or "BIS™≥60". A multivariate logistic regression model and a classification and regression tree algorithm were used to evaluate the HRV variables' ability to detect "BIS™≥60". All the models gave high specificity but poor sensitivity. Excluding T0 from the classification increased the sensitivity for all the models and gave AUCROC>0.7. In conclusion, we found that HRV analysis provided encouraging results to predict arousal at the end of general anaesthesia.

<figure><img src="image_4.png" /></figure>

### OPEN ACCESS

**Citation:** Wojtanowski A, Hureau M, Ternynck C, Tavernier B, Jeanne M, de Jonckheere J (2024) Heart rate variability analysis for the prediction of pre-arousal during propofol-remifentanil general anaesthesia: A feasibility study. PLoS ONE 19(10): e0310627. https://doi.org/10.1371/journal.pone.0310627

**Editor:** Ernesto Iadanza, University of Siena; Universita degli Studi di Siena, ITALY

**Received:** March 6, 2024

**Accepted:** September 4, 2024

**Published:** October 31, 2024

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pone.0310627

**Copyright:** © 2024 Wojtanowski et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the manuscript and its Supporting Information files.

## Introduction

Accidental awareness during general anaesthesia (AAGA) is a major complication of anaesthesia and can cause long-term psychiatric problems. The incidence of AAGA ranges from 1 per 15,000 general anaesthesia (GA) procedures (when reported by anaesthesiologist) to 1–2 per 1000 (when reported by patients) [1–3]. The risk factors for AAGA include the type of surgery, demographic characteristics, the specific phase of anaesthesia (induction and emergence) and,
**Funding:** The author(s) received no specific funding for this work.

**Competing interests:** I have read the journal's policy and the authors of this manuscript have the following competing interests: Julien De Jonckheere and Mathieu Jeanne are shareholders and scientific consultants for MDoloris Medical Systems. The remaining authors have disclosed that they do not have any potential conflicts of interest.

most clearly, the use of neuromuscular blocking agents [4]. Several strategies for avoiding AAGA have been suggested. Clinical monitoring can be used to detect the occurrence of light anaesthesia, although neuromuscular blockade clearly interferes with this approach [5]. Another option is electroencephalographic (EEG) monitoring, in order to track the patient's level of consciousness via the Bispectral Index™ (BIS™, Medtronic, Minneapolis, MN, USA), the SedLine® (Masimo, Irvine, CA, USA) or entropy monitoring (GE Healthcare, Chicago, IL, USA). EEG monitors are considered as the gold standard to measure the level of consciousness and to detect AAGA. However, the increasing use of neuromuscular blocking agents in modern anaesthesia significantly impacts these measures. A study of myorelaxed-only volunteers showed that the use of neuromuscular blocking agents led to a decrease in EEG indexes like the BIS™ [6]. Similar results have been published with other EEG monitors [7]. Simplified EEG signal used for depth of hypnosis monitoring during general anaesthesia rely on a limited number of frontal electrodes, which capture the electrical potential of cortical and subcortical layers [8]. Propofol mediated unconsciousness has been related to an alpha/delta pattern and slow oscillations in the EEG [9, 10], but some authors found evidence that connected consciousness sometimes occurs in spite of these patterns being present [11]. Sensory disconnection and loss of consciousness may have different EEG spectral markers, with anterior and posterior cingulate regions probably playing a major role in consciousness [12]. Although a BIS™ below 60 is typically targeted during GA, one cannot be absolutely sure that awareness (i.e. AAGA) is precluded in this context. Although several studies have shown that the use of the BIS™ significantly reduced the occurrence of AAGA [13], recent evidence shows that awareness remains possible even when the EEG indexes are on target during anaesthesia [14].

Blood pressure monitoring and electrocardiography (ECG) are mandatory during GA. Acquisition of the ECG signal enables the evaluation of heart rate variability (HRV), which is a valid, well-established, non-invasive index of autonomous nervous system (ANS) responses during GA. For example, the Analgesia Nociception Index (ANI) is an HRV-based index of the nociception/antinociception (NAN) balance during GA [15, 16]. This technology is now used routinely to personalize analgesia: ANI monitoring has been shown to reduce or dynamically modify the use of opioids during and after surgery [17–19].

Many ANS monitors are used in clinical practice to assess the NAN balance and some of these ANS monitor are based, as a part, on HRV analysis [20]. On the other hand, hypnotic drugs may also affect HRV notably through its sympatholytic effect [21, 22]. In their study, R Huhle et al. demonstrated that propofol induction significantly reduced HRV features [23]. Base on this statement, Zhan J. et al. developed a HRV-derived system to evaluate the level of consciousness during GA and demonstrated good agreement between HRV and depth of anaesthesia [24]. We hypothesized that HRV is correlated with EEG signals of cortical awakening (i.e. when the BIS™ exceeds 60 but the patient is not yet conscious) at the end of GA. Hence, the objectives of the present ancillary, retrospective analysis of data collected during a prospective, single-centre, randomized clinical trial were to measure various HRV variables during propofol-remifentanil anaesthesia and to test the ability of an HRV-based model to predict “BIS™≥60” periods.

# Materials and methods

## Ethics

All methods were carried out in accordance with relevant guidelines and regulations. The study data were generated during the randomized ANI REMI LOOP clinical study (“Automatic control of ANI-guided remifentanil administration during propofol general anaesthesia”). This study aimed to demonstrate the efficacy and safety of an original closes-loop
remifentanil administration medical device (ANI-loop) in comparison to the standard prac-
tice. This study demonstrated that the use of the ANI-loop device significantly reduced the
total amount of remifentanil administered. The study was registered at clinicaltrials.gov
(NCT03556696) and with the French National Agency for Medicines and Health Product
Safety (IDRCB 2017-A00858-45) and was conducted at Lille University Hospital between 27th
June 2018 and 27th June 2022 for which patients gave a written informed consent. The present
study is a retrospective analysis of the ANI REMI LOOP data. Under French regulations, stud-
ies concerning data reuse do not require additional authorisation by an ethical committee and
only require informing participants. All patients were contacted and none of them objected to
the use of their data for the study. The data were accessed after the completion of the ANI
REMI LOOP study on September 2022. Authors had access to information that could identify
participants during data collection but not during this data analysis.

Procedures

The participants in the ANI REMI LOOP trial had undergone skin graft surgery under GA because of an accidental burn. Only adult patients with a sinus rhythm and no history of dysautonomia or diabetes mellitus were included in the ANI REMI LOOP study. Patients who presented with a non-sinus rhythm, had a pacemaker, heart transplant, baseline haemodynamic measures considered as hypotension or bradycardia or the opposite, heart rate baseline >120min⁻¹ or systolic blood pressure baseline >160 mmHg, were excluded. All patients were monitored throughout the procedure with standard techniques, including ECG, plethysmography, non-invasive blood pressure monitoring, and spirometry. Furthermore, all patients were monitored with the PhysioDoloris® device (V1.4.1.0, MDoloris, Loos, France) and the BIS™ VISTA system (Medtronic, Minneapolis, USA). The standardized anaesthetic protocol comprised midazolam (0.08 mg.kg⁻¹), remifentanil and propofol for the induction and maintenance of GA using a target-controlled infusion device (Orchestra® Base, Fresenius Kabi AG, Bad Homburg, Germany). The remifentanil and propofol targets were initially 6.0 ng.ml⁻¹ and 5.0 µg.ml⁻¹, respectively. Propofol was adjusted during GA so that the BIS™ stayed in the range [40–60] and remifentanil was adjusted so that the ANI stayed in the range [50–70]. Once myorelaxation had been achieved (cisatracurium, 0.15 mg.kg⁻¹), the trachea was intubated, and the propofol and remifentanil targets were lowered to 3 µg.ml⁻¹ and 3 ng.ml⁻¹, respectively. All patients were ventilated using assisted controlled ventilation with tidal volume set at 8 ml.kg⁻¹ of ideal body weight, and respiratory rate was set of 12 min⁻¹ (GE Aisys anaesthesia machine). The fraction of end-tidal (Fet) CO₂ was maintained in the 30–35 mmHg range. Propofol effect-site target was adjusted manually using the Schiller model. The remifentanil effect-site target was adjusted manually using the Minto model in the ANI REMI LOOP trial's control group and automatically in the interventional group. During the procedure, neuromuscular blockade was monitored by a train of four (TOF) stimulation [25]. At the end of the intervention, propofol and remifentanil were stopped only if the TOF ratio was up to 90%. Pre-emptive analgesia consisted of the local administration of ropivacain 2 mg.ml⁻¹ during the procedure and paracetamol 1 g and sufentanil 0.12 µg.kg⁻¹ (maximum: 10 µg) at the end of the procedure. Patients from whom the BIS™ recordings did not include the awakening phase were excluded for this retrospective analysis. Recordings showing synchronization problems between HRV and BIS™ data were not analyzed.

ECG signal processing

The study database included the BIS™ and PhysioDoloris® recordings. Post-hoc analysis of
the RR data exported from the PhysioDoloris® monitor enabled us to compute the HRV
variables. Ectopic beats and RR series artefact were detected and replaced by a linear interpolation using a specific algorithm [26, 27]. After RR series artefact filtering, the various HRV variables computed were:

*   **Time analysis:**

o The root mean square of successive differences (RMSSD) between adjacent RR intervals was computed over 100 RR intervals (approximately 2 minutes.). The RMSSD corresponds to short-term variations in the parasympathetic nervous system's activity [28].

○ The standard deviation of normal to normal (SDNN) RR intervals was computed over 100 RR intervals, which correspond to approximately 2 minutes. The SDNN corresponds to the overall variability of the ANS's activity.

○ Short-term variability (STV) was computed over 1 minute. After 4-Hz resampling of the RR series, the 1-minute window was divided into sixteen 3.75 s epochs. The average RR value was computed for each epoch. STV was defined as the mean difference between successive 3.75 s periods. STV corresponds to short-term variations in the parasympathetic nervous system's activity.

○ Long-term variation (LTV) was computed over 1 min. After 4-Hz resampling of the RR series, the 1-minute window is divided into sixteen 3.75 s epochs (i.e. 16 epochs). The average RR value was computed for each epoch. LTV was defined as the difference between the maximum and minimum values of the 16 epochs. LTV corresponds to the overall variability in the ANS's activity.

Although STV and LTV are mainly used for fetal heart rate variability analysis, they use different methods than SDNN and RMSSD to assess the overall and short term variability of the ANS and could therefore provide additional information in the present clinical context.

*   **Spectral analysis:**

o Low frequency (LF) was based on spectral analysis of the RR series: for an adult, the standard bandwidth was defined as [0.05–0.15 Hz]. LF was computed via a wavelet transform (e.g. Daubechies 4-wavelet) on a 64-second window resampled at 8 Hz using a linear interpolation. LF reflects both sympathetic and parasympathetic modulations and is mainly influenced by baroreflex activity.

○ High frequency (HF) was based on spectral analysis of the RR series: for an adult, the standard bandwidth was defined as [0.15–0.4 Hz]. HF was computed in the same way as LF. HF reflects parasympathetic modulations only and is mainly influenced by respiratory sinus arrhythmia.

*   **Graphical analysis:**

o The ANI is a graphical index of the magnitude of HF oscillations. The RR series were resampled at 8 Hz, centred, and normalized in 64-second windows. The signal was band-pass filtered a wavelet transform in order to retain HF content only. The magnitude of these oscillations was measured as the surface area between the upper and lower envelopes. The ANI gives the relative magnitude of the HF spectral content on a 0-to-100 scale. The ANI is therefore a measure of parasympathetic activity and is mainly influenced by respiratory sinus arrhythmia.

Each HRV measure was computed over a 60-second moving window with a 1-second sliding period and then averaged over 3 minutes.
Using the BIS signal as the existing gold standard, we defined T0 as the time at which the BIS™ exceeded 60 at the end of the surgery, i.e. after the remifentanil and propofol targets had been set to zero in order to awake the patient. We collected HRV data 5, 10, 20 and 30 minutes before T0 and 2 and 5 minutes after T0 defined respectively as T-5, T-10, T-20, T-30, T0, T+2 and T+5. The periods between these time points were rated as “BIS™<60” or “BIS™≥60”.

## Statistical analysis

The normality of distribution was assessed using histograms and the Shapiro-Wilk test. Quantitative variables were expressed as the mean ± standard deviation when normally distributed and as median [interquartile range (IQR)] when not. Qualitative variables were expressed as the number (percentage).

In order to identify highly correlated variables, a principal component analysis (PCA) followed by a varimax rotation was applied on all the available variables. Pearson's correlation coefficients between quantitative variables were described. Based on these results and clinical/physiological expertise, a set of uncorrelated variables was also selected. We first assessed the associations between all uncorrelated candidates predictors and the cortical awakening using univariable logistic regressions and secondly, we used multivariable method to develop a predictive model of cortical awakening. For each continuous predictor, the log-linearity assumption was assessed using restricted cubic spline functions. The absence of collinearity between variables was checked by calculating the variance inflation factors. When the log-linearity was rejected, variables were log-transformed. To account for the number of candidate predictors and limit the risk of over-optimism, this model was built using bootstrap resampling (n = 500) with, in each sample, automated stepwise backward selection procedure (with a removal criteria of p>0.05) from among all the uncorrelated candidate predictors [29]. A variable was kept in the final model if it was selected in at least 70% of these 500 analyses. The final predictive model was obtained by logistic regression of the retained variables. We computed the odds ratio (OR) and its 95% confidence intervals (CIs). The model's predictive ability was assessed as the area under the receiver operating characteristic curve (AUROC) [95%CI]. In parallel, we also developed a decision tree by using a classification and regression tree (CART) algorithm (Rpart package, R software) [30]. The models' respective levels of performance were assessed by calculating the sensitivity (Se), specificity (Sp), positive predictive value (PPV) and negative predictive value (NPV). The analyses were conducted on the complete set of cases. All tests were two-tailed, and the threshold for statistical significance was set to p<0.05. The data were analysed using the SAS software (version 9.4, SAS Institute Inc., Cary, NC, USA), and R software (version 3.6.1) [30] and Sipina research software (version 3.13).

## Results

A total of 52 patients were included in the initial clinical trial. Eleven of the 52 datasets lacked reliable synchronisation between the BIS™ and HRV signals, and 7 BIS™ recordings did not include the awakening phase (the electrode was removed before the end of the procedure). Hence, 34 datasets were available for the analyses and generated 208 measurements: 70 with “BIS™≥60” and 138 with “BIS™<60” (for details, see Table 1). Table 2 describes the characteristics of the study participants, and Table 3 describes the characteristics of the HRV indices.
<table><caption>Table 1. Time points and measurements.</caption><thead><tr><td>Time point</td><td>Number of available measurements</td><td>Number of measurement with a BIS<sup>TM</sup> ≥ 60</td></tr></thead><tbody><tr><td>T-30</td><td>34</td><td>1</td></tr><tr><td>T-20</td><td>34</td><td>0</td></tr><tr><td>T-10</td><td>34</td><td>0</td></tr><tr><td>T-5</td><td>34</td><td>1</td></tr><tr><td>T0</td><td>34</td><td>34</td></tr><tr><td>T+2</td><td>23</td><td>21</td></tr><tr><td>T+5</td><td>15</td><td>13</td></tr></tbody></table>

Time points definition: T-30: 30 minutes before BIS<sup>TM</sup> ≥ 60; T-20: 20 minutes before BIS<sup>TM</sup> ≥ 60; T-10: 10 minutes before BIS<sup>TM</sup> ≥ 60; T-5: 5 minutes before BIS<sup>TM</sup> ≥ 60; T0: the moment that the BISTM ≥ 60; T+2: 2 minutes after BIS<sup>TM</sup> ≥ 60; T+5: 5 minutes after BIS<sup>TM</sup> ≥ 60

https://doi.org/10.1371/journal.pone.0310627.t001

<table><caption>Table 2. Demographic and clinical characteristics of the study participants.</caption><thead><tr><td>Characteristic</td><td></td></tr></thead><tbody><tr><td>Sex, male</td><td>13/34 (38.23)</td></tr><tr><td>Age, y</td><td>41.50 [34.25–50.00]</td></tr><tr><td>Weight, kg</td><td>75.59 ± 13.49</td></tr><tr><td>Height, cm</td><td>171.6 ± 8.16</td></tr><tr><td>ASA score</td><td></td></tr><tr><td>ASA I</td><td>14/34 (41.18)</td></tr><tr><td>ASA II</td><td>19/34 (55.88)</td></tr><tr><td>ASA III</td><td>1/34 (2.94)</td></tr><tr><td>Duration of surgery, min</td><td>48.09 ± 14.89</td></tr><tr><td>Duration of anaesthesia, min</td><td>91.76 ± 15.80</td></tr><tr><td>Time between T0 and extubation, min</td><td>5.43 [3.277–8.932]</td></tr><tr><td>Anaesthetic management</td><td></td></tr><tr><td>propofol, mg</td><td>682.5 ± 239.90</td></tr><tr><td>remifentanil, µg</td><td>757.0 [508.5–1006.8]</td></tr><tr><td>atropine</td><td>2/34 (5.88)</td></tr><tr><td>ephedrine</td><td>7/34 (20.58)</td></tr><tr><td>noradrenalin</td><td>3/34 (8.88)</td></tr><tr><td>Analgesia management</td><td></td></tr><tr><td>sufentanil, µg</td><td>7.956 ± 1.13</td></tr><tr><td>ANI REMI LOOP group</td><td>18/34 (52.94)</td></tr><tr><td>Gabapentin</td><td>6/33 (18.18)</td></tr><tr><td>Tramadol</td><td>12/32 (37.5)</td></tr><tr><td>Paracetamol-opium</td><td>4/31 (12.90)</td></tr><tr><td>Short-acting oral morphine</td><td>3/30 (10)</td></tr><tr><td>Long-acting oral morphine</td><td>3/30 (10)</td></tr><tr><td>Amitriptyline</td><td>1/30 (3.33)</td></tr><tr><td>Ketoprofen</td><td>0/32 (0)</td></tr></tbody></table>

Values are expressed as the number/total number (%), the mean ± standard deviation, or the median [IQR].

ASA: American Society of Anaesthesiologists Physical Status Classification System

https://doi.org/10.1371/journal.pone.0310627.t002
<table><caption>Table 3. HRV variables, as a function of BIS<sup>TM</sup> < 60 and BIS<sup>TM</sup> ≥ 60 periods.</caption><thead><tr><th>Variable</th><th>BIS<sup>TM</sup> < 60 (n = 138)</th><th>BIS<sup>TM</sup> ≥ 60 (n = 70)</th></tr></thead><tbody><tr><td>BIS</td><td>41.6 [34.4–47.2]</td><td>64.7 [61.7–74.2]</td></tr><tr><td>HR</td><td>64.79 ± 13.1</td><td>65.82 ± 13.86</td></tr><tr><td>ANI</td><td>68.78 ± 17.12</td><td>61.65 ± 16.13</td></tr><tr><td>STV</td><td>6.81 [4.55–9.85]</td><td>9.03 [5.81–13.50]</td></tr><tr><td>SDNN</td><td>24.27 [16.07–32.99]</td><td>30.89 [20.65–44.49]</td></tr><tr><td>LTV</td><td>64.52 [44.98–89.47]</td><td>90.69 [60.46–127.14]</td></tr><tr><td>RMSSD</td><td>16.94 [8.69–26.74]</td><td>21.35 [13.61–27.69]</td></tr><tr><td>LF</td><td>0.17 [0.05–0.35]</td><td>0.29 [0.13–0.59]</td></tr><tr><td>HF</td><td>0.03 [0.01–0.07]</td><td>0.06 [0.02–0.13]</td></tr></tbody></table>

Values are expressed as the number/total number (%), the mean ± standard deviation, or the median [IQR].

BIS: Bispectral Index<sup>TM</sup>, HR: heart rate, ANI: Analgesia Nociception Index, STV: short-term variation, SDNN: standard deviation of RR intervals, LTV: long term variation, RMSSD: root mean square of successive difference in RR intervals, LF: low frequency, HF: high frequency

https://doi.org/10.1371/journal.pone.0310627.t003

## Identification of uncorrelated variables

The PCA with varimax rotation demonstrated that SDNN, RMSSD, LTV and STV were highly correlated with the first factor, HF and LF were highly correlated with the second factor, and heart rate and ANI were highly correlated with the third factor (Fig 1 and Table 4).

Three sets of variables were strongly correlated with the three factors. Both the ANI and the heart rate were correlated with factor 3. Pearson's correlations coefficients between the variables were computed (Fig 1) and given that the ANI and the heart rate were weakly correlated with each other (Pearson's coefficient = -0.28), both variables were retained. For variables correlated with factor 1 (the SDNN, RMSSD, LTV and STV), we retained the SDNN and RMSSD. LTV and STV were highly correlated with the SDNN (Pearson's coefficient = 0.93 and 0.9, respectively). We retained the SDNN (rather than LTV and STV) because it is a well-known index of overall ANS activity; in contrast, LTV and STV were initially developed for foetal diagnosis and are only used in this particular clinical context [31, 32]. Although the RMSSD was also correlated with the SDNN (Pearson's coefficient = 0.82), we decided to keep it out of the model because it represents another autonomic pathway (the parasympathetic nervous system). Both LF and HF were correlated with factor 2 but we retained LF only; HF is solely influenced by parasympathetic activity, which is also accounted for by the ANI and the RMSSD. Hence, five variables were retained (given in bold in Table 4).

## Construction of cortical awakening prediction models

Shape of associations between each parameter and awakening were shown in Fig 2. There is no major deviation in log-linearity assumptions for all parameters, except for SDNN and LF where log-transformations were needed.

In univariate logistic regressions, the heart rate was the only variable that did not predict cortical awakening (Fig 3). The SDNN, RMSSD, LF and ANI were significant for the prediction of BIS<sup>TM</sup> ≥ 60, with AUROC of 0.65, 0.57, 0.64 and 0.64, respectively.

A multivariate logistic regression model with bootstrap resampling and stepwise backward selection (five variables) led to selection of the SDNN and the ANI. The AUROC of the model with these two variables was 0.74, which corresponds to a good level of discriminant power. By considering a threshold of 0.39, we obtained a sensitivity of 0.53, a specificity of 0.85, a PPV of 0.64, and an NPV of 0.78. However, when considering the 33 (out of 70) periods misclassified as
<figure><img src="image_5.png" /><figcaption>**Fig 1. Pearson coefficients for the correlations between the eight variables.** Variables were considered to highly correlated if Pearson's coefficient was > 0.70.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0310627.g001

"BISTM≥60", 18 were related to T0 (Table 5). Excluding T0 from the classification increased the sensitivity to 0.58. For T2, 7 of the 21 "BISTM≥60" were misclassified, giving a specificity of 0.67.

Next, a CART decision tree algorithm was applied (Fig 4). The five uncorrelated variables were used in the algorithm. Using this algorithm, the variable that first classified the recordings as "BISTM≥60" or "BISTM<60" was the SDNN, followed by the ANI, LF, and the SDNN again (Fig 4). This model correctly classified 38 of the 70 "BISTM≥60" periods and 124 of the 138 "BISTM<60" periods, giving a sensitivity of 0.54, a specificity of 0.90, a PPV of 0.73, and an NPV of 0.79. However, when considering the 32 (out of 70) "BISTM≥60" periods misclassified as "BISTM<60", 20 were related to T0 (Table 6). Excluding T0 from the classification led to a sensitivity of 0.68. For T2, 16 of the 21 "BISTM≥60" periods were correctly classified, giving a sensitivity of 0.76.
<table><caption>Table 4. Rotated factor loadings.</caption><thead><tr><th></th><th>Factor 1</th><th>Factor 2</th><th>Factor 3</th></tr></thead><tbody><tr><td>SDNN</td><td><strong>0.94</strong></td><td>0.26</td><td>0.10</td></tr><tr><td>RMSSD</td><td><strong>0.81</strong></td><td>0.17</td><td>0.28</td></tr><tr><td>LTV</td><td><strong>0.61</strong></td><td>0.03</td><td>-0.26</td></tr><tr><td>STV</td><td><strong>0.93</strong></td><td>0.20</td><td>0.00</td></tr><tr><td>LF</td><td>0.21</td><td><strong>0.97</strong></td><td>0.03</td></tr><tr><td>HF</td><td>0.19</td><td><strong>0.97</strong></td><td>-0.05</td></tr><tr><td>HR</td><td>-0.15</td><td>0.01</td><td><strong>-0.69</strong></td></tr><tr><td>ANI</td><td>-0.12</td><td>0.00</td><td><strong>0.84</strong></td></tr></tbody></table>

Absolute values greater than 0.6 are given in bold. Variables retained in the multivariable models are also given in bold. SDNN: standard deviation of RR intervals, RMSSD: root mean square of successive differences in RR intervals, LF: low frequency HR: heart rate, ANI: Analgesia Nociception Index.

https://doi.org/10.1371/journal.pone.0310627.t004

Discussion

This exploratory ancillary study of a number of HRV variables during awakening from propofol/remifentanil- GA was based on data recorded for 34 patients in the periods before and after the BIS<sup>TM</sup> crossed the threshold of 60; this usually occurs several minutes before the patient awakes. Our study's objective was to test the ability of an HRV-based model to predict "BIS<sup>TM</sup>≥60" periods. We developed two multivariate models: the first using multivariate logistic regression, and the second using a CART algorithm.

HRV analysis results interpretation

The regression analysis enabled us to evaluate each HRV variable's ability to predict "BIS<sup>TM</sup>≥60" periods (Fig 3). Unsurprisingly, the heart rate did not predict the early stages of awakening (OR [95%CI]: 1.01 [0.98-1.03]); this was probably because many factors influence the sinus node during this phase of GA. In contrast, the SDNN, RMSSD, LF and ANI highlighted significant differences between "BIS<sup>TM</sup>≤60" and "BIS<sup>TM</sup>≥60" periods but were poor predictors (AUROC of 0.65, 0.57, 0.64 and 0.64 for SDNN, RMSSD, LF, and ANI, respectively). The poor performances of ANI and RMSSD are in line with our current understanding of parasympathetic monitoring, which is considered, in the particular context of surgery under GA, to mostly reflect the NAN balance. The results are also in agreement with Boselli et al.'s finding that the ANI was a good predictor of the hemodynamic response to surgical stimulation but a poor predictor of the level of sedation [16]. In their study, Aragón-Benedí, C et al. demonstrated an association between ANI and the Richmond Agitation-Sedation Scale (RASS) in sedated intensive care unit patients [33]. However, although significant, the correlation was weak (r = 0.246, p = 0.009).

The results of the PCA and Pearson's correlation test showed that less specific HRV measures (such as STV and LTV) were strongly correlated with the SDNN (Pearson coefficients of 0.90 and 0.93, respectively (Fig 1)); hence we selected the SDNN alone for the construction of a predictive model. From a clinical point of view, the SDNN is a well-known HRV variable but is usually computed over several hours; in contrast, we used very short time windows (approx. 2 minutes for computation and 3 minutes for averaging). The SDNN during GA has been linked to the general decrease in ANS variability under the influence of (in most cases) hypnotic drugs [34]. It is therefore not surprising that the SDNN increases once propofol administration has been discontinued, leading to cortical awakening in the absence of other hypnotic drugs (as was the case in our series). From a physiopathological standpoint, the SDNN
<figure><img src="image_6.png" /><figcaption>**Fig 2. Relationship between awakening and HRV variables.** P-values of the likelihood ratio test comparing the full model (including both nonparametric component and linear terms to the model including a linear term only.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0310627.g002

measured over 2 minutes is likely to be influenced by both sympathetic and parasympathetic activities. Hence, the SDNN might represent the overall modulation of the ANS, whereas LF and HF reflect the HRV's information content in detail.

The RMSSD was strongly correlated with the SDNN (Fig 1), even though the latter gave high ORs during the awakening period; this observation shows how complex HRV measurements are. Interestingly, the LF spectral content was correlated strongly with the RMSSD (Pearson's correlation coefficient: 0.87 (Fig 1)). Like the SDNN, LF had a significant OR [95% CI] of 1.47 [1.17–1.84] (Fig 3).

Although the univariate analysis did not accurately predict BIS<sup>TM</sup>≥60 periods, the multivariate logistic regression model had a AUROC of 0.74 and good specificity (0.85) but poor sensitivity (0.53). Decision trees generated with CART algorithm also gave high specificity (0.90) but gave poor sensitivities (0.54). Interestingly, the SDNN and the ANI were the most discriminant indexes in the two models. In the early stages of arousal, this discriminant ability might be related to activation of subcortical structures previously inhibited by hypnotic drugs.
<figure><img src="image_7.png" /><figcaption>Fig 3. Results of the univariate analysis for the prediction of BIS<sup>TM</sup>≥60.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0310627.g003

By combining statistically significant, independent measurements and applying the CART algorithm, we produced a decision tree that described the study variables' respective abilities to predict awakening status (as measured by a BIS<sup>TM</sup> below or above 60). SDNN values appear on both the first and fourth branches of the decision tree. The ANI appears in the tree's second line. An important marker of the decision tree's performance is shown in Table 6: most of the misclassified values occurred at T0, when the BIS<sup>TM</sup> just exceeded 60. One must bear in mind that HRV measurements are averaged over three minutes; this necessarily induces a delay and might explain why the variables did not "see" the arousal. This hypothesis is strengthened by the greater sensitivity observed at T2 (vs. T0). The CART model's performance has clinical value because the patient was extubated a median of 5.4 [3.3–8.9] minutes after T0; hence, HRV measurements could alert the anaesthesiologists to the patient's imminent arousal.

To the best of our knowledge, the present feasibility study is one of the first to have combined exploratory statistical methods, the drastic selection of HRV variables, and clinical knowledge about the ANS in the context of GA. The study's main strength was the use of a well-documented dataset from a recent, standardized, prospective clinical trial in which pro-pofol-remifentanil GA was guided by continuous monitoring of the BIS<sup>TM</sup> and ANI. Heart rate, blood pressure, motor blockade, NAN balance, and ventilator settings were recorded continuously. Pre-emptive analgesia was provided by a single bolus of sufentanil at the end of surgery. Lastly, in this ancillary study, we only assessed the end of the anaesthetic procedure in

<table><caption>Table 5. Values misclassified by the regression model including the ANI and the SDNN.</caption><thead><tr><th>Time</th><th>-30</th><th>-20</th><th>-10</th><th>-5</th><th>0</th><th>2</th><th>5</th></tr></thead><tbody><tr><td>BIS<sup>TM</sup>≥60 classified as BIS<sup>TM</sup>&lt;60</td><td>1</td><td>0</td><td>0</td><td>0</td><td>18</td><td>7</td><td>7</td></tr><tr><td>BIS<sup>TM</sup>&lt;60 classified as BIS<sup>TM</sup>≥60</td><td>7</td><td>7</td><td>2</td><td>5</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0310627.t005
<figure><img src="image_8.png" /><figcaption>**Fig 4.** The decision tree generated by the CART algorithm. For each leaf, the predicted class is given in bold, the number of individuals in the BIS<sup>TM</sup>≥60 class is given on the left, and the number of individuals in the BIS<sup>TM</sup><60 class is given on the right.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0310627.g004

order to ensure the absence of neuromuscular blockage, which can negatively affect BIS<sup>TM</sup> measurements. Therefore, despite the known pitfalls of BIS monitoring, we are confident that HRV measurements were made at a time when connectedness was gradually restored in sub-cortical structures of the patients.

## Heart rate variability as a prognosis tool for unconsciousness

Researchers have used various strategies to try to determine whether HRV measurement can distinguish between various anaesthesia states (ranging from an awake state to unconsciousness): comparison of time-domain and frequency-domain variables in two anaesthesia states, construction of new index based on HRV, etc. The Similarity Index was described in 2007 as a time-domain marker of the similarity between consecutive data and was well able to distinguish between two states of isoflurane-induced anaesthesia (probability of correct prediction: 0.91) [35]. Similarly, two anaesthesia states differ significantly in their LF and HF spectra [35]. The non-rhythmic-to-rhythmic ratio (NRR) was defined in 2014 as the ratio between the momentary non-rhythmic power and the momentary rhythmic power; it provides insight into the magnitude of controlled, ventilatory-induced arrhythmia. The NRR performed better than

<table><caption>**Table 6.** Values misclassified by the CART algorithm.</caption><thead><tr><th>Time</th><th>-30</th><th>-20</th><th>-10</th><th>-5</th><th>0</th><th>2</th><th>5</th></tr></thead><tbody><tr><td>BIS<sup>TM</sup>≥60 classified as BIS<sup>TM</sup><60</td><td>1</td><td>0</td><td>0</td><td>0</td><td>20</td><td>5</td><td>6</td></tr><tr><td>BIS<sup>TM</sup><60 classified as BIS<sup>TM</sup>≥60</td><td>0</td><td>3</td><td>5</td><td>5</td><td>0</td><td>1</td><td>0</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0310627.t006
a gold standard in predicting a response to a nociceptive stimulus [36]. The Similarity and Dis-
tribution Index (SDI) described in 2017 is based on the same mathematical method as the Sim-
ilarity Index and to which an artificial neural network model is additionally fitted. The SDI
was strongly correlated with the state of deep anaesthesia, as defined by a panel of experienced
anaesthesiologists. However, the application of the SDI is limited by the very heterogeneous
anaesthetic protocols administered, which blurs the specific relationship with the depth of
anaesthesia [37]. Other researchers have computed several HRV variables and fed them into a
deep neural network; compared with an expert assessment, the overall accuracy for the depth
of anaesthesia was 90.1%. However, the specific accuracy for awakening after GA was only
56.6%. Again, different anaesthesia protocols led to different conclusions [24].

The literature methods featured a variety of anaesthesia protocols (e.g. total intravenous anaesthesia, halogenated ethers, or even mixed protocols) and ran the models before the induction of anaesthesia and then through steady-state anaesthesia until awakening. We believe that inclusion of the induction period and the endotracheal tube placement during the modelled period is extremely challenging because the apnoea induced by induction has a strong influence on HRV variables in general and those related to parasympathetic activity (on which the influence of respiratory sinus arrhythmia is very strong) in particular [38].

Limitations

The present study had some limitations. Firstly, the prospective ANI REMI LOOP clinical trial was well documented but was not designed for the specific objective of the post-hoc analysis presented here. Particularly, the sample size was small and not computed in accordance to the primary objective of this study. These results need to be confirmed on a larger population in a prospective clinical trial. Secondly, the amount of data was smaller during the arousal period (after T0) than during the time between T-30 and T0; this might have overloaded the models with "deep anaesthesia" data, relative to "pre-arousal" data. Lastly, we only analyzed simple, well-known time-domain and spectral HRV variables. The inclusion of more complex, non-linear indexes (e.g. the approximate entropy or the sample entropy) [39] would perhaps have led to better results. However, in order to facilitate physiopathological interpretation of our results, we decided to focus on variables known to be related to various aspects of the ANS. In this study we only considered HRV features whereas taking into account the blood pressure variability would allow obtaining a much more realistic model of the baroreflex response to anaesthetic agents [21, 22]. In future work, the addition of invasive or non-invasive continuous blood pressure monitoring could allow testing this approach.

Conclusion

We found that HRV analysis provided accurate clinically useful insights into the "pre-arousal"
recovery period of GA. Further prospective large scale clinical studies are necessary to validate
these results and to confirm the ability of HRV measurements to predict AAGA.

**Supporting information**

**S1 Table. Heart rate variability measures.** BIS: Bispectral Index™, HR:heart rate, ANI: Analgesia Nociception Index, STV: short-term variation, SDNN: standard deviation of RR intervals, LTV: long term variation, RMSSD: root mean square of successive difference in RR intervals, LF: low frequency, HF: high frequency. (XLSX)
## Author Contributions

**Formal analysis:** Camille Ternynck.

**Investigation:** Anne Wojtanowski.

**Supervision:** Julien de Jonckheere.

**Validation:** Julien de Jonckheere.

**Writing – original draft:** Anne Wojtanowski, Maxence Hureau, Mathieu Jeanne, Julien de Jonckheere.

**Writing – review & editing:** Anne Wojtanowski, Maxence Hureau, Camille Ternynck, Benoit Tavernier, Mathieu Jeanne, Julien de Jonckheere.

## References

1. Pandit JJ, Andrade J, Bogod DG, Hitchman JM, Jonker WR, Lucas N, et al. 5th National Audit Project (NAP5) on accidental awareness during general anaesthesia: summary of main findings and risk factors. Br J Anaesth. 2014 Oct; 113(4):549–59. https://doi.org/10.1093/bja/aeu313 PMID: 25204697

2. Walker EMK, Bell M, Cook TM, Grocott MPW, Moonesinghe SR. Patient reported outcome of adult peri-operative anaesthesia in the United Kingdom: a cross-sectional observational study. Br J Anaesth. 2016 Dec; 117(6):758–766r. https://doi.org/10.1093/bja/aew381 PMID: 27956674

3. Kim MC, Fricchione GL, Akeju O. Accidental awareness under general anaesthesia: Incidence, risk factors, and psychological management. BJA Educ. 2021 Apr; 21(4):154–61. https://doi.org/10.1016/j.bjae.2020.12.001 PMID: 33777414

4. Tasbihgou SR, Vogels MF, Absalom AR. Accidental awareness during general anaesthesia—a narrative review. Anaesthesia. 2018 Jan; 73(1):112–22. https://doi.org/10.1111/anae.14124 PMID: 29210043

5. Macgregor K. A Waking Nightmare: How can we Avoid Accidental Awareness during General Anaesthesia? J Perioper Pract. 2013 Sep; 23(9):185–90. https://doi.org/10.1177/175045891302300902 PMID: 24245360

6. Messner M, Beese U, Romstöck J, Dinkel M, Tschaikowsky AK. The Bispectral Index Declines During Neuromuscular Block in Fully Awake Persons: Anesth Analg. 2003 Aug; 97(2):488–91. https://doi.org/10.1213/01.ANE.0000072741.78244.C0 PMID: 12873942

7. Christ B, Guerci P, Baumann C, Meistelman C, Schmartz D. Influence of neuromuscular block and reversal on bispectral index and NeuroSense values: Eur J Anaesthesiol. 2014 Aug; 31(8):437–9. https://doi.org/10.1097/EJA.0b013e32836394df PMID: 23867777

8. Akeju O, Hamilos AE, Song AH, Pavone KJ, Purdon PL, Brown EN. GABAA circuit mechanisms are associated with ether anesthesia-induced unconsciousness. Clin Neurophysiol Off J Int Fed Clin Neurophysiol. 2016 Jun; 127(6):2472–81. https://doi.org/10.1016/j.clinph.2016.02.012 PMID: 27178867

9. Soplata AE, Adam E, Brown EN, Purdon PL, McCarthy MM, Kopell N. Rapid thalamocortical network switching mediated by cortical synchronization underlies propofol-induced EEG signatures: a biophysical model. J Neurophysiol. 2023 Jul 1; 130(1):86–103. https://doi.org/10.1152/jn.00068.2022 PMID: 37314079

10. Mulvey DA, Klepsch P. Use of Processed Electroencephalography in the Clinical Setting. Curr Anesthesiol Rep. 2020 Dec; 10(4):480–7. https://doi.org/10.1007/s40140-020-00424-3 PMID: 33110400

11. Gaskell AL, Hight DF, Winders J, Tran G, Defresne A, Bonhomme V, et al. Frontal alpha-delta EEG does not preclude volitional response during anaesthesia: prospective cohort study of the isolated fore-arm technique. Br J Anaesth. 2017 Oct; 119(4):664–73. https://doi.org/10.1093/bja/aex170 PMID: 29121278

12. Casey CP, Tanabe S, Farahbakhsh Z, Parker M, Bo A, White M, et al. Distinct EEG signatures differentiate unconsciousness and disconnection during anaesthesia and sleep. Br J Anaesth. 2022 Jun; 128 (6):1006–18. https://doi.org/10.1016/j.bja.2022.01.010 PMID: 35148892

13. Myles P, Leslie K, McNeil J, Forbes A, Chan M. Bispectral index monitoring to prevent awareness during anaesthesia: the B-Aware randomised controlled trial. The Lancet. 2004 May; 363(9423):1757–63. https://doi.org/10.1016/S0140-6736(04)16300-9 PMID: 15172773

14. Avidan MS, Zhang L, Burnside BA, Finkel KJ, Searleman AC, Selvidge JA, et al. Anesthesia Awareness and the Bispectral Index. N Engl J Med. 2008 Mar 13; 358(11):1097–108. https://doi.org/10.1056/NEJMoa0707361 PMID: 18337600
15. Jeanne M, Delecroix M, De Jonckheere J, Keribedj A, Logier R, Tavernier B. Variations of the Analgesia Nociception Index During Propofol Anesthesia for Total Knee Replacement. Clin J Pain. 2014 Dec; 30 (12):1084-8. https://doi.org/10.1097/AJP.0000000000000083 PMID: 24525906

16. Boselli E, Bouvet L, Bégou G, Torkmani S, Allaouchiche B. Prediction of hemodynamic reactivity during total intravenous anesthesia for suspension laryngoscopy using Analgesia/Nociception Index (ANI): a prospective observational study. Minerva Anestesiol. 2015 Mar; 81(3):288-97. PMID: 25014481

17. Sabourdin N, Burey J, Tuffet S, Thomin A, Rousseau A, Al-Hawari M, et al. Analgesia Nociception Index-Guided Remifentanil versus Standard Care during Propofol Anesthesia: A Randomized Controlled Trial. J Clin Med. 2022 Jan 11; 11(2):333. https://doi.org/10.3390/jcm11020333 PMID: 35054027

18. Dundar N, Kus A, Gurkan Y, Toker K, Solak M. Analgesia nociception index (ani) monitoring in patients with thoracic paravertebral block: a randomized controlled study. J Clin Monit Comput. 2018 Jun; 32 (3):481-6. https://doi.org/10.1007/s10877-017-0036-9 PMID: 28631050

19. Upton HD, Ludbrook GL, Wing A, Sleigh JW. Intraoperative "Analgesia Nociception Index"-Guided Fentanyl Administration During Sevoflurane Anesthesia in Lumbar Discectomy and Laminectomy: A Randomized Clinical Trial. Anesth Analg. 2017 Jul; 125(1):81-90. https://doi.org/10.1213/ANE. 0000000000001984 PMID: 28598927

20. De Jonckheere J, Bonhomme V, Jeanne M, Boselli E, Gruenewald M, Logier R, et al. Physiological Signal Processing for Individualized Anti-nociception Management During General Anesthesia: a Review. Yearb Med Inform. 2015 Aug; 24(01):95-101. https://doi.org/10.15265/IY-2015-004 PMID: 26293855

21. Ebert TJ, Muzi M, Berens R, Goff D, Kampine JP. Sympathetic Responses to Induction of Anesthesia in Humans with Propofol or Etomidate. Anesthesiology. 1992 May 1; 76(5):725-33. https://doi.org/10.1097/00000542-199205000-00010 PMID: 1575340

22. Bari V, Vaini E, Pistuddi V, Fantinato A, Cairo B, De Maria B, et al. Comparison of Causal and Non-causal Strategies for the Assessment of Baroreflex Sensitivity in Predicting Acute Kidney Dysfunction After Coronary Artery Bypass Grafting. Front Physiol. 2019 Oct 18; 10:1319. https://doi.org/10.3389/fphys.2019.01319 PMID: 31681021

23. Huhle R, Burghardt M, Zaunseder S, Wessel N, Koch T, Malberg H, et al. Effects of awareness and nociception on heart rate variability during general anaesthesia. Physiol Meas. 2012 Feb 1; 33(2):207-17. https://doi.org/10.1088/0967-3334/33/2/207 PMID: 22260880

24. Zhan J, Wu Z xi, Duan Z xin, Yang G ying, Du Z yong, Bao X hang, et al. Heart rate variability-derived features based on deep neural network for distinguishing different anaesthesia states. BMC Anesthesiol. 2021 Mar 2; 21(1):66. https://doi.org/10.1186/s12871-021-01285-x PMID: 33653263

25. Thilen SR, Hansen BE, Ramaiah R, Kent CD, Treggiari MM, Bhananker SM. Intraoperative Neuromuscular Monitoring Site and Residual Paralysis. Anesthesiology. 2012 Nov 1; 117(5):964-72. https://doi.org/10.1097/ALN.0b013e31826f8fdd PMID: 23001053

26. Logier R, De Jonckheere J, Dassonneville A. An efficient algorithm for R-R intervals series filtering. In: The 26th Annual International Conference of the IEEE Engineering in Medicine and Biology Society [Internet]. San Francisco, CA, USA: IEEE; 2004 [cited 2024 Jul 23]. p. 3937-40. Available from: http://ieeexplore.ieee.org/document/1404100/

27. Logier R., Dassonneville A. Method and device for filtering a series of cardiac rhythm signals (rr) derived from a cardiac signal, and more particularly an ecg signal. WO2002069178A3.

28. Heart rate variability: standards of measurement, physiological interpretation and clinical use. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electro-physiology. Circulation. 1996 Mar 1; 93(5):1043-65.

29. Sauerbrei W, Schumacher M. A bootstrap resampling procedure for model building: Application to the cox regression model. Stat Med. 1992 Jan; 11(16):2093-109. https://doi.org/10.1002/sim.4780111607 PMID: 1293671

30. R Core Team (2013). R: A language and environment for statistical computing [Internet]. Vienne, Austria; Available from: http://www.R-project.org/

31. Lu K, Holzmann M, Abtahi F, Lindecrantz K, Lindqvist PG, Nordstrom L. Fetal heart rate short term variation during labor in relation to scalp blood lactate concentration. Acta Obstet Gynecol Scand. 2018 Oct; 97(10):1274-80. https://doi.org/10.1111/aogs.13390 PMID: 29799630

32. Wretler S, Holzmann M, Graner S, Lindqvist P, Falck S, Nordström L. Fetal heart rate monitoring of short term variation (STV): a methodological observational study. BMC Pregnancy Childbirth. 2016 Dec; 16(1):55. https://doi.org/10.1186/s12884-016-0845-8 PMID: 26984160

33. Aragón-Benedí C, Caballero-Lozada AF, Perez-Calatayud AA, Marulanda-Yanten AM, Oliver-Fornies P, Boselli E, et al. Prospective multicenter study of heart rate variability with ANI monitor as predictor of mortality in critically ill patients with COVID-19. Sci Rep. 2022 Dec 16; 12(1):21762. https://doi.org/10.1038/s41598-022-25537-z PMID: 36526646
34. Zickmann B, Hofmann HC, Pottkämper C, Knothe C, Boldt J, Hempelmann G. Changes in heart rate variability during induction of anesthesia with fentanyl and midazolam. J Cardiothorac Vasc Anesth. 1996 Aug; 10(5):609-13. https://doi.org/10.1016/S1053-0770(96)80138-8 PMID: 8841868

35. Huang HH, Lee YH, Chan HL, Wang YP, Huang CH, Fan SZ. Using a short-term parameter of heart rate variability to distinguish awake from isoflurane anesthetic states. Med Biol Eng Comput. 2008 Oct; 46(10):977-84. https://doi.org/10.1007/s11517-008-0342-y PMID: 18414913

36. Lin YT, Wu HT, Tsao J, Yien HW, Hseu SS. Time-varying spectral analysis revealing differential effects of sevoflurane anaesthesia: non-rhythmic-to-rhythmic ratio: Non-rhythmic-to-rhythmic ratio. Acta Anaesthesiol Scand. 2014 Feb; 58(2):157-67.

37. Liu Q, Ma L, Chiu RC, Fan SZ, Abbod MF, Shieh JS. HRV-derived data similarity and distribution index based on ensemble neural network for measuring depth of anaesthesia. PeerJ. 2017; 5:e4067. https://doi.org/10.7717/peerj.4067 PMID: 29158992

38. Nakatsuka I, Ochiai R, Takeda J. Changes in heart rate variability in sevoflurane and nitrous oxide anaesthesia: effects of respiration and depth of anesthesia. J Clin Anesth. 2002 May; 14(3):196-200. https://doi.org/10.1016/s0952-8180(01)00384-1 PMID: 12031752

39. Shaffer F, Ginsberg JP. An Overview of Heart Rate Variability Metrics and Norms. Front Public Health. 2017; 5:258. https://doi.org/10.3389/fpubh.2017.00258 PMID: 29034226