Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Conf Proc IEEE Eng Med Biol Soc. Author manuscript; available in PMC 2010 January 11.
Published in final edited form as:
PMCID: PMC2804255

Linear and Nonlinear Quantification of Respiratory Sinus Arrhythmia during Propofol General Anesthesia


Quantitative evaluation of respiratory sinus arrhythmia (RSA) may provide important information in clinical practice of anesthesia and postoperative care. In this paper, we apply a point process method to assess dynamic RSA during propofol general anesthesia. Specifically, an inverse Gaussian probability distribution is used to model the heartbeat interval, whereas the instantaneous mean is identified by a linear or bilinear bivariate regression on the previous R-R intervals and respiratory measures. The estimated second-order bilinear interaction allows us to evaluate the nonlinear component of the RSA. The instantaneous RSA gain and phase can be estimated with an adaptive point process filter. The algorithm’s ability to track nonstationary dynamics is demonstrated using one clinical recording. Our proposed statistical indices provide a valuable quantitative assessment of instantaneous cardiorespiratory control and heart rate variability (HRV) during general anesthesia.

I. Introduction

Heart rate variability (HRV), the beat-to-beat variation in heart rate (HR), is considered to be a reflection of autonomic nervous system activity. Respiratory sinus arrhythmia (RSA) is thought by many to reflect the parasympathetic component of the autonomic nervous system during spontaneous ventilation [1]. HR increases with inhalation and decreases with exhalation. Therefore, RSA can be considered to be an important index of vagal control of HR and an indicator of autonomic nervous system activity processed in the nucleus solitarius of the brainstem. In clinical practice, RSA has been shown to be related to clinical signs of depth of anesthesia and has been proposed to have potential as a diagnostic monitoring measure for the depth of anesthesia [11],[17],[4] as well as cardiac or neurologic dysfunction [14],[13]. It was found in animals that morphine with α-chloralose-urethane significantly increased RSA, whereas thiopental and halothane both significantly decreased RSA [12] Studies on children during the induction of general anesthesia with halothane and nitrous oxide showed that RSA decreased in a significant manner, corresponding from awake, to loss of pharyngeal tone, to stage 3 or deep anesthesia [4].

Quantitive assessment of RSA poses a challenging statistical signal processing problem in biomedical engineering. In the literature, time domain-based sequence methods and frequency domainbased spectral methods have been proposed to evaluate RSA [10]. However, all of these methods are based on the assumption that signals are stationary or quasi-stationary within a window. Previously, we have proposed a point process method to evaluate instantaneous RSA during a non-stationary environment [5],[7]. Our previous investigation of RSA in a pharmacological blockade protocol [7] validated the linear RSA gain estimates and suggested that suppression of parasympathetic activity using atropine reduces RSA, whereas the blockade of sympathetic activity with propranolol has little effect on the RSA gain. Here, we extend the previous model to include a 2nd-order Volterra bilinear expansion to characterize the nonlinear interaction between heartbeat interval and respiratory signals for a healthy subject during general anesthesia [8] where the anesthetic drug propofol was administered with different concentrations [18]. Our preliminary study provides a firsthand quantitative assessment of RSA under propofol anesthesia for healthy subjects using a rigorous statistical estimation method.

II. Heartbeat Interval Point Process Model

Given a set of R-wave events {uj}j=1J detected from the electrocardiogram (ECG), let RRj = ujuj−1 > 0 denote the jth R-R interval. By treating the R-waves as discrete events, we may develop a probabilistic point process model in the continuous-time domain. Assuming history dependence, the waiting time tuj (as a continuous random variable) until the next R-wave event can be modeled by an inverse Gaussian model [2],[3],[5]:


where uj denotes the previous R-wave event occurred before time t, θ > 0 denotes the shape parameter, and μRR (t) denotes the instantaneous R-R mean. When the mean is much greater than the variance, the inverse Gaussian can be well approximated by a Gaussian model with a variance equal to μRR3/θ. In point process theory, the inter-event probability p(t) is related to the conditional intensity function (CIF) λ(t) by a one-to-one transformation: λ(t)=p(t)1ujtp(τ)dτ. The estimated CIF can be used to evaluate the goodness-of-fit of the probabilistic heartbeat model.

In modeling the instantaneous R-R interval mean μt, we propose to use the following two models:

  • A linear bivariate autoregressive (AR) model that reflects a bivariate discrete-time linear system:
    where the first two terms represent the AR model of the past R-R intervals (RR), and RPt j denotes the value of the respiratory (RP) measure, which is synchronously sampled at the beats prior to time t.
  • A nonlinear bivariate Volterra model that reflects a bivariate bilinear system:

