Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Ann Biomed Eng. Author manuscript; available in PMC 2012 January 1.
Published in final edited form as:
PMCID: PMC3010293

Dynamic Assessment of Baroreflex Control of Heart Rate during Induction of Propofol Anesthesia Using a Point Process Method


In this paper, we present a point process method to assess dynamic baroreflex sensitivity by estimating the baroreflex gain as focal component of a simplified closed-loop model of the cardiovascular system. Specifically, an inverse Gaussian probability distribution is used to model the heartbeat interval, whereas the instantaneous mean is identified by linear and bilinear bivariate regressions on both the previous R-R intervals (RR) and blood pressure (BP) beat-to-beat measures. The instantaneous baroreflex gain is estimated as the feedback branch of the loop with a point-process filter, while the RR→BP feedforward transfer function representing heart contractility and vasculature effects is simultaneously estimated by a recursive least-squares (RLS) filter. These two closed-loop gains provide a direct assessment of baroreflex control of heart rate. In addition, the dynamic coherence, cross-bispectrum, and their power ratio can also be estimated. All statistical indices provide a valuable quantitative assessment of the interaction between heartbeat dynamics and hemodynamics. To illustrate the application, we have applied the proposed point process model to experimental recordings from eleven healthy subjects in order to monitor cardiovascular regulation under propofol anesthesia. We present quantitative results during transient periods, as well as statistical analyses on steady state epochs before and after propofol administration. Our findings validate the ability of the algorithm to provide a reliable and fast-tracking assessment of baroreflex sensitivity (BRS), and show a clear overall reduction in baroreflex gain from the baseline period to the start of propofol anesthesia, confirming that instantaneous evaluation of arterial baroreflex control of heart rate may yield important implications in clinical practice, particularly during anesthesia and in postoperative care.

Keywords: Baroreflex control, Baroreflex sensitivity, Heart rate variability, Hemodynamics, Point processes, Adaptive filters, Volterra series, Closed-loop feedback control, Cardiovascular system


A change of heart rate (HR), or beat-to-beat interval, in response to a change in arterial blood pressure (BP), provides a noninvasive measure of the sensitivity of the baroreceptor-cardiac reflex (baroreflex), which is essential in characterizing cardiovascular autonomic control and explaining the complex interactions between heartbeat dynamics and hemodynamics.27 Baroreflex sensitivity (BRS) assessment provides a valuable measure of cardiovascular regulation in normal or diseased states.26, 11, 44 The evaluation of baroreflex control of heart rate during or after general anesthesia has important implications for clinical practice and patient safety. 9, 12, 58, 60, 41, 52, 53 General anesthetic drugs are known to have direct effects on vascular tone and myocardial contractility, but little is known about how they influence cardiovascular regulation. The cardiovascular side-effects of general anesthetic agents can be serious and potentially life-threatening, particularly in very critically ill patients. Despite the importance of understanding the underlying physiological mechanism and its clinical value, few studies in the literature quantify the effect of anesthetic drugs to on the alteration of cardiovascular control under general anesthesia. Propofol (2,6-diisopropylphenol) is a lipid soluble intravenous anesthetic agent. It is widely used for the induction and maintenance of general anesthesia, as well as sedation.

One of the earliest studies quantifying the baroreceptor reflex under propofol anesthesia demonstrated a “resetting” of the baroreflex. In this study, patient volunteers were anesthetized with intravenous infusions of propofol combined with nitrous oxide. It was found that in the steady state the subjects maintained both low blood pressure and low heart rate. The authors concluded that this was the result of a “resetting” of the baroreflex, but that there was no impairment of baroreflex sensitivity.22 In contrast, in a study that investigated surgical subjects stimulated with mircrolaryngoscopy, it was found that under propofol anesthesia, in addition to the inhibition of sympathetic nervous activity in the periphery, the sensitivity of the baroreflex was decreased.58 Another study found that the sensitivity of the baroreceptors was depressed by propofol infusion during general anesthesia, lasting for up to 60 minutes after the discontinuation of the propofol infusion in 13 healthy human volunteer subjects.53 It should be emphasized that in these studies, BRS was estimated using a simple sequence method based on linear regression analysis; however, the baroreceptor gain function is known to be frequency-dependent, and the influence of vasculature and heart contractility on BP has to be simultaneously considered for a correct identification of the closed-loop cardiovascular control.3 Quantifying BRS under general anesthesia is an even more challenging statistical signal processing problem, because BRS can change rapidly in time as a result of anesthetic drug effects, compensatory maneuvers or pharmacological interventions, with not easily predictable physiological responses.

The idea of estimating time-varying transfer function or coherence function of physiological systems is not new (see reviews63, 49, 39). In the literature, several methods have been proposed to estimate time-varying transfer or coherence function, such as the Kalman filter, recursive leasts-quares (RLS), or time-varying optimal parameter search (TVOPS) algorithms.3, 67, 65, 64 In a large body of reported work of our own,5, 13, 17, 18 we have previously applied probabilistic point process models for estimating instantaneous indices of HR, HR variability (HRV), as well as respiratory sinus arrhythmia (RSA). By “instantaneous”, we mean that the statistics can be estimated in principle at any time point with arbitrarily fine time resolution, without resorting to approximation by interpolation. The point process framework enabled us to estimate these physiological indices in a dynamic fashion at a millisecond timescale. Since the cardiovascular system presents several closed-loop interactions between many variables, including R-R interval and BP, research efforts have been devoted to quantifying BRS by estimating the baroreflex gain with a closed-loop system identification approach.6, 38, 2, 3, 29, 30 In this paper, we extend the point process model of the heartbeat interval to include BP as a covariate, allowing for a dynamic assessment of the baroreflex gain within the feedback BP→RR transfer function. At the same time, a discrete-time RLS filter is used to track the parameters for estimating the RR→BP frequency response in the feedforward cardiovascular pathway. This two-filter estimation structure provides a direct assessment of the baroreflex control of heart rate, as well as an instantaneous quantification of beat-to-beat cardiovascular variability, in an online fashion.14 It is noted that the main focus of this paper is to illustrate the strength of the point process method as potential application in clinical anesthesiology, a thorough and systematic comparison between our method and other linear/nonlinear ARX methods is beyond the scope of the current paper. Nevertheless, as an illustration, we compare the tracking performance of the point process filter and the standard RLS filter in analyzing known transient dynamics of BRS. Another important issue of cardiovascular modeling is goodness-of-fit. The current point process approach provides a rigorous statistical framework for model evaluation, which is not available for the standard linear/nonlinear AR or ARX type methods (based on R-R intervals).

In healthy subjects, the heartbeat interval dynamics are known to be nonlinear or even possibly chaotic.45, 48 In a previous investigation,18, we have modeled the nonlinear heartbeat dynamics within the point process framework using the beat intervals alone. In the present study, in order to characterize a potential nonlinear interaction between the beat intervals and blood pressure measures, we model the heartbeat interval mean using a bilinear system. The use of the bilinear system identification also allows us to estimate the dynamic cross-bispectrum between the R-R and BP series, as well as the power ratio between the cross-spectrum and cross-bispectrum. We apply our point process model to experimental physiological recordings of eleven healthy subjects during induction of propofol anesthesia,51 and we conduct quantitative assessment of baroreflex control during both transient periods as anesthesia is initiated, and performing statistical analyses on steady-state epochs before and after propofol administration.


