Fractal Properties of Heart Rate Dynamics: A New Biomarker for Anesthesia—Biphasic Changes in General Anesthesia and Decrease in Spinal Anesthesia

Processed electroencephalogram (EEG) has been considered a useful tool for measuring the depth of anesthesia (DOA). However, because of its inability to detect the activities of the brain stem and spinal cord responsible for most of the vital signs, a new biomarker for measuring the multidimensional activities of the central nervous system under anesthesia is required. Detrended fluctuation analysis (DFA) is a new technique for detecting the scaling properties of nonstationary heart rate (HR) behavior. This study investigated the changes in fractal properties of heart rate variability (HRV), a nonlinear analysis, under intravenous propofol, inhalational desflurane, and spinal anesthesia. We compared the DFA method with traditional spectral analysis to evaluate its potential as an alternative biomarker under different levels of anesthesia. Eighty patients receiving elective procedures were randomly allocated different anesthesia. HRV was measured with spectral analysis and DFA short-term (4–11 beats) scaling exponent (DFAα1). An increase in DFAα1 followed by a decrease at higher concentrations during propofol or desflurane anesthesia is observed. Spinal anesthesia decreased the DFAα1 and low-/high-frequency ratio (LF/HF ratio). DFAα1 of HRV is a sensitive and specific method for distinguishing changes from baseline to anesthesia state. The DFAα1 provides a potential real-time biomarker to measure HRV as one of the multiple dimensions of the DOA.