The above two models can be derived within the general framework of Volterra-Wiener theory in nonlinear system identification [16].

A. Instantaneous Indices of HR and HRV

Heart rate is defined as the reciprocal of the R-R intervals. For time t measured in seconds, the new variable r = c(tuj)−1 (where c = 60 s/min) can be defined in beats per minute (bpm). By the change-of-variables formula, the HR probability p (r) = p (c(tuj)−1) is given by p(r)=[mid ]dtdr[mid ]p(t), and the mean and the standard deviation of HR r can be derived [2]:


where [mu] = c−1μRR and [theta w/ tilde] = c−1θ. Essentially, the instantaneous indices of HR and HRV are characterized by the mean μHR and standard deviation σHR, respectively.

B. Adaptive Point Porcess Filtering

Let ξ=[{ai}i=0p,{bj}j=1p,{hij},θ]T denote the vector that contains all unknown parameters in the probabilistic model, we can recursively estimate them via adaptive point process filtering [3]:

ξk[mid ]k1=ξk1[mid ]k1Pk[mid ]k1=Pk1[mid ]k1+Wξk[mid ]k=ξk[mid ]k1+Pk[mid ]k1([nabla]logλk)[nkλkΔ]Pk[mid ]k=[Pk[mid ]k11+[nabla]λk[nabla]λkTΔλk[nabla]2logλk[nkλkΔ]]1

where P and W denote the parameter and noise covariance matrices, respectively; and Δ=5 ms denotes the time bin size. Diagonal noise covariance matrix W that determines the level of parameter fluctuation at the timescale of Δ can be initialized either empirically from a random-walk theory or estimated by a maximum likelihood estimate. Symbols [nabla]λk=[partial differential]λk[partial differential]ξk$ and [nabla]2λk=[partial differential]2λk[partial differential]ξk[partial differential]ξkT denote the first- and second-order partial derivatives of the CIF w.r.t. ξ at time t = kΔ, respectively. The indicator variable nk = 1 if a heart beat occurs in time ((k − 1)Δ, kΔ] and 0 otherwise.

C. Quantification of RSA

From equations (1) or (2), we can derive the transfer function and frequency response between RP (input) and RR (output). Since the RSA effect is frequency dependent, we propose the following measure to quantify the RSA in the frequency domain [7]:

H12(f)=j=1qbj(k)zj[mid ]z=ej2πf1i=1pai(k)zi[mid ]z=ej2πf,

where f denotes the rate for the RR and RP measurements (the samples of both series are assumed to be synchronized). With the estimated time-varying AR coefficients {ai(k)} and {bj(k)} at time t = kΔ, we may evaluate the dynamic frequency response of (4) at high-frequency (HF) range (0.15–0.5 Hz, which is the range of respiration rhythm). The RSA gain, characterized by |H12 (f)|, represents the effect of RP on heartbeat. Given the baroreflex gain, we can estimate the cross-spectrum between RP (u) and RR (y) as An external file that holds a picture, illustration, etc.
Object name is nihms162133ig1.jpg (f) = H12 (f)An external file that holds a picture, illustration, etc.
Object name is nihms162133ig2.jpg (f). When the coefficients {ai(t)} and {bj(t)}are iteratively updated, the point process filter produces an assessment of instantaneous (parametric) RSA gain, as well as the cross-spectrum, at a very fine temporal resolution and without using the window technique.

In addition, in order to characterize the nonlinear coupling between RP (u) and RR (y) in the frequency domain, we can compute the cross-bispectrum [8]:


where H(f1,f2)=k=1ql=1qhklej2kπf1ej2lπf2 denotes the Fourier transform of the 2nd-order kernel coefficients {hkl}, and the instantaneous R-R spectrum is

QRR(f,t)=σRR2(t)[mid ]1i=1pai(t)zi|z=ej2πf.

Let h(t) denote a vector that contains all of 2nd-order coefficients {hkl(t)}; in light of (5), we may compute an instantaneous index that quantifies the fractional contribution between the cross-spectrum and the cross-bispectrum [8]:

ρt=[mid ]Cuy(f,t)[mid ][mid ]Cuy(f,t)[mid ]+[mid ]Cuuy(f1,f2,t)[mid ]11+2[mid ]h(t)[mid ]·[mid ]QRR(f,t)[mid ],

where |·| denotes either the norm of a vector or the modulus of a complex variable. The “≈” is due to a Gaussian assumption used in deriving (5). A small value of ρ implies a presence of significant (nonzero) values in {hkl} (i.e. nonlinearity), whereas a perfect linear Gaussian model would imply ρ = 1.

