Simulation

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 , respectively.

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 (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 . 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 system^{59}) 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 , 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 , 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 for a clear illustration.

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

*h*_{ij} 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 study^{4}).

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 . 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.

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 . 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 (more ...) |

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 ; 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 .

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 (). 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 . 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 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 (), 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.

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 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 ). 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).

Discussion

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.