A Probability Model for the Heartbeat Interval

Given a set of R-wave events {uj}j=1J detected from the electrocardiogram (ECG) waveform, 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, where t > uj) until the next R-wave event can be modeled by an inverse Gaussian model:4, 5, 13


where uj denotes the previous R-wave event occurred before time t, θ > 0 denotes the shape parameter (which might also be time-varying), and µRR(t) denotes the instantaneous R-R mean parameter. Note that when the mean µRR(t) is much greater than the variance, the inverse Gaussian can be well approximated by a Gaussian model with a variance equal to µRR3(t)/θ. The Gaussian approximation will be used for some derivations presented later. 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:10 λ(t)=p(t)1ujtp(τ)dτ. The estimated CIF can be used to evaluate the goodness-of-fit of the probabilistic heartbeat interval model.

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 is a constant) can be defined in beats per minute (bpm). By virtue of the change-of-variables formula, from Eq. (1) the HR probability p(r) = p(c(tuj)−1) is given by p(r)=|dtdr|p(t), and the mean and the standard deviation of HR r can be derived:4

μHR=μ˜1+θ˜1,  σHR=(2μ˜+θ˜)/μ˜θ˜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. In a non-stationary environment, where the probability distribution of HR is possibly slowly changing over time, we aim to dynamically estimate the instantaneous mean μRR(t) and instantaneous shape parameter θt in Eq. (1) so that the evolution of the probability density p(r) can be tracked in an online fashion.

Modeling the Instantaneous Heartbeat Interval’s Mean

A common methodological way to understand a biological or physiological system is through system identification.40 In general, let us consider a causal, continuous-time nonlinear mapping F between an output variable y(t) and two input variables x(t) and u(t). Expanding the Wiener-Volterra series of function F (up to the second order) with respect to inputs x(t) and u(t) yields:55


where F(·) : R2R, and a(·),b(·),h1(·,·),h2(·,·), and h3(·,·) are Volterra kernels with appropriate orders. In our case of system identification, y(t) will be replaced by μRR(t), x(t) will be replaced by previous R-R intervals, u(t) will be replaced by BP, and the continuous-time integral will be approximated by a finite and discrete summation.

Let us consider two specific cases of discrete-time Volterra series expansion:

  • Case (a): Dropping off all of second-order terms in the Volterra series expansion (3), we obtain a bivariate discrete-time linear system:
    where the first two terms represent a linear autoregressive (AR) model of the past R-R intervals, a0 compensates the nonzero mean effect of the R-R measurements, and BPtj denotes the previous jth BP value prior to time t. The BP in Eq. (4) can be either the systolic or the diastolic value, and it will be preprocessed to have zero mean (since the DC component is of no interest to model the oscillation). Eq. (4) can also be viewed as a linear autoregressive moving average (ARMA) model (without the noise term).37 Note that here we have used RRti instead of μRR(ti) as regressors since this would require a higher order p due to the long-range dependence of μRR(ti) under a very small timescale.
  • Case (b): Dropping off the last two quadratic terms in the Volterra series expansion (3), we obtain
    where <RR>t=1/k=1RRtk(=max{p,r}) denotes the local mean of the past [ell] R-R intervals. Eq. (5) is a bivariate bilinear system,62 which can also be viewed as a nonlinear ARMA or nonlinear ARX model.37

In the current paper, with the goal to study HR responses to BP, Eqs. (4) and (5) will be used for modeling the instantaneous mean μRR(t) of the inverse Gaussian distribution defined in Eq. (1).

Adaptive Point Process Filtering

Let ξ = [{ai}i=0p, {bj}j=1q,{hkl},θ]T denote the vector that contains all unknown parameters in the heartbeat interval probability model. In order to adapt the model in a nonstationary environment, we can recursively estimate the parameters via adaptive point process filtering.5 A state-space formulation of the discrete-time point process filtering algorithm is described here:

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

where P and W denote the parameter and noise covariance matrices, respectively; and Δ denotes the time bin size. The choice of bin size reflects the timescale of estimation interest, we often use Δ=5 ms. Diagonal noise covariance matrix W, which determines the level of parameter fluctuation at the timescale of Δ, can be either initialized empirically from a random-walk theory or estimated from a maximum likelihood estimate. Symbols λk=λkξkand  2λk=2λkξkξkT denote the first- and second-order partial derivatives of the CIF with respect to ξ at time t = kΔ, respectively. The indicator variable nk = 1 if a heart beat occurs in time ((k − 1)Δ,kΔ] and 0 otherwise. The point process filtering equations can be viewed as a point process analog of the Kalman filtering equations (for continuous-valued observations).28 Given a predicted (a priori) estimate ξk|k−1, the innovations [nk − λkΔ] is weighted by Pk|k−1([nabla]log λk) (viewed as an adaptation gain) to further produce the filtered (a posteriori) estimate ξk|k. Since the innovations term is likely to be nonzero in the absence of a beat, the parameters are always updated at each time step k. The a posteriori covariance Pk|k is derived based on a Gaussian approximation of the log-posterior.

Closed-loop Cardiovascular Control

The cardiovascular control on arterial blood pressure can be modeled as a closed-loop system,23, 6, 3 and it can be illustrated by the simplified diagram shown in Fig. 1. In the feedforward pathway (RR→BP), the R-R intervals influence the forthcoming BP measure as defined by the H12 transfer function (either in the time or frequency domains), including the effects of heart contractility and vasculature tone on arterial pressure. In the feedback pathway (BP→RR), the baroreceptors modulate the beat-to-beat interval through a feedback control mechanism, where the H21 transfer function includes both the baroreflex and the autonomic control to the heart. It shall be pointed out that a complete model of the cardiovascular control system can be far more complex. However, in order to focus our attention on the mentioned closed-loop interaction (in the dashed box in Fig. 1), we have purposely ignored all other possible effects, including the important respiratory influence directly affecting the considered loop.

Figure 1
A simplified diagram of closed-loop model of cardiovascular system. The beat-to-beat R-R interval is modulated by the blood pressure (BP) through the feedback baroreflex loop. The dashed box is the closed-loop model we aim to identify. In addition, the ...

Modeling the Baroreflex Gain (BP→RR)

First, we aim to assess the BP→RR feedback pathway, which is directly related to the HR baroreflex. In the literature, baroreflex sensitivity was usually either estimated from a time-domain based sequence method by regression on series sequences,44, 35 or estimated from a frequency-domain based coherence method.21 These two methods are purely data-driven, and therefore might be sensitive to the outliers within the selected samples. Here, we estimate the BRS through a closed-loop parametric bivariate AR model.3, 43 Moreover, our point process model is adaptive and sidesteps the local stationarity assumption, therefore it is capable of capturing the non-stationary nature of the physiological signals due to the drastic cardiovascular control compensatory changes. Specifically, in light of Eq. (4) we can compute the frequency response for the baroreflex (BP→RR)