III. Experimental Protocol and Setup

The present pilot study and experimental protocol was approved by the Massachusetts General Hospital (MGH). Any subject whose medical evaluation was not classified as American Society of Anesthesiologists Physical Status I was excluded from the study. Intravenous and arterial lines were placed in each subject. Propofol was infused intravenously using a previously validated computercontrolled delivery system running STANPUMP [19] connected to a Harvard 22 syringe pump (Harvard Apparatus, Holliston, MA). Six effect-site target concentrations (0–5 mcg/ml) were each maintained for 15 minutes respectively. Capnography, pulse oximetry, ECG, and arterial blood pressure (BP) were recorded (at 1 kHz sampling rate) and monitored continuously by an anesthesiologist throughout the study. Bag-mask ventilation with 30% oxygen was administered as needed in the event of propofol-induced apnea. Since propofol is a potent peripheral vasodilator, phenylephrine was administered intravenously to maintain mean arterial BP within 20% of the baseline value [18]. The respiratory signals were simultaneously sampled and recorded at 1 kHz. The quantitative analysis of the BP and baroreflex sensitivity has been reported elsewhere [8].

In the present study, as acquisition sessions are still ongoing, we focus on one specific subject. Figure 1 shows representative R-R and respiratory recordings from this subject across 6 epochs. Specifically, during the experiment, the phenylephrine was administered at around 2960 s after the recording onset (which is at the drug effect-site target concentration of 2 mcg/ml), and it was turned off at around 7580 s. In addition, hand ventilation started at around 3130 s and ended at around 7555 s.

Fig. 1
Selected snapshots of raw R-R and non-calibrated respiratory (RP, arbitrary unit) recordings from one subject during 5 consecutive epochs (with gradually increasing levels of drug concentration from 0 to 5 mcg/ml propofol).

For the linear model (2), the bivariate orders p and q were fitted from 2 to 8 and the optimal order was chosen according to the Akaike information criterion (AIC). For the bilinear model (3), the order r = 2 was empirically chosen to avoid demanding computation burden, and the initial hij was estimated by fitting the residual error via least-squares. After computing the CIF, the goodness-of-fit of the probabilistic model for the heartbeat interval is evaluated with the Kolmogorov-Smirnov (KS) test and autocorrelation independence test [2]. To assess the RSA, we computed the frequency response (4) within the range of the respiratory frequency±0.15 Hz. For the majority of epochs, the respiratory frequency stays around 0.25~0.3 Hz, which is also the frequency band in which the RR and BP time series achieved the maximum coherency. Note that, since the respiratory measures are not calibrated, the RSA unit is arbitrary; however, we still can measure the relative change in RSA as compared to the awake baseline (level 0) for any specific subject.

IV. Results and Discussion

The main goal of this paper is to demonstrate that our proposed point process method is well suited for analyzing data from clinical studies under general anesthesia, and that it provides a valid model capable of tracking instantaneous HR, HRV, and RSA indices in a non-stationary environment. The model was validated by well established goodness-of-fit tests [3]. Specifically for the considered subject, the KS plot lies almost entirely within the 95% confidence bounds; meanwhile, the autocorrelation function has most autocorrelation lag estimates within the 95% confidence bounds, suggesting that their transformed times are nearly independent.

During general anesthesia, temporary interventions are required to stabilize vital signs. Clinical interventions such as phenylephrine administration and artificial ventilation elicit important transient responses and recurrent nonstationarities along the experiment. These kinds of procedures, which are routine in the operating room environment, require an instantaneous dynamic assessment such as the one we propose. More clearly, we show two explicit instances in Fig. 2. In the top panel we are between level 0 and level 1, and at around 2010 s propofol administration was initiated. Parallel to a sudden drop in mean RR we observe an increase in variability accompanied by a drop in RSA. A second sharper drop (both in mean RR and RSA) occurs 40 s later, then RSA values recover and stabilize around a lower level than without anesthetic effect. In the second bottom panel, at around 2960 s, phenylephrine was administered to restore blood pressure closer to baseline values. Consequently, a gradual drop in HRV is accompanied by a more rapid drop in RSA. As values tend to stabilize, hand ventilation was started at around 3125 s to compensate for a reduced respiration (see RP trace), eliciting a sudden increase in HRV accompanied by a sharp decrease in RSA. Note the marked oscillatory variations of RSA during hand ventilation, which could not have been observed with stationary window-based estimates. Also note that even in such a highly non-stationary environment, the overall goodness-of-fit of our model is still excellent (Fig. 3), although we also observed a slight mismatch in tracking μRR right after certain transient effects (e.g., ventilation in Fig. 2).