Introduction
Delivering adequate anesthetics to achieve a suitable depth of anesthesia (DOA) is an ongoing challenge for anesthesiologists. Researchers have focused on analyzing electroencephalograms (EEG) as reliable noninvasive ways to monitor the DOA [1]. However, measures derived from EEGs mainly reflect anesthetic effects on the cortex and consciousness state but not the brainstem and hemodynamic stability. The latter two are controlled by the autonomic nervous system (ANS) nuclei and are included in the multiple dimensions of DOA. Heart rate variability (HRV) has been used to assess the ANS [2][3][4][5][6] and may reflect the effects of the anesthetic on the brainstem component [7]. Propofol infusion [8][9][10][11], desflurane inhalation [12][13][14], and spinal anesthesia [15][16][17] have been reported to modulate the sympathovagal effects using spectral analysis of HRV. Because nonlinear phenomena are involved in the genesis of human heart rate fluctuations [18] and limitations of traditional linear measurement are observed, a new analytical method has been developed to evaluate cardiac regulation and to characterize the feature of heart rate (HR) dynamics [7,[19][20][21][22].

Subjects and Study Protocol
After Institutional and Ethical investigational Committee approval (NTUH IRB: 200902012R) and obtaining informed consent, we enrolled 80 healthy (ASA 1) adult patients scheduled for orthopedic or general surgery under either spinal or general anesthesia in our study. Patients with a history of ischemic heart disease, congestive heart failure, diabetes mellitus, or any other disorders known to affect autonomic nervous functions were excluded. None of the patients were taking medications that may affect cardiovascular functions.

General Patient Management
Each patient was fasted at least 8 h prior to testing. Vigorous exercise, alcohol, or coffee intake were prohibited for 48 h before the scheduled surgery. Without any premedication, on the day of surgery, the patient lay in a supine position in a quiet room at least 5 min prior to preanesthetic baseline data collection [37]. A baseline electrocardiogram (ECG) was recorded and stored on a personal computer by ECG Holter (E30-8010, Micro-Star International Co., New Taipei City, Taiwan). The sampling rate was 1000 Hz. An AS/3 Anesthesia Monitor (Datex-Engstrom, Finland) was used to collect blood pressure (BP) data, ECG, and pulse oximetry (SpO 2 ) signals from the patients. Randomization was achieved using an opaque sealed envelope technique that had been sorted by a computergenerated random allocator. Patients undergoing elective trunk or upper limb surgery were randomly allocated to two groups receiving total intravenous infusion with propofol (Group P, n = 20) or inhalational desflurane anesthesia (Group D, n = 20). DOA was continuously monitored using the A-line Index (AAI) generated by an auditory evoked potential (AEP) monitor (Danmeter, Odense, Denmark). AAI value below 35 is considered deep-plane anesthesia, and between 35 and 60, light-plane anesthesia. AEP was elicited with a binaural click stimulus of 65 dB intensity, 2 ms duration, and a repetition rate of 9 Hz. Patients undergoing elective lower limb surgery were randomly allocated to two groups receiving intrathecally either low-dose hyperbaric bupivacaine alone (Group LM, n = 20) or bupivacaine supplemented with intrathecal fentanyl (Group LMf, n = 20).

Induction and Maintenance of Intravenous Propofol Anesthesia
All patients received 100% oxygen via a facemask for 2 to 3 min prior to induction of anesthesia. Anesthesia was induced and maintained with continuous infusion of propofol at a rate of 300 µg·kg −1 min −1 via an indwelling peripheral venous catheter. Spontaneous respiration was maintained and assisted with gentle intermittent positive-pressure ventilation (IPPV) via a mask if required to maintain End-tidal carbon dioxide (ETCO 2 ) between 30 and 40 mmHg.

Induction and Maintenance of Desflurane Anesthesia
To avoid possible bronchial irritation caused by desflurane, all patients received intravenous fentanyl 1µg/kg before induction. Anesthesia was induced with an incremental increase in desflurane of 3, 6, 9, and 12% in a gas mixture of 50:50% O 2 and N 2 O. Spontaneous ventilation was maintained and assisted with IPPV via mask if required to maintain ETCO 2 between 35 and 40 mmHg.

Spinal Anesthesia
Spinal anesthesia was performed with patients placed in a lateral decubitus position with a 27-G needle inserted via L3-4 or L4-5 interspace. In Group LM, 1.2 mL of 0.5% hyperbaric bupivacaine was injected intrathecally; in Group LMf, 1.2 mL of 0.5% hyperbaric bupivacaine combined with 20 µg (0.4 mL) of intrathecal fentanyl was injected (total injectant was 1.6 mL). The block height was determined by pinprick and cold sensation test examined every 10 min until desired level T6-7 was reached. Patients who developed severe bradycardia (HR < 50/min) or hypotension (systolic BP < 80 mmHg) after spinal anesthesia requiring treatment with atropine or ephedrine were excluded from the final statistical analysis.

Data Collection
The recorded ECG signals were retrieved after surgery to measure the consecutive R-R intervals using LabChart v5 (ADInstruments, Colorado Springs, CO, USA), and all R-R intervals were edited manually to exclude all premature beats and noise. The last 660 stationary R-R intervals were obtained for DFA and spectral analysis. If the percentage of deletion was greater than 5%, the subject was excluded from the study.

Time and Frequency Domain Analysis
The parameters of HRV were calculated according to the Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology [37]. The time domain measures of HRV, which included mean R-R intervals, the standard deviation of normal-to-normal (N-N) interbeat intervals (SDNN), and the root mean square of the successive difference in N-N intervals (RMSSD), were calculated. The direct current component was excluded before the calculation of power spectra. The areas under the spectra peaks within the ranges of 0.04-0.15 Hz and 0.15-0.4 Hz were defined as lowfrequency power (LFP) and high-frequency power (HFP), respectively. Normalized highfrequency power (nHFP = 100 × HFP/(LFP + HFP)) was used as the index of cardiac vagal modulation, and the normalized low-frequency power (nLFP = 100 × LFP/(LFP + HFP)) as the index of combined sympathetic and vagal modulation ( Figure 1).

Detrended Fluctuation Analysis (DFA)
The total time series was first integrated and divided into segments of length n; each segment was then detrended by subtracting the best linear fit. The fluctuation function F(n) was then calculated as the root mean square of the detrended time series as a function of the segment size n. The α value represents the correlation properties of the time series ( Figure 2) [30,38].   segment was then detrended by subtracting the best linear fit. The fluctuation function F(n) was then calculated as the root mean square of the detrended time series as a function of the segment size n. The α value represents the correlation properties of the time series ( Figure 2) [30,38].  The scaling exponent α can be estimated by a linear fit on the log-to-log plot of F(n) versus n. The α value represents the correlation properties of the time series. The global scaling exponent a value was calculated within the range of n, between n = 4 and n = 165 space. Short-term α1 was calculated between n = 4 and n = 11.

Statistical Analysis
Statistical analysis was performed with STATA 8.2 (StataCorp, College Station, Texas, USA). We decided that a 10% difference in percentage changes of HRV parameters relative to baseline between the groups is important; therefore, n = 15 patients in each group would be necessary to detect such a difference if α = 0.05 and β = 0.1. The normality of distribution was tested using Shapiro-Wilk test. Because of the data's skewed distribution, SDNN, low-frequency/high-frequency power ratio (LF/HF ratio), nLFP, and nHFP were transformed by calculating their neutral logarithm. The differences within groups were analyzed by multivariate analysis of variance (MANOVA) with Bonferroni correction of multiple measurements. All data are expressed as the mean ± SD.
A Pearson product-moment correlation analysis was used if the data analyzed were normally distributed, or a Spearman's rank order correlation analysis was used if the data analyzed were distribution-free data in the correlation analysis between traditional measures and DFA measures.
A ROC curve was used to estimate the sensitivity and specificity of the different analysis methods in individual anesthesia states. The ROC curve reflects relative true positive, true negative, false positive, and false negative values termed specificity and sensitivity. The area under the ROC curve measures discrimination, that is, the ability of the analytical methods to correctly classify HR series data before and after anesthesia. Thus, the purpose of this comparison is to examine whether DFAα1 and LF/HF ratio would exhibit different sensitivities and specificities under different types of anesthesia.

Statistical Analysis
Statistical analysis was performed with STATA 8.2 (StataCorp, College Station, TX, USA). We decided that a 10% difference in percentage changes of HRV parameters relative to baseline between the groups is important; therefore, n = 15 patients in each group would be necessary to detect such a difference if α = 0.05 and β = 0.1. The normality of distribution was tested using Shapiro-Wilk test. Because of the data's skewed distribution, SDNN, low-frequency/high-frequency power ratio (LF/HF ratio), nLFP, and nHFP were transformed by calculating their neutral logarithm. The differences within groups were analyzed by multivariate analysis of variance (MANOVA) with Bonferroni correction of multiple measurements. All data are expressed as the mean ± SD.
A Pearson product-moment correlation analysis was used if the data analyzed were normally distributed, or a Spearman's rank order correlation analysis was used if the data analyzed were distribution-free data in the correlation analysis between traditional measures and DFA measures.
A ROC curve was used to estimate the sensitivity and specificity of the different analysis methods in individual anesthesia states. The ROC curve reflects relative true positive, true negative, false positive, and false negative values termed specificity and sensitivity. The area under the ROC curve measures discrimination, that is, the ability of the analytical methods to correctly classify HR series data before and after anesthesia. Thus, the purpose of this comparison is to examine whether DFAα1 and LF/HF ratio would exhibit different sensitivities and specificities under different types of anesthesia.

Baseline Characteristics of the Study Subjects
In the general anesthesia groups, two patients in Group D were excluded from the final statistical analysis because of excitement and arrhythmia following the induction of anesthesia. Demographic data were listed in Table 1. Under deep anesthesia with propofol (AAI < 35), BP decreased significantly, while HR remained stable (Table 1), whereas, with desflurane, both HR and BP decreased significantly from baseline (Table 1). In the spinal anesthesia groups, no patient received sedatives nor required atropine treatment. No patient had arrhythmia >5% throughout the study period. All 40 patients enrolled in the study were included in the final statistical analysis. Thirty minutes after induction of spinal anesthesia, HR decreased significantly in both groups, and BP decreased significantly only in group LMf, but not in group LM (Table 1). Unless otherwise noted, values represent the number or mean ± SD; * p < 0.05 vs. baseline; p values were calculated based on paired t-test; Group LM is spinal anesthesia with only bupivacaine; Group LMf is spinal anesthesia supplemented with intrathecal fentanyl; Group P is intravenous anesthesia with propofol infusion; Group D is inhalational anesthesia with desflurane; SA-spinal anesthesia; AAI-A-Line ARX Index; T-thoracic dermatome.

HR Variability and DFAα1 between Baseline and after Anesthesia
Spectral and DFAα1 analysis of HRV parameters at baseline and postanesthesia are shown in Table 2.
Because of the data's skewed distribution, SDNN, LF/HF ratio, nLFP, and nHFP were transformed by calculating their neutral logarithm. The differences within groups were analyzed by multivariate analysis of variance with Bonferroni correction of multiple measurements.
Because of the data's skewed distribution, SDNN, LF/HF ratio, nLFP, and nHFP were transformed by calculating their neutral logarithm. The differences within groups were analyzed by multivariate analysis of variance with Bonferroni correction of multiple measurements.
General anesthesia: in both propofol and desflurane groups, DFAα1 increased significantly under light-plane anesthesia (35 < AAI < 60) but decreased significantly under the deep plane (AAI < 35) (Figures 3 and 4). LF/HF ratio, nLFP, and LFP decreased, whereas SDNN and nHFP increased significantly from baseline in Group D under the deep-plane anesthesia (AAI < 35). RMSSD, HFP, nHFP, and LFP decreased significantly in Group P under the deep-plane anesthesia when AAI < 35. (Table 2) Figure 3. LF/HF ratio at baseline and post-general anesthesia for Group P (propofol); Group D (desflurane). Box plots show the median, 10th, 25th, 75th, and 90th percentiles as vertical boxes with error bars and plot all data points that lie outside the 10th and 90th percentiles as black circle dots. The significance levels with paired t-tests between baseline and successive measurement periods are as follows: * p < 0.05.  . LF/HF ratio at baseline and post-general anesthesia for Group P (propofol); Group D (desflurane). Box plots show the median, 10th, 25th, 75th, and 90th percentiles as vertical boxes with error bars and plot all data points that lie outside the 10th and 90th percentiles as black circle dots. The significance levels with paired t-tests between baseline and successive measurement periods are as follows: * p < 0.05. Spinal anesthesia: following spinal anesthesia, nLFP, LF/HF ratio, and DFAα1 decreased significantly, whereas nHFP increased significantly from baseline in both groups ( Figures 5 and 6). SDNN decreased significantly only in Group LMf, and HFP decreased significantly only in Group LM. (Table 2).  (desflurane). Box plots show the median, 10th, 25th, 75th, and 90th percentiles as vertical boxes with error bars and plot all data points that lie outside the 10th and 90th percentiles as black circle dots. The significance levels with paired t-tests between baseline and successive measurement periods are as follows: * p < 0.05.

Correlation between α1 and Other HR Variability Variables
DFAα1 value did not correlate with the traditional LF/HF ratio in all groups following induction of anesthesia (Table 3). RMSSD in the time domain, nLFP, nHFP, and HFP in the frequency domain correlated significantly with DFAα1 in Group P under deepplane anesthesia (AAI < 35). DFAα1 also significantly correlated with nLFP and nHFP in showing the median, 10th, 25th, 75th, and 90th percentiles as vertical boxes with error bars and plot all data points outside the 10th and 90th percentiles as black circle dots. The significance levels with paired t-tests between baseline and successive measurement periods are as follows: *** p < 0.001.

Correlation between α1 and Other HR Variability Variables
DFAα1 value did not correlate with the traditional LF/HF ratio in all groups following induction of anesthesia (Table 3). RMSSD in the time domain, nLFP, nHFP, and HFP in the frequency domain correlated significantly with DFAα1 in Group P under deep-plane anesthesia (AAI < 35). DFAα1 also significantly correlated with nLFP and nHFP in Group LMf after spinal anesthesia. * Correlation coefficient: Pearson product-moment correlation analysis for normally distributed data and Spearman rank order correlation analysis for distribution-free data; Abbreviations: see Table 2.

The ROC Curve
In the area under the curve (AUC) in spinal anesthesia Group LM and LMf, the AUC of DFAα1 were 0.84 and 0.98, respectively; the AUC of LF/HF ratio was 0.65 and 0.73, respectively ( Figure 7A). In general anesthesia Groups P and D, the AUC of DFAα1 was 0.68 and 0.85, respectively; the AUC of LF/HF ratio was 0.55 and 0.72, respectively ( Figure 7B,C).
In the area under the curve (AUC) in spinal anesthesia Group LM and LMf, the A of DFAα1 were 0.84 and 0.98, respectively; the AUC of LF/HF ratio was 0.65 and 0 respectively ( Figure 7A). In general anesthesia Groups P and D, the AUC of DFAα1 0.68 and 0.85, respectively; the AUC of LF/HF ratio was 0.55 and 0.72, respectively ( Fig  7B,C).

Discussion
In this study, we demonstrated the biphasic effects of HRV after general anesth induction with desflurane or propofol using DFA analysis compatible with the clin "paradoxical excitation" phenomenon [39,40]. Biphasic effects, that is, the increas DFAα1 followed by a decrease in DFAα1 at higher concentrations during the transi from awake to deep anesthesia, were not observed by traditional spectral HRV metho A biphasic reaction during induction of propofol anesthesia was also observe spectral EEG analyses [41,42], i.e., a paradox excitation of EEG alpha and beta power e if unconsciousness was observed. Interestingly, nonlinear EEG methods do not follo biphasic course but adequately reflect the hypnotic state [43]. Intravenous propofol i sion is known to cause a reduction in BP and HR and inhibit sympathetic [8,9

Discussion
In this study, we demonstrated the biphasic effects of HRV after general anesthesia induction with desflurane or propofol using DFA analysis compatible with the clinical "paradoxical excitation" phenomenon [39,40]. Biphasic effects, that is, the increase in DFAα1 followed by a decrease in DFAα1 at higher concentrations during the transition from awake to deep anesthesia, were not observed by traditional spectral HRV methods.
A biphasic reaction during induction of propofol anesthesia was also observed in spectral EEG analyses [41,42], i.e., a paradox excitation of EEG alpha and beta power even if unconsciousness was observed. Interestingly, nonlinear EEG methods do not follow a biphasic course but adequately reflect the hypnotic state [43]. Intravenous propofol infusion is known to cause a reduction in BP and HR and inhibit sympathetic [8,9] or parasympathetic [10,11] activity. This inhibition of autonomic activity is believed to be one of the major mechanisms underlying propofol-induced hemodynamic depression. Although it is generally agreed that propofol anesthesia is associated with a reduction in HRV [11,44,45], results regarding the effect on cardiac sympathetic or parasympathetic tone are conflicting. In the present study, we found that DFAα1 of HRV increased from the awake state to light general anesthesia (strong fractal property) and decreased from light to deep anesthesia (breakdown of fractal property). Tulppo et al. demonstrated changes in HR dynamics toward stronger short-term fractal correlation properties (α1) when ANS was activated, i.e., by decreasing vagal tone while increasing sympathetic activity, by physical stress [46]. Therefore, under light general anesthesia, increased DFAα1 is most likely due to predominant inhibition of vagal tone (decreased HFP) and concomitant inhibition but of slight, sympathetic tone. Such similar behavior of α1 has also been reported by many other physical stress studies, such as passive head-up tilt test [47], light-intensity exercise [48], and cold hand test [46]. As anesthesia becomes deeper, further inhibition of vagal and sympathetic tone occurs, causing DFAα1 to decrease. Traditional spectral analysis using LF/HF ratio did not reflect these significant changes in HRV induced by the deepening of anesthesia.
Biphasic changes, increase then decrease, in DFAα1 was also observed in desflurane anesthesia. Desflurane can cause neurocirculatory excitation manifested as an increase in HR and BP (increase in sympathetic activity) during induction and transitional stage as desflurane concentration increases [12][13][14]. In our study, under light-plane anesthesia, we did not see evidence of activation of sympathetic tone either by direct physical signs or by traditional spectral analysis of HRV; DFAα1, however, seemed able to detect the subtle changes in R-R interval dynamics better than traditional spectral analysis (Figure 3). Under deep anesthesia, both DFAα1 and LF/HF ratio decrease, which would indicate further inhibition of sympathetic tone (reduction in LFP).
The study showed DFAα1 decrease by the breakdown of the short-term fractal correlation properties of HR dynamics towards a more random pattern in spinal anesthesia. We found DFAα1 decrease is accompanied by the reciprocal increase in nHFP and a decrease in nLFP and LF/HF ratio. It also indicates that reciprocal change is associated with the breakdown of HR fractal property. The results showed enhanced vagal activity and withdrawal of sympathetic activity after spinal anesthesia had been observed in HRV [49]. Decreased DFAα1 has been observed in various disease conditions or under certain physical states, such as heart failure [23], before the onset of atrial fibrillation [29] and under cold face tests in healthy subjects [46]. However, when the local anesthetic was supplemented with fentanyl, we did not observe the synergetic effect of the sympathetic blockade as previously reported when fentanyl was intravenously administered [50]. We added 20 µg fentanyl (0.4 mL) to 1.2 mL of bupivacaine, increasing the injecting volume to 1.6 mL in patients in Group LMf. This could have made the local anesthetic less hyperbaric or even become hypobaric. However, Patterson et al. [51] have reported that the addition of fentanyl does not alter the extent of the spread of intrathecal isobaric bupivacaine. Furthermore, the blockade levels reached were equal between the two groups; we believe the increased volume and change in baricity of the injectate would have a very minimal impact on the result of our findings. The effect of intrathecal fentanyl might be confounded by bupivacaine or by its central nervous system sedative effects [52,53]. The difference in heart fractal property demonstrated by intrathecal and intravenous administration of fentanyl needs to be further investigated.
The present study did not show parameters derived from traditional spectral analysis of HRV correlate well with the DFAα1 value. The traditional spectral analysis measures the magnitude of each frequency component (oscillation mode), but the DFA method focuses on how this oscillation mode is arranged. Thus theoretically, DFA scaling exponents will not correlate with most HRV parameters derived from spectral analysis. The weak correlation during anesthesia may partially be due to uncontrolled breathing rate since some studies have reported a high correlation between DFA and LF/HF ratio in controlled ventilation [47,48]. However, this is still unclear and is a subject of controversy, as some studies are reporting a low correlation even with controlled ventilation [54] or a high correlation without controlled ventilation [55,56]. In our study, the rates of spontaneous breathing were similar between spinal and general anesthesia groups, 10-15 breaths/min under spinal and 8 −15 breaths/min under general anesthesia. This seemed to suggest other factors, such as differences in frequency characteristics, stationary assumption of spectral analysis, and BP may be involved in the inconsistency. In addition to respiratory rate, other respiratory parameters such as tidal volume [57] and possible mental stress during controlled breathing could also be important factors that affect the HRV. Several nonlinear methods for HRV have tried to distinguish awake from sleep state under anesthesia [58] and detect subtle sympathetic and parasympathetic-mediated alteration [7]. In our study, we suggest that DFAα1 in HRV is a sensitive and specific index to be used in clinical practice. Though this method could not provide direct information regarding the modulation of ANS, it may reflect the response of ANS, which may be considered one of the multiple dimensions of DOA. Moreover, our method provides real-time measurement of DOA (less than 1 s), which is an indispensable requisite for monitoring.
There are several limitations of our study. First, we did not include patients with heart disease or heart failure because we wanted to investigate how anesthetic drugs interact with HRV without avoidable confounding factors [23,25,28,29]. Second, anesthesia-induced changes in respiratory rate and tidal volume could have influenced HRV. Thus, we would have underestimated the effects of anesthetics on HRV should induction of anesthesia have resulted in a decrease in both respiratory rate and tidal volume. Third, spinal anesthesia is accompanied by significant sedation [49,59,60]. Thus, both direct blockade of cardiac sympathetic fibers and sedation caused by spinal anesthesia might have exaggerated the sympathovagal effect of spinal anesthesia. Fourth, despite equivalent AAI value numbers, it is uncertain whether the depths of anesthesia achieved were equal between propofol and desflurane anesthesia. However, there is no convincing evidence that the AAI index number was independent of the DOA in our study. Since the AAI value falls in a similar pattern as anesthesia deepened with these two agents, we can reasonably believe that the DOA achieved was equal.
In the future, we plan to employ the DFA method to study HRV behavior in heart disease patients and other anesthetics, such as ketamine or opioids, as both agents are known to have very different effects on the central nervous system and/or autonomic nervous system modulation of the cardiovascular system.

Conclusions
Biphasic effects of HRV followed by deepening propofol or desflurane anesthesia were observed. DFAα1 decreased after spinal anesthesia, even with intrathecal bupivacaine alone or supplemented with fentanyl. Light general anesthesia resulted in the change in HR dynamics toward stronger short-term fractal correlation properties, while deep general anesthesia resulted in the breakdown of the fractal property. DFAα1 change may reflect the indirect response of ANS and provide a potential real-time nonlinear method to measure HR dynamic and DOA.