where f1 and f2 denote the rate (beat/s) for the R-R and BP-BP intervals, respectively; here we assume f1f2 [equivalent] f (namely, the heartbeat period is about the same as the BP-event period). The order of the AR model also determines the number of poles, or oscillations, in the frequency range. Modifying the AR coefficients is equivalent to changing the positions of the poles and reshaping the frequency response curve. With the time-varying AR coefficients {ai(t)} and {bj(t)} estimated from the point-process filter, we may evaluate the dynamic frequency response of Eq. (6) at different ranges (LF, 0.04–0.15 Hz; HF, 0.15-min{0.5,0.5/RR} Hz, where 0.5/RR denotes the Nyquist sampling frequency). The frequency-dependent baroreflex gain, characterized by |H12(f)|, represents the effect of BP on heartbeat, mediated by the neural autonomic reflex. Since the R-R interval is also influenced by respiratory input (Fig. 1) at the HF range, it is more common and meaningful to examine the baroreflex gain at the LF range.

Modeling the Feedforward Gain (RR→BP)

Simultaneous to baroreflex assessment, we aim to model the RR→BP feedforward pathway, which enables us to evaluate the impact of the heartbeat durations on the hemodynamics. Similar to Eq. (4), BP is also modeled by a bivariate linear AR model:


where μRR(ki) represents the estimated instantaneous R-R mean value at the time bin when BP-events occur. The coefficients {ci(k)}i=0p and {di(k)}i=1p are dynamically tracked by a RLS filter (viewed as a special kind of Kalman filter).34 Unlike the point-process filter, the update occurs only at the occurrence time of BP-events. Similarly, we also estimate the frequency response of the RR→BP pathway:


where f denotes the sampling rate (beat/s) for BP-BP intervals. Likewise, we can estimate the dynamic gain and phase of H21(f) at each single BP-event (whereas during between-events period, the coefficient estimates remain unchanged).

Estimating the Dynamic Coherence

Let QBP (f) and QRR(f) denote the power spectra of the BP and R-R series, respectively. To estimate the cross-spectrum in the closed-loop system, we assume that the noise variance and the nonlinear interactions in the feedforward and feedback loops are sufficiently small. Given the BRS, we can estimate the cross-spectrum (between BP and RR) in the feedback loop as Cuy (f) ≈ H12 (f)QBP (f). As the coefficients {ai(t)} and {bj(t)} are iteratively updated in time, the point-process filter produces a direct assessment of instantaneous (parametric) baroreflex gain as well as cross-spectrum at a very fine temporal resolution without using a window-based identification. Similarly, we can estimate the cross-spectrum in the feedforward pathway: Cyu (f) ≈ H21 (f)QRR (f), where QRR (f) can be estimated from


from which we can also compute the time-varying LF/HF power ratio. Note that the unit of the time-varying R-R spectrum in Eq. (9) is cycles/beat, and we have assumed that the variance σBP2(t) (estimated from the feedforward pathway) remains unchanged between two consecutive systolic values. Furthermore, the instantaneous normalized cross-spectrum (i.e., coherence) can be computed as


where |·| denotes the modulus of a complex variable. The second equality holds due to the fact that 𝒞yu(f)=𝒞uy*(f)=H12(f)𝒬BP(f), where * denotes the Hermitian operator (note that | Cyu |=| Cuy | and they are anti-phase against each other). The third equality indicates the fact that the time-varying coherence function can be expressed by the multiplication of two (feedback and feedforward) time-varying transfer functions,65 computed from Eqs. (6) and (8), respectively. Our time-varying closed-loop coherence function can be computed within very fine timescale, several studies have examined its properties (e.g., stability, numerical bound) in detail.50, 65, 66

Estimating the Dynamic Cross-Bispectrum

If μRR(t) is identified by a bilinear system as in Eq. (5), we can estimate the higher-order statistics from the nonlinear system. For simplicity of computation, we further assume that the RR and BP measurements are all Gaussian (such that their own third and higher-order cumulant statistics are zeros), and the cross-bispectrum between input RR (denoted as x) and output μRR (denoted as y) is given by:


where (f1,f2)=k=1ql=1qhklej2kπf1ej2lπf2 denotes the Fourier transform of the 2nd-order kernel coefficients {hkl}.

Proof: see the APPENDIX for details of deriving Eq. (11).

Furthermore, let h(t) denote a vector that contains all of 2nd-order coefficients {hij(t)}; in light of Eq. (10) and the Parseval’s theorem, it follows that |(f1,f2)|2df1df2=|(f1,f2)|2df1df2=i,j|hij|2. In light of this property, we can define an instantaneous index that quantifies the fractional contribution between the cross-spectrum and the cross-bispectrum:


where |h(t)|=klhkl2(t). The spectrum norm defines the area integrated over the frequency range under the spectral density curve. Since the norm units of spectral and bispectral density are the same, their ratio ρ(t) is dimensionless. The “≈” in Eq. (12) is due to the approximation of a jointly Gaussian assumption used in the proof (see APPENDIX)—in our case, u(t) and x(t) are approximately marginally Gaussian. A small value of ρ implies a presence of significant (nonzero) values in {hkl} (namely, indirect evidence of nonlinearity or bilinear interactions), whereas a perfect linear Gaussian model would imply that ρ = 1. Therefore, the ρ-statistic can be used to characterize the linear/nonlinear frequency interactions of two signals. Because ρ(t) is a function of the estimated parameters, this index can be estimated at each moment in time, and is updated at the beat as well as in-between beats.



Before application to real data, we first tested our proposed point-process filter on simulated data, and compared its performance with the standard method. Specifically, we generated the R-R and RP time series of 10 min using a closed-loop bivariate AR model. For simplicity and illustration purposes, we present a linear case where the BRS has two oscillatory frequencies, one at 0.1 Hz (LF range), and the other at 0.3 Hz (HF range). Data were generated to simulate a baroreflex gain at 0.1 Hz, dropping abruptly from 1.75 to 0.96 at 300 seconds, and a gain at 0.3 Hz (HF) always constant at 1. The synthetic RR and BP traces (from one simulation) are shown in Figs. 2a and 2b, respectively.

Figure 2
(a) Simulated trace of the R-R intervals. (b) Simulated trace of the BP signal. (c) The BRS estimate computed from the adaptive point process filter (among one of 20 simulations) at both LF (solid line) and HF (dashed line) range. (d) Comparison of the ...