Fig. 2
Two snapshot examples of dynamic tracking (using a linear model) for instantaneous HR, HRV, and RSA (dashed, dash-dot and solid lines mark the onset time of propofol anesthesia, phenylephrine, and ventilation, respectively). The red trace in the top panel ...
Fig. 3
Goodness-of-fit tests by KS plot and autocorrelation plot. The line or dots falling within 95% confidence bounds (dashed line) indicate a good fit of the probability model for heartbeat intervals.

In Table I, the HR, HRV, LF/HF ratio, RSA gain, and the coherence (around respiration frequency), time-averaged over each 15- minute propofol level, are reported. Specifically, the time-averaged HR initially increased up to level 3 and then decreased; the time-averaged HRV gradually decreased as the drug concentration level increased; the LF/HF ratio did not correlate with drug level and showed a maximum during level 3 and a minimum during level 1; and the RSA gain generally showed a decreasing trend from level 1 to level 5. It should be stressed that these trends were not statistically significant. Furthermore, transient changes such as the one observed in our example may significantly affect reliability of the averaged statistics within the 15-min epochs of interest. In our case, this issue becomes even more relevant as a consequence of the transient effects elicited by phenylephrine administration and hand ventilation (Fig. 2).

Table I
Comparison of Mean Statistics at 6 Different Levels of Propofol Drug Concentrations.

Next, we investigate the nonlinear component of the RSA effect. Our investigation is motivated by a previous study, which suggested that RSA fluctuations during anesthesia contain nonlinear dynamic mechanisms [15]. Essentially, an increase/decrease of coherencevalues reflects the increase/decrease of linear dependence in frequency domain between two time series, which is also reflected in the cross-spectrum/cross-bispectrum ratio ρ.

In Table I, we can see that both the coherence and crossspectrum/cross-bispectrum ratio ρ show similar mixed trends across the different drug levels, suggesting that the degree of RSA nonlinearity is not a function of propofol drug level. It is interesting to point out that in our previous investigation [8] we observed an increase of nonlinearity in heartbeat interval dynamics from baseline to anesthesia, where the nonlinearity involved the bilinear interactions between RR and systolic blood pressure (SBP) accompanied by a significant decrease in linear coherence between these two series. This seems to suggest that the nonlinear component of heartbeat interval dynamics during anesthesia is mainly contributed from the cardiovascular (baroreflex) loop, whereas the linear interaction within the cardiorespiratory loop roughly remains unchanged. It is also possible that the respiratory system indirectly influences HR by modulating the baroreceptor and chemoreceptor input to cardiac vagal neurons. However, in our experimental condition, it is difficult to separate the influence of SBP from the influence of respiration on HRV.

V. Conclusion

We apply and validate a point process method for dynamic assessment of RSA during propofol general anesthesia. The proposed point process method enables us to simultaneously estimate instantaneous HR, HRV, cross-spectrum, and cross-bispectrum, all of which may serve as useful noninvasive indicators in clinical practice. It should be emphasized that since the 15-min time-averaged statistics may not give reliable estimates, our main point in this study is to demonstrate the ability to provide instantaneous assessments of cardio-respiratory coupling non-stationary events due to transient interventions as well as sudden changes in physiological state that often occur in the operating room environment. Overall, our dynamic estimates suggest that (i) RSA gradually decreases from awake baseline after administration of propofol anesthesia, although whether the RSA gain may accurately measure efferent vagal activity is still controversial [9]; (ii) the RSA effect is suppressed by the phenylephrine; and (iii) the nonlinear interactions within the cardiorespiratory control remain relatively stable. It is known that different anesthetic drugs have different impacts on autonomic control and cardiac vagal tone [12]. Complementary to our previous study [8], it would be interesting to further investigate the joint interactions between the cardiovascular and cardiorespiratory couplings and their impact on HRV. Specifically, RSA is likely to be mediated by phasic withdrawal of vagal efferent activity resulting from the mechanisms of baroreflex response to spontaneous BP fluctuations or respiratory gating of central arterial baroreceptor and chemoreceptor afferent inputs [20].

To conclude, our statistical method and quantitative analysis may also serve as a potential measure for diagnosis of a variety of cardiovascular diseases, such as hypertension, myocardial infarction, and heart failure, all of which typically have abnormal RSA.


This work was supported by NIH Grants R01-HL084502, K25-NS05758, DP1-OD003646, and R01-DA015644.

The authors are with the Neuroscience Statistics Research Laboratory, Massachusetts General Hospital, Harvard Medical School, Boston, MA 02114, USA. E. N. Brown is also with the Harvard-MIT Division of Health Science and Technology and the Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA.

The authors thank L. Citi for assistance in preprocessing the data used in our experiment.


1. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. Heart rate variability. Circulation. 1996;93(5):1043–1065. [PubMed]
2. Barbieri R, Matten EC, Alabi AA, Brown EN. A pointprocess model of human heartbeat intervals: new definitions of heart rate and heart rate variability. Am J Physiol Heart Cicr Physiol. 2005;288:424–435. [PubMed]
3. Barbieri R, Brown EN. Analysis of heart beat dynamics by point process adaptive filtering. IEEE Trans Biomed Eng. 2006;53(1):4–12. [PubMed]
4. Blues CM, Pomfrett CJD. Respiratory sinus arrhythmia and clinical signs of anaesthesia in children. British J Anaesthesia. 1998;81(3):333–337. [PubMed]
5. Chen Z, Brown EN, Barbieri R. A study of probabilistic models for characterizing human heart beat dynamics in autonomic blockade control. Proc. ICASSP’08; Las Vegas, USA. pp. 481–484. [PMC free article] [PubMed]
6. Chen Z, Brown EN, Barbieri R. Characterizing nonlinear heartbeat dynamics within a point process framework. Proc. IEEE EMBC’08; Vancouver, Canada. pp. 2781–2784. [PMC free article] [PubMed]
7. Chen Z, Brown EN, Barbieri R. Assessment of autonomic control and respiratory sinus arrhythmia using point process models of human heart beat dynamics. IEEE Trans Biomed Eng. 2009;56(7) [PMC free article] [PubMed]
8. Chen Z, Purdon PL, Pierce ET, Harrell G, Brown EN, Barbieri R. Assessment of baroreflex control of heart rate during general anesthesia using a point process method. Proc. ICASSP’09; April 19–24; Taiwan. [PMC free article] [PubMed]
9. Dexter F, Ben-Haim S. Theoretical analysis predicts that respiratory sinus arrhythmia does not accurately measure efferent vagal activity during anesthesia. J Theo Biol. 1994;169:133–141. [PubMed]
10. Dinh TP, Perrault H, Calabrese P, Eberhard A, Benchetrit G. New statistical method for detection and quantification of respiratory sinus arrhythmia. IEEE Trans Biomed Eng. 1999;46(9):1161–1165. [PubMed]
11. Donchin Y, Feld JM, Porges SW. Respiratory sinus arrhythmia during recovery from isoflurane-nitrous oxide anesthesia. Anesth Analg. 1985;64:811–815. [PubMed]
12. Halliwill JR, Billman GE. Effect of general anesthesia on cardiac vagal tone. Am J Physiol Heart Circ Physiol. 1992;262:H1719–H1724. [PubMed]
13. Huang CJ, Kuok CH, Kuo TBJ, et al. Pre-operative measurement of heart rate variability predicts hypotension during general anesthesia. Acta Anaesthesiol Scand. 2006;50:542–548. [PubMed]
14. Loula P, Jäntti V, Yli-Hankala A. Respiratory sinus arrhythmia during anaesthesia: assessment of respiration related beat-to-beat heart rate variability analysis methods. Int J Clinical Monitoring and Computing. 1997;14(4):241–249. [PubMed]
15. Loula P, Jäntti V, Yli-Hankala A. Nonlinear interpretation of respiratory sinus arrhythmia in anesthesia. Methods of Information in Medicine. 1994;33(1):52–57. [PubMed]
16. Nikias C, Petropulu AP. Higher Order Spectra Analysis: A Non- Linear Signal Processing Framework. Prentice Hall; 1993.
17. Pomfrett CJD, Barrie JR, Healy TEJ. Respiratory sinus arrhythmia: an index of light anaesthesia. British J Anaesthesia. 1993;71:212–217. [PubMed]
18. Purdon PL, Pierce ET, Bonmassar G, et al. Simultaneous electroencephalography and functional magnetic resonance imaging of general anesthesia. Annals New York Acad Sci. 2009 (in press) [PMC free article] [PubMed]
19. Shafer A, Doze VA, Shafer SL, et al. Pharmacokinetics and pharmacodynamics of propofol infusions during general anesthesia. Anesthesiology. 1988;69:348–356. [PubMed]
20. Tzeng YC, Sin PYW, Galletly DC. Human sinus arrhythmia: inconsistencies of a teleological hypothesis. Am J Physiol Heart Circ Physiol. 2009;296:H65–70. [PubMed]