We conducted 20 Monte Carlo runs (with different realizations) and in each run estimated the feedback pathway BRS using the adaptive point-process filter. A correct model size was used (thereby with no model selection), and the model parameters were estimated based on the first 200 samples using the least-squared method. To compare tracking ability and estimation accuracy, we also applied an adaptive RLS filter using a forgetting factor 0.98.17 The instantaneous BRS estimates (LF and HF) are shown in Fig. 2c (from one Monte Carlo run). For the point-process filter, the mean±SD statistics of the BRS during the first and second time periods are 1.75±0.09 and 0.97±0.05, respectively; whereas for the RLS filter, the mean±SD statistics of the BRS during the first and second time periods are 1.71±0.05 and 0.90±0.01, respectively. Mainly due to the exponential smoothing window, the RLS filter produces a less variable BRS trajectory, as reflected by the lower variance in the estimate. However, the point process algorithm achieves a better precision, as reflected by the mean statistics, despite the small difference (rank-sum test, statistically non-significant). The averaged traces of the BRS estimate from both algorithms (LF only) are shown in Fig. 2d. Clearly, the point-process filter exhibits a faster tracking ability than the RLS filter during the transient period. Furthermore, the point-process filter reaches the steady state around 360 s, which is also faster than the RLS filter (around 415 s).

Experimental Protocol

A total of fifteen healthy volunteer subjects (mean age 24 ± 4), gave written consent to participate in this study approved by the Massachusetts General Hospital (MGH) Department of Anesthesia and Critical Clinical Practices Committee, the MGH Human Research Committee and the MGH General Clinical Research Center. Subjects were evaluated with a detailed review of his/her medical history, physical examination, electrocardiogram, chest X-ray, a urine drug test, hearing test, and for female subjects, a pregnancy test. Any subject whose medical evaluation did not allow him or her to be classified as American Society of Anesthesiologists (ASA) Physical Status I was excluded from the study. Other exclusion criteria included neurological abnormalities, hearing impairment, and use of either prescribed or recreational psychoactive drugs. Intravenous and arterial lines were placed in each subject. Propofol was infused intravenously using a previously validated computer-controlled delivery system running STANPUMP (a computer controlled delivery system59) connected to a Harvard 22 syringe pump (Harvard Apparatus, Holliston, MA), using the well-established pharmacokinetic and pharmacodynamic models.56, 57 In Subject 1, five effect-site target concentrations (0.0, 1.0, 2.0, 3.0 and 4.0 µg/ml) were each maintained for about 15 minutes respectively, where concentration level-0 corresponds to the conscious and wakefulness baseline. In the remainder of subjects, an additional effect-site target concentration of 5.0 µg/ml was administered. Capnography, pulse oximetry, ECG, and arterial blood pressure (ABP) were monitored continuously by an anesthesiologist team throughout the study. Bag-mask ventilation with 30% oxygen was administered as needed in the event of propofol-induced apnea. Because propofol is a potent peripheral vasodilator, phenylephrine was administered intravenously to maintain mean arterial pressure within 20% of the baseline value. ECG and ABP were recorded at a sampling rate of 1 kHz using a PowerLab ML795 data acquisition system (ADInstruments, Inc., Colorado Springs, CO). Electroencephalogram (EEG) data were also recorded, and the EEG analysis was reported elsewhere.51 Four recordings (Subjects #6, 11, 12, 14) were excluded for analysis either because the subjects fell asleep during the experimental behavioral protocol or because of poor quality of the data recordings. In the present study, we only study the epochs of level-0 (baseline) and level-1 (1.0 µg/ml) drug concentration, and focus on the comparison of these two levels (in both transient period and steady state period). This comparison is also aimed to reduce any possible bias caused by the intravenous infusion of Phenylephrine during the administration of higher concentrations of propofol.

In Fig. 3, the raw R-R interval series and the systolic BP (SBP) measures are shown for five selected subjects during baseline and level-1 drug concentration. Note that under induction of propofol anesthesia, there were moments where SBP recordings were corrupted (see Fig. 3, level-1 SBP traces). Such artifacts, caused by periodic collection of blood samples, were excluded from the analysis. A preliminary empirical statistical analysis of the R-R intervals indicated a better fit by the inverse Gaussian distribution (P < 0.05, Kolmogorov-Smirnov test) as compared with a Gaussian distribution, thus validating our probability model used in Eq. (1). The inverse Gaussian distribution bears a similar shape as the Gaussian distribution, but with a longer tail for modeling data outliers. See Fig. 4 for a clear illustration.

Figure 3
The raw R-R intervals and systolic blood pressure (SBP) recordings at baseline and level-1 drug concentration for selected five subjects. Note that during administration of propofol anesthesia, in each epoch, there were moments where SBP recordings were ...
Figure 4
(a) Histogram of the R-R intervals from one recording epoch. (b) Maximum likelihood fitted Gaussian and inverse Gaussian probability density functions (pdfs) for R-R intervals. Note that Gaussian pdf is symmetric, whereas the inverse Gaussian pdf is slight ...

Model Selection and Goodness-of-Fit

In the present experiment, the SBP measure was used for baroreflex evaluation. For the linear model (Eq. 4), we assumed p = q, the bivariate order was fitted from 2 to 8, and the optimal order was chosen according to the Akaike information criterion (AIC). A preliminary analysis on model order concluded that increasing model order beyond 8 did not yield significant improvement in model fit (sometimes results were even worse). For the nonlinear model (Eq. 5), the order r = 2 was chosen empirically to avoid demanding computation burden, and the initial hij was estimated by fitting the residual error via least-squares.7 During online estimation, the stability of the AR model was examined and always assured (i.e., the poles are within unit circle; if the condition is violated, previous estimates are retained), such that the time-varying coherence function is always bounded by 1. Upon estimating the CIF at every time step, the goodness-of-fit of the probabilistic heartbeat model is evaluated by a one-sided Kolmogorov-Smirnov (KS) test, as well as the autocorrelation independence test (for details, see previous references 4, 5, 17, 18).

For all the considered data, our point process model achieves a quite satisfactory performance in model goodness-of-fit. Among a total of 22 epochs, our probability model with a linear predictor is able to provide estimates within 95% confidence intervals in the KS test for 17 epochs, and the average percentage of the time-rescaled points that fall within the 95% confidence bounds is 89% among all of fits for all epochs. This implies that roughly 77% of the beats are correctly predicted by our probability model, which achieved much better results (in terms of KS distance) than any other window-based heartbeat interval model.4 If the criterion is relaxed to 90%, the success percentage is increased to 91%. In contrast, other standard time-varying models produced very poor goodness-of-fit (none of them within 95% confidence bounds), mainly due to either the beat-to-beat updating structure, or to interpolation biases (see previous study4).

Quantitative Assessment

We first examined the tracking performance of the point-process filter during transient dynamics from level 0 to level 1, when the baroreflex gain may undergo a sudden shift. We use Subject 1 to illustrate this point in Fig. 5. As seen from the figure, the baroreflex responses were triggered by injection of the anesthetic drug around 1890 s, which was indeed accompanied by a drop in the baroreflex gain (a rough decrease of about 55%), with the point-process filter quickly capturing the sudden change. In addition to estimating instantaneous (HR, HRV, and baroreflex) statistics, the adaptive filter also tracks the temporal evolution of the probability distribution (beyond first and second order cumulant statistics) of the heartbeat interval. The KS plot and autocorrelation plot (data not shown) both indicate that the goodness-of-fit of the model for these data is very good.

Figure 5
A snapshot of estimated dynamic statistical indices ranging from the baseline to level-1 (Subject 1). The dashed line marks the start of anesthesia at level-1 drug concentration, and the blank region represents the transient period after the baseline ...

For each subject and each epoch, we computed the following instantaneous indices: μHR, σHR, SBP, coherence, ρ, LF/HF, and BRS. For each recording, to compute corresponding time-averaged indices, one epoch of 10 min was selected before propofol administration (level 0), whereas a second epoch of 10 min was selected after drug administration (level 1). The mean±SD statistics for each of the considered indices from all 11 subjects and 22 epochs are reported in Table 1. The level-1 values in bold indicate decreasing trends from the level-0 baseline to the start of administration of propofol anesthesia (level 1). To elucidate the table, we observed a nearly ubiquitous (subject ratio: 10/11) decrease in the SBP value. However, results in HR and HRV did not show a consistent trend across subjects. The mean statistics of the LF/HF index, which is believed to reflect the sympathetic-vagal balance 1, also did not show dramatic or consistent trends. Conversely, we did observe a clear reduction of baroreflex gain, more dominant in LF (subject ratio: 10/11) than in HF (subject ratio: 8/11). These observations are consistent with published results in the literature.12, 60 This may indicate that under propofol anesthesia the parasympathetic activity is attenuated, and that the baroreflex responses are reset to allow HR modulation at a lower BP than during wakefulness. Furthermore, the reason why the reduction in BRS is less significant in HF, and why we have not found significant shifts in the sympatho-vagal balance as quantified by the LF/HF index, might have been additionally influenced by the fact that the oscillatory rhythms associated with this frequency range are highly influenced by changes in the respiratory patterns elicited by the involuntary respiratory mechanism during propofol anesthesia.

Table 1
The mean±SD statistics (averaged over time in each epoch) of instantaneous HR, HRV, SBP, coherence, and ρ-value during level-0 and level-1 drug concentrations across 11 subjects. The statistic in bold font at level-1 indicates a decrease ...

To conduct a group study analysis on BRS, we computed the “relative” BRS ratio for all eleven subjects, by which the BRS at level-1 drug concentration is normalized with respect to the corresponding gain value at the level-0 baseline. The group statistics are shown in Fig. 6a; clearly, the median BRS ratios at the LF and HF ranges are smaller than 1, indicating an overall decrease (highly significant in LF where the 95 percentile box is entirely below 1) of the barorefeflex gain with anesthesia. In addition, we examined the relationship between SBP and relative BRS ratio. It is worth noting that, although there is a decreasing trend in the SBP value from level-0 to level-1 drug concentration (i.e., the relative SBP ratio is smaller than 1), there is no correlation between the relative SBP ratio (level-1/level-0) and the relative BRS (LF) ratio (level-1/level-0), as shown in Fig. 6b.

Figure 6
(a) Box plot of relative BRS ratio (level-1/level-0) values of 11 subjects at the LF and HF ranges. (b) Scatter plot of relative BRS ratio vs. relative SBP ratio (level-1/level-0).

We also investigated the role of nonlinearity as assessed by the bilinear model before and during induction of propofol anesthesia. Specifically, we compared the time-averaged ρ-ratio statistic between the wakefulness baseline and the level-1 propofol concentration. In all but three subjects (Subjects #1, 2 and 9), the ρ value (LF) is slightly greater (statistically non-significant in paired Mann-Whitney test for 11 subjects, P = 0.0582) in the baseline (level-0) condition than in the level-1 anesthesia condition, which suggests that the bilinear interaction between BP and RR in the LF range became more active or strengthened under propofol anesthesia, where the parasympathetic activity may be suppressed or attenuated.31 The phenomenon of increasing nonlinearity in heartbeat interval dynamics in the presence of attenuated parasympathetic activity is also in agreement with our previous observations in another experimental protocol.17, 18 Meanwhile, the mean coherence was slightly decreased or remained unchanged after the induction of propofol anesthesia in the majority of the subjects (Table 1). A low value of coherence also suggests that the BP-RR relation may have either nonlinear components, or that the two signals are less linearly correlated in the frequency domain. A statistical summary is provided by the box plots of coherence and ρ statistics from all 11 subjects, as shown in Fig. 7. Here the group median ρ-values (level-0 to level-1, LF: 0.56→0.43; HF: 0.35→0.34) point at a mild increase of nonlinear dynamics, accompanied by a slight decrease in linear coherence (LF: 0.83→0.78; HF: 0.75→0.70), but a statistical (unpaired Mann-Whitney) test on the population level did not reveal significant differences. Another noteworthy point from Fig. 7 is that our results indicated that there is more nonlinear interaction in HF (lower median ρ-ratio as well as lower median coherence) than in the LF range, which is somewhat in disagreement with previous findings pointing at more prominent HRV nonlinearities within the VLF and LF ranges. This could be because the current model has not considered the influence of the respiratory input (which is acting within the HF range) in the feedback pathway (Fig. 1), with consequent estimation bias in the HF range, or more simply because the linear coupling is higher in LF, thus reducing the role of nonlinearity as percentage of the total variability; furthermore, the presence of nonlinear effects of RSA, mainly in the HF range, may enhance nonlinear RR-BP dynamics at these frequencies, thus explaining our observations. However, further analysis is required to examine this issue with more controlled data and the inclusion of respiratory measures, so to provide clearer interpretations and give meaning to the pronounced nonlinear dynamics beyond just their association with a higher degree of decoupling between the variables of interest. For example, measuring and analyzing additional variables directly related to vascular tone, resistance dynamics, or blood flow autoregulation, where nonlinear effects play an important role, could be included in our nonlinear model and would help shed more light into such mechanisms, as well as provide important insights into more complex cardiovascular control dynamics.

Figure 7
Box plots of coherence and ρ statistics between level-0 and level-1 anesthetic drug concentration levels.

Further insights on this final issue may be gained by looking at the time-varying coherence trends within the level-0 and level-1 epochs considered for one subject. In Fig. 8 the instantaneous coherence between SBP and RR shows that these two series are strongly correlated at the HF range, first staying around 0.3 Hz along the wakefulness baseline and then shifting around 0.25 Hz at level-1 of propofol concentration. This is likely due to a slowdown of the respiratory rhythm, which influences both SBP and RR (see Fig. 1). In the LF range, as generally lower and less consistent values of peak coherence are observed, it is important to focus the attention on the sudden drop in coherence at 500–600 ms since the start of propofol administration. Allegedly, the instantaneous identification indicates the exact moment in time where the autonomic-hemodynamic linear coupling generating the Meyer waves may be completely lost, and SBP and RR variability may be attributed solely to the respiratory influence (in HF).

Figure 8
For Subject 1, the comparison of the linear coherence between RR and SBP series during the wakefulness (level-0) baseline (top two panels) and during the level-1 of anesthetic drug concentration (bottom two panels). The coherence amplitude is color coded, ...


In cardiovascular signal processing, growing efforts have been proposed in devising transfer function models for analyzing cardiovascular regulation.54, 46, 3 The assumption for using transfer function analysis for estimating baroreflex gain depends on the linear relationship between RR and BP through the autonomic feedback gain to the heart.7, 50, 47 From a more mathematical point of view, using a closed-loop bivariate AR model and through system identification, the bivariate AR coefficients are able to provide a compact description of the autospectrum and cross-spectrum between the RR and BP series.7, 50, 66

Although the adaptive point-process filter overcomes critical issues related to the inherent non-stationarity of the data, a sensible initialization of model parameters is still very important for determining the correct estimate for baroreceptor sensitivity. Currently, our initial model parameters were identified by least-squares estimation using short segments of recordings during each epoch. However, if the steady state has not been reached or there is a strong trend shift during the epoch, it is likely that the initialized parameters are insufficient or inaccurate to characterize the remaining data recordings. A practical way to alleviate this issue is to estimate the parameters from time to time or based on other observed physiological markers.

In general, the baroreflex LF gain is thought to be more relevant to HR control, which is responsible for translating specific BP oscillations (the Meyer waves) into automatic outflow modulation to the SA node. Nevertheless, the baroreflex HF gain is believed to be of vagal origin and has sufficiently high coherence during positive pressure ventilation.32 In our current model, the baroreflex gain is estimated through the frequency response of a linear bivariate system. Therefore, we have implicitly assumed that baroreflex sensitivity, the relation of HR responses to changing BP, is purely linear. In a complex physiological system, this assumption is certainly a subject of debate and might deviate from the truth of physiology.25 Indeed, the baroreflex responses normally exhibit a sigmoid or S-shape curve.23 However, we believe that the response curve is approximately linear in the region around the BP set point. Shall BP oscillate in ranges far from set point, our model would likely yield a biased estimate of the baroreceptor gain. Nevertheless, at the same time the nonlinear structures incorporated into the model would allow for quantification of specific markers (e.g., the ρ index), which point at the present of relevant nonlinear dynamics and invite a more cautious interpretation of the linear assessment.

Another limitation of our current probability model is that we have not considered the influence of respiratory input on the R-R intervals and hemodynamics, while it has been well known that respiration also modulates HR through the RSA mechanism;17 consequently, the change of respiratory patterns during propofol anesthesia is likely to bias the baroreflex gain estimate (specially at the HF range). This simplification was purposely chosen in order to make our mathematical analysis clear and simpler, but the possibility of inclusion of other physiological covariates into the model is attainable and will be explored in future studies.


In cardiology, as BRS is known to be variable and non-stationary at various physiological or pharmacological conditions,27 it would be informative to monitor it continuously as a marker of autonomic-mediated hemodynamic control. Modeling and monitoring non-stationary dynamics of physiological systems continues to be an active research topic.4, 5, 65, 66. In this paper, we have proposed a point process method for assessing the baroreflex control of heart rate before and during induction of propofol anesthesia, as applied to 11 healthy subjects. The investigation focuses on the comparison of BRS between the periods during awake baseline and after an initial concentration of drug is administered. In a non-stationary scenario where the physiological state may change dramatically, dynamic assessment of baroreflex control during the transient period is of vital importance. The study of the transient periods due to physiological or pharmacological changes in cardiovascular system has demonstrated the capacity of the point-process filter to quickly capture fast physiologic changes, when baroreflex responses are supposedly triggered, and consequently accompanied by a significant drop in the instantaneous baroreflex gain. Our method potentially provides a real-time indicator for monitoring the baroreflex control of HR during induction of general anesthesia. The analysis of steady-state periods showed a clear reduction of baroreflex gain in the LF region from baseline to the start of anesthesia, accompanied by a decrease in SBP. This might suggest that the baroreflex responses were reset during propofol anesthesia to control HR at a lower BP than during wakefulness and that the quantitative baroreflex gain decreased after administration of propofol anesthesia. The shift in the HR/BP set point may reflect propofol’s systemic vasodilatory effect, whereas the reduction in baroreflex gain is most likely the result of disruption of cardiac control within the central nervous system.

In a closed-loop cardiovascular control system, interactions between the heartbeat interval (or HR) and BP dynamics can be complex and state-dependent. We have proposed novel ways to assess their interactions at arbitrarily fine timescales, as well as in the frequency domain. The statistical study of the linear vs. bilinear cross-spectrum ratio shows greater nonlinear dynamics within the LF range in the awake baseline condition than in the level-1 propofol concentration when the parasympathetic activity starts to attenuate, suggesting that the bilinear interactions between the two variables become more active under propofol anesthesia.

In summary, we have developed a point process model for assessing the baroreflex control of heart rate during induction of propofol anesthesia using clinical recordings. The probability point process model used here can be viewed as one specific example within a more general framework.13, 16 The proposed point process method also enables us to estimate instantaneous HR, HRV, coherence and cross-bispectrum. All of these statistical indices may serve as potential indicators for ambulatory monitoring in clinical practice, and may particularly provide a valuable quantitative assessment of the interaction between heartbeat dynamics and hemodynamics during general anesthesia. More importantly, these quantitative indices could be monitored intraoperatively in order to improve drug administration and reduce the side-effects of anesthetic drugs.


The research was supported by NIH Grants R01-HL084502 (R.B.), K25-NS05758 (P.L.P.), DP2-OD006454 (P.L.P.), T32NS048005 (G.H.), DP1-OD003646 (E.N.B.), and R01-DA015644 (E.N.B.), as well as a CRC UL1 grant RR025758 (P.L.P.). The authors thank L. Citi, K. Habeeb, R. Merhar, A. Salazar and C. Tavares for assistance in collecting and preprocessing the data used in our experiments. We also thank the valuable comments from three reviewers that help to improve the manuscript. Preliminary results of this work have been reported in Proceedings of IEEE ICASSP’09, Taiwan.


For clarity of proof of Eq. (11), we assume that inputs u(t) and x(t) are both stationary and have zero means (in our case, we only model the “de-meaned” RR and BP signals with {ai} and {bj}). The line of logic is similar to the one that was presented in the literature,42 which only considered a univariate input. We first decompose the output y(t) into three (two linear and one bilinear) terms


From the assumption that E[x(t)] = 0 and E[u(t)] = 0, it further follows that


where the expectation operation is averaged on the argument over time.

Next, we compute the cross third-order cumulant statistic between x(t) and y(t) (viz. cross bicovariance):


where we have used the fact that if x(t) is zero-mean Gaussian, then the third cumulant statistic Cxxx12) is zero, such that iai𝔼{x(t+τ1)x(t+τ2)x(ti)}=0. Furthermore, if u(t) and x(t) are zero-mean jointly Gaussian distributed, and the odd-moment statistic Cxxu12) ≈ 0 (because of the symmetry of Gaussian distribution), then the following relationship holds:42, 55




In light of Eqs. (A.2) and (A.3), we obtain

Cxxy(τ1,τ2)=2klh(kτ1,lτ2)Cxx(k)Cxu(l)  =2(2π)2(f1,f2)ej2πτ1f1ej2πτ2f2𝒬RR(f1)𝒞xu(f2)df1df2

Finally, we compute the two-dimensional Fourier transform of Cxxy12) to obtain the cross-bispectrum Cxxy(f1,f2):


which then completes the proof of Eq. (11).


1. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. Heart rate variability. Circulation. 1996;93:1043–1065. [PubMed]
2. Barbieri R, Waldmann RA, Di Virgilio V, Triedman JK, Bianchi AM, Cerutti S, Saul JP. Continuous quantification of baroreflex and respiratory control of heart rate by use of bivarate autoregressive techniques. Ann. Noninvasive Electrocardiology. 1996;3:264–277.
3. Barbieri R, Parati G, Saul JP. Closed- versus open-loop assessment of heart rate baroreflex. IEEE Eng. Med. Biol. 2001;20:33–42. [PubMed]
4. Barbieri R, Matten EC, Alabi AA, Brown EN. A point-process model of human heartbeat intervals: new definitions of heart rate and heart rate variability. Am J. Physiol. Heart Cicr. Physiol. 2005;288:H424–H435. [PubMed]
5. Barbieri R, Brown EN. Analysis of heart beat dynamics by point process adaptive filtering. IEEE Trans. Biomed. Eng. 2006;53:4–12. [PubMed]
6. Baselli G, Cerutti S, Civardi S, Malliani A, Pagani M. Cardiovascular variability signals: towards the identification of a closed-loop model of the neural control mechanisms. IEEE Trans. Biomed. Eng. 1988;35:1033–1046. [PubMed]
7. Baselli G, Porta M, Rimoldi O, Pagani M, Cerutti S. Spectral decomposition in multichannel recordings based on multivariate parametric identification. IEEE Trans. Biomed. Eng. 1997;44:1092–1101. [PubMed]
8. Betzel J, Mukkamala R, Baselli G, Chon KH. Modeling and disentangling physiological mechanisms: Linear and nonlinear identification techniques for analysis of cardiovascular regulation. Philo. Trans. R. Soc. A. 2009;367:1377–1391. [PMC free article] [PubMed]
9. Bristow JD, Prys-Roberts C, Fisher A, Pickering TG, Sleight P. Effects of anesthesia on baroreflex control of heart rate in man. Anesthesiology. 1969;31:422–428. [PubMed]
10. Brown EN, Barbieri R, Eden UT, Frank LM. Likelihood methods for neural data analysis. In: Feng J, editor. Computational Neuroscience: A Comprehensive Approach. CRC Press; 2003. pp. 253–286.
11. Carlson JT, Hedner JA, Sellgren J, Elam M, Wallin BG. Depressed baroreflex sensitivity in patients with obstructive sleep apnea. Am. J. Respir. Crit. Care Med. 1996;154:1490–1496. [PubMed]
12. Carter JA, Clarke TNS, Prys-Roberts C, Spelina KR. Restoration of baroreflex control of heart rate during recovery from anaesthesia. Br. J. Anaesth. 1986;58:415–421. [PubMed]
13. Chen Z, Brown EN, Barbieri R. A study of probabilistic models for characterizing human heart beat dynamics in autonomic blockade control. Proc. IEEE ICASSP. 2008:481–484. [PMC free article] [PubMed]
14. Chen Z, Brown EN, Barbieri R. A point process approach to assess dynamic baroreflex gain. Proc. Computers in Cardiology. 2008:805–808. [PMC free article] [PubMed]
15. 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. IEEE ICASSP. 2009:333–336. [PMC free article] [PubMed]
16. Chen Z, Brown EN, Barbieri R. A unified point process framework for assessing heartbeat dynamics and cardiovascular control; Proc. IEEE 35th Northeast Bioengineering Conference; 2009. pp. 1–2.
17. 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:1791–1802. [PMC free article] [PubMed]
18. Chen Z, Brown EN, Barbieri R. Characterizing nonlinear heartbeat dynamics within a point process framework. IEEE Trans. Biomed. Eng. 2010;57:1335–1347. [PMC free article] [PubMed]
19. Chon KH, Mullen TJ, Cohen RJ. A dual-input nonlinear system analysis of autonomic modulation of heart rate. IEEE Trans. Biomed. Eng. 1995;43:530–540. [PubMed]
20. Chung OY, Bruehl S, Diedrich L, Diedrich A, Chont M, Robertson D. Baroreflex sensitivity associated hypoalgesia in healthy states is altered by chronic pain. Pain. 2008;138:87–97. [PubMed]
21. Clayton RH, Bowman AJ, Murray A. Measurement of baroreflex gain from heart rate and blood pressure spectra. Physiol. Meas. 1995;16:131–139. [PubMed]
22. Cullen PM, Turtle M, Prys-Roberts C, Way WL, Dye J. Effect of propofol anesthesia on baroreflex activity in humans. Anesth. Analg. 1987;66:115–120. [PubMed]
23. De Boer RW, Karemaker JM, Strackee J. Relationships between short-term blood-pressure fluctuations and heart-rate variability in resting subjects: a spectral analysis approach. Med. Biol. Eng. Comput. 1985;23:352–358. [PubMed]
24. De Ceccoa M, Angrillib A. Measurement of human baroreceptor reflex sensitivity by means of parametric identification. Measurement. 1998;24:187–196.
25. Eckberg DL. Nonlinearities of the human carotid baroreceptor-cardiac reflex. Circ. Res. 1980;47:208–216. [PubMed]
26. Eckberg DL, Harkins SW, Fritsch JM, Musgrave GE, Gardner DF. Baroreflex control of plasma norepinephrine and heart period in healthy subjects and diabetic patients. J. Clin. Invest. 1986;78:366–374. [PMC free article] [PubMed]
27. Eckberg DL. Arterial baroreflexes and cardiovascular modeling. Cardiovasc. Eng. 2008;8:5–13. [PubMed]
28. Eden UT, Frank LM, Solo V, Brown EN. Dynamic analyses of neural encoding by point process adaptive filtering. Neural Comput. 2004;16:971–998. [PubMed]
29. Faes L, Porta A, Cucino R, Cerutti S, Antolini R, Nollo G. Causal transfer function analysis to describe closed loop interactions between cardiovascular and cardiorespiratory variability signals. Biol. Cybern. 2004;90:390–399. [PubMed]
30. Faes L, Nollo G, Chon KH. Assessment of Granger causality by nonlinear model identification: application to short-term cardiovascular variability. Ann. Biomed. Eng. 2007;36:381–395. [PubMed]
31. Feld J, Hoffman W, Paisansathan C, Park H, Ananda RC. Autonomic activity during dexmedetomidine or fentanyl infusion with desflurane anesthesia. J. Clinical Anesthesia. 2003;19:30–36. [PubMed]
32. Fietze I, Romberg D, Glos M, Endres S, Theres H, Witt C, Somers VK. Effects of positive-pressure ventilation on the spontaneous baroreflex in healthy subjects. J. Appl. Physiol. 2004;96:1155–1160. [PubMed]
33. Guyton AC. Textbook of Medical Physiology. 8th ed. Philadelphia, PA: Harcourt Brace; 1991.
34. Haykin S. Adaptive Filter Theory. 4th ed. NJ: Prentice Hall; 2001.
35. Hughson RL, Quintin L, Annat G, Yamamoto Y, Gharib C. Spontaneous baroreflex by sequence and power spectral methods in humans. Clin. Physiol. 1993;13:663–676. [PubMed]
36. Jo JA, Blasi A, Valladares EM, Juarez R, Baydur A, Khoo MCK. A nonlinear model of cardiac autonomic control in obstructive sleep apnea syndrome. Ann. Biomed. Eng. 2007;35:1425–1443. [PubMed]
37. Lu S, Ju KH, Chon KH. A new algorithm for linear and nonlinear ARMA model parameter estimation using afne geometry. IEEE Trans. Biomed. Eng. 2001;48(10):1116–1124. [PubMed]
38. Mainardi LT, Bianchi AM, Furlan R, Piazza S, Barbieri R, de Virgilio V, Malliani A, Cerutti S. Multivariate time-variant identification of cardiovascular variability signals: a beat-to-beat spectral parameter estimation in vasovagal syncope. IEEE Trans. Biomed. Eng. 1997;44(10):978–988. [PubMed]
39. Mainardi LT. On the quantification of heart rate variability spectral parameters using time-frequency and time-varying methods. Phil. Trans. R. Soc. A. 2009;367:255–275. [PubMed]
40. Marmarelis VZ. Nonlinear Dynamic Modeling of Physiological Systems. New York: New York; 2004.
41. Nagasaki G, Tanaka M, Nishikawa T. The recovery profile of baroreflex control of heart rate after isoflurane or sevoflurane anesthesia in humans. Anesth. Analg. 2001;93:1127–1131. [PubMed]
42. Nikias C, Petropulu AP. Higher Order Spectra Analysis: A Non-Linear Signal Processing Framework. Prentice Hall; 1993.
43. Nollo G, Porta A, Faes L, Del Greco M, Disertori M, Ravelli F. Causal linear parametric model for baroreflex gain assessment in patients with recent myocardial infarction. Am. J. Physiol. Heart Circ. Physiol. 2001;280:H1830–H1839. [PubMed]
44. Parati G, DiRienzo M, Mancia G. Dynamic modulation of baroreflex sensitivity in health and disease. Ann. NY Acad. Sci. 2001;940:469–487. [PubMed]
45. Peng C-K, Havlin S, Stanley HE, Goldberger AL. Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos. 1995;5:82–87. [PubMed]
46. Pinna GD, Maestri R. New criteria for estimating baroreflex sensitivity using the transfer function method. J. Med. Biol. Eng. Comput. 2001;40:79–84. [PubMed]
47. Pinna GD. Assessing baroreflex sensitivity by the transfer function method: what are we really measuring? J. Appl. Physiol. 2007;102:1310–1311. [PubMed]
48. Poon C-S, Merrill CK. Decrease of cardiac chaos in congestive heart failure. Nature. 1997;389:492–495. [PubMed]
49. Porta A, Aletti F, Vallais F, Baselli G. Multimodal signal processing for the analysis of cardiovascular variability. Phil. Trans. R. Soc. A. 2009;367:391–409. [PubMed]
50. Porta A, Furlan R, Rimoldi O, Pagani M, Malliani A, van de Borne P. Quantifying the strength of linear causal coupling in closed loop interacting cardiovascular variability signals. Biol. Cybern. 2002;86:241–251. [PubMed]
51. Purdon PL, Pierce ET, Bonmassar G, Walsh J, Harrell G, Kwo J, Deschler D, Barlow M, Merhar RC, Lamus C, Mullaly CM, Sullivan M, Maginnis S, Skoniecki D, Higgins H, Brown EN. Simultaneous electroencephalography and functional magnetic resonance imaging of general anesthesia. Ann. NY Acad. Sci. 2009;1157:61–70. [PMC free article] [PubMed]
52. Peden CJ, Cloote AH, Stratford N, Prys-Roberts C. The effect of intravenous dexmedetomidine premedication on the dose requirement of propofol to induce loss of consciousness in patients receiving alfentanil. Anaesthesia. 2001;56:408–413. [PubMed]
53. Sato M, Tanaka M, Umehara S, Nishikawa T. Baroreflex control of heart rate during and after propofol infusion in humans. Br. J. Anaesth. 2005;94:577–581. [PubMed]
54. Saul JP, Berger RD, Albrecht P, Stein SP, Chen MH, Cohen RJ. Transfer function analysis of the circulation: unique insights into cardiovascular regulation. Am. J. Physiol. Heart. Circ. Physiol. 1991;261:H1231–H1245. [PubMed]
55. Schetzen M. The Volterra and Wiener Theories of Nonlinear Systems. New York: Wiley; 1980.
56. Schnider TW, Minto CF, Gambus PL, Andresen C, Goodale DB, Shafer SL, Youngs EJ. The influence of method of administration and covariates on the pharmacokinetics of propofol in adult volunteers. Anesthesiology. 1998;88:1170–1182. [PubMed]
57. Schnider TW, Minto CF, Shafer SL, Gambus PL, Andresen C, Goodale DB, Youngs EJ. The influence of age on propofol pharmacodynamics. Anesthesiology. 1999;89:67–72. [PubMed]
58. Sellgren J, Ejnell H, Elam M, Pontén J, Wallin BG. Sympathetic muscle nerve activity, peripheral blood flows, and baroreceptor reflexes in humans during propofol anesthesia and surgery. Anesthesiology. 1994;80:534–544. [PubMed]
59. Shafer A, Doze VA, Shafer SL, White PF. Pharmacokinetics and pharmacodynamics of propofol infusions during general anesthesia. Anesthesiology. 1988;69:348–356. [PubMed]
60. Tanaka M, Nagaski G, Nishikawa T. Moderate hypothermia depresses arterial baroreflex control of heart rate during, and delays it recovery after, general anesthesia in humans. Anesthesiol. 2001;95:51–55. [PubMed]
61. Tanaka M, Nishikawa T. The concentration-dependent effects of general anesthesia on spontaneous baroreflex indices and their correlations with pharmacological gains. Anesth. Analg. 2005;100:1325–1332. [PubMed]
62. Tsoulkas V, Koukoulas P, Kalouptsidis N. Identification of input output bilinear systems using cumulants. IEEE Trans. Signal Process. 2001;49:2753–2761.
63. Xiao X, Mullen TJ, Mukkamala R. System identication: a multi-signal approach for probing neural cardiovascular regulation. Physiol. Meas. 2005;26:R41–R71. [PubMed]
64. Wang H, Ju K, Chon KH. Closed-loop nonlinear system identification via the vector optimal parameter search algorithm: Application to heart rate baroreflex control. Med. Eng. & Phy. 2007;29:505–515. [PubMed]
65. Zhao H, Lu S, Zou R, Ju K, Chon KH. Estimation of time-varying coherence function using time-varying transfer functions. Ann. Biomed. Eng. 2005;33:1582–1594. [PubMed]
66. Zhao H, Cupples WA, Ju K, Chon KH. Time-varying causal coherence function and its application to renal blood pressure and blood flow data. IEEE Trans. Biomed. Eng. 2007;54:2142–2150. [PubMed]
67. Zou R, Wang H, Chon KH. A robust time-varying identification algorithm using basis functions. Ann. Biomed. Eng. 2003;31:840–853. [PubMed]