PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Biomed Eng. Author manuscript; available in PMC 2010 July 1.
Published in final edited form as:
PMCID: PMC2804879
NIHMSID: NIHMS130493

Assessment of Autonomic Control and Respiratory Sinus Arrhythmia Using Point Process Models of Human Heart Beat Dynamics

Zhe Chen, Member, IEEE, Emery N. Brown, Fellow, IEEE, and Riccardo Barbieri, Senior Member, IEEE

Abstract

Tracking the autonomic control and respiratory sinus arrhythmia (RSA) from electrocardiogram and respiratory measurements is an important problem in cardiovascular control. We propose a point process adaptive filter algorithm based on an inverse Gaussian model to track heart beat intervals that incorporates respiratory measurements as a covariate and provides an analytic form for computing a dynamic estimate of RSA gain. We use Kolmogorov-Smirnov tests and autocorrelation function analyses to assess model goodness-of-fit. We illustrate the properties of the new dynamic estimate of RSA in the analysis of simulated heart beat data and actual heart beat data recorded from subjects in a four-state postural study of heart beat dynamics: control, sympathetic blockade, parasympathetic blockade, and combined sympathetic and parasympathetic blockade. In addition to giving an accurate description of the heart beat data, our adaptive filter algorithm confirms established findings pointing at a vagally mediated RSA, and it provides a new dynamic RSA estimate that can be used to track cardiovascular control between and within a broad range of postural, pharmacological and age conditions. Our paradigm suggests a possible framework for designing a device for ambulatory monitoring and assessment of autonomic control in both laboratory research and clinical practice.

Index Terms: Adaptive filters, autoregressive (AR) processes, heart rate variability (HRV), point processes, respiratory sinus arrhythmia (RSA)

I. Introduction

Heart rate (HR) and heart rate variability (HRV) are important quantitative markers of cardiovascular control, as regulated by the autonomic nervous system [1]. It has long been understood that the healthy heart is influenced by multiple neural and hormonal inputs that result in variations of duration in the interbeat intervals (R-R intervals). The synergic interaction between the two branches of the autonomic nervous system to the heart has a major influence in regulating the cardiac dynamics and physiological mechanism of HRV. In particular, parasympathetic influences decrease the firing rate of pacemaker cells in the heart's sinus-atrial (SA) node, whereas sympathetic influences have the opposite effect [9]. In cardiovascular physiology, it is known that lung volume tends to be correlated with variations in the timing of heart beat, or HRV. Typically, HR slows down during expiration and speeds up during inspiration. This phenomenon is known as the respiratory sinus arrhythmia (RSA) [28]. RSA is primarily mediated by modulation of vagal outflow to the SA node. Quantification of RSA provides important information about some of the mechanisms involved in cardio-respiratory coupling [20], [28]. In clinical practice, RSA is often treated as an indirect and noninvasive measure of parasympathetic cardiac control [24], [20], even in the presence of paced breathing [29], and may also be considered as a reliable indicator of cardiac dysfunction [26]. A quite comprehensive review on RSA may be found in [10].

A central goal in biomedical engineering applied to cardiovascular control is to develop quantitative measures and informative indices that can be extracted from physiological measurements. Specifically, a major challenge in cardiovascular engineering is to develop statistical models and apply signal processing tools to investigate various cardiovascular-respiratory functions [8], such as HRV, RSA, and baroreflex.

In the literature, numerous methods have been proposed for quantitative HRV analysis [1], [25], including point process analysis [2], [3], frequency-domain analysis [7], and nonlinear dynamics analysis [19]. In [22], RSA was defined as “the difference between the maximum HR rate after the onset of inspiratory flow and the immediately minimum HR"; whereas in [35], RSA was calculated using the formula: 100 x (mean longest R-R - mean shortest R-R) / mean R-R interval. Saul and colleagues [27] proposed a transfer function analysis approach for evaluating the RSA, which requires to directly model the SA node. In [4], [5], a bivariate autoregressive model was proposed to evaluate a time-varying index of RSA (within a temporal window). However, none of these RSA indices provide a truly instantaneous evaluation of the cardiorespiratory dynamics.

Several issues in RSA assessment from R-wave events have yet to be addressed [10]. First, estimates of RSA have been derived from either HR or heart period data. The former is more commonly computed in clinical practice, whereas the latter would be preferred on biometric grounds, especially when the interest is in indexing parasympathetic control because of the relative linearity between vagal frequency and heart period. Second, the R-R interval series are unevenly spaced in time. Direct application to these data to spectral analysis is not appropriate and is usually solved by use of interpolating filters. In addition, longer heart periods may significantly decrease the Nyquist frequency under fast respiratory oscillation, giving rise to possible aliasing effects. Third, standard time-series analysis usually assumes that the data show at least weak sense stationarity, thus requiring particular care in choosing appropriate data segments for analysis, or requiring removal of nonstationary trends.

To address these issues, we investigate different probability models for human heart beat intervals with an adaptive point process filtering paradigm [3], and illustrate the analysis with both synthetic data and electrocardiogram (ECG) and lung volume recordings from a previous study [31], [32], [33] under an autonomic blockade assessment protocol. Furthermore, we extend the inverse Gaussian probability model to take into account the influence of respiration on HRV, based on which we derive an analytic form for computing the frequency response of RSA. Modeling accuracy is evaluated via goodness-of-fit tests, including the Kolmogorov-Smirnov (KS) test. To our knowledge, this is the first effort in the literature that investigates the dynamic RSA effect within a point process framework. In line with the most conventional guidelines [1], [10], our paradigm resolves the old conflict between heart period and HR, obviates the need of interpolation with the potential to solve possible aliasing problems, allows for instantaneous measures at virtually any time resolution offering (in contrast to the interpolated R-R interval values) a more rigorous frequency analysis, and overcomes nonstationarity issues associated with window-based estimation models.

The manuscript is organized as follows. Section II first introduces the point process framework, then presents several probability models of the heart beat dynamics, and finally proposes the extended bivariate probability model by inclusion of respiration (e.g., lung volume) measurements for the purpose of quantifying instantaneous RSA. In Section III, both synthetic and real experimental recordings are used to illustrate and validate the instantaneous RSA gain as computed by the novel point process algorithm. In addition, statistical tests are conducted on the autonomic blockade protocol to evaluate inter-subject statistical trends of RSA gain across different posture and pharmacological conditions. Finally, discussions and conclusion are given in Section IV.

II. Methods

In this section, we present the heart beat interval and the HR probability models, the instantaneous estimates of heart rate and HRV from the heart beat interval (RR) model parameters, the extension to a bivariate model with respiration (RP) as covariate to derive RSA measures, the point process adaptive filtering algorithm for instantaneous assessment of the HR and RSA indices, and finally goodness-of-fit tests to evaluate how well these estimates describe the stochastic structure of the wave events extracted from an ECG. A detailed description of the conceptual framework of the history-dependent point process model and the adaptive paradigm can be found in our previous publications [2], [3].

A. Point Process Probabilistic Models of the Heart Beat Interval

The R-wave events mark the electrical impulses from the heart’s conduction system that represent ventricular contractions. Hence, they are a sequence of discrete occurrences in continuous time, and as such, form a point process. Suppose that we are given a set of discrete R-wave events {uj}j=1J detected from the ECG, let RRj = uj − uj−1 > 0 denote the jth R-R interval, or equivalently, the waiting time until the next R-wave event. By treating the R-waves as discrete events, we propose different parametric point process probability models (Table I) in the continuous-time domain.

TABLE I
Comparison of 4 two-parameter probability models for the heart beat R-R interval and the heart rate (HR). Here τ = tuj is treated as the waiting time random variable; θ12 are two parameters that characterize the respective ...

As an example, assuming history dependence, the waiting time until t − ujthe next R-wave event may be modeled as an inverse Gaussian model as follows [18]:

p(t)=(θ2πt3)12exp(θ(tujμt)22μj2(tut))(t>uj)
(1)

where uj denotes the previous R-wave event occurred before time t, θ > 0 denotes the shape parameter, and μt denotes the instantaneous R-R mean value. Because the parasympathetic and sympathetic inputs to the SA node can occur on a millisecond timescale, but their effects can last for several seconds, the intervals must be modeled as dependent on the recent history of the SA node inputs:

μtμRR(t)=a0+i=1paiRRti
(2)

Namely, the mean value is modeled by a univariate p-order autoregressive (AR) process, which is assumed (approximately) to be influenced by the past p R-R interval values. In order to track the nonstationary behavior of heart beat dynamics, all of parameters involving in (1) and (2) are both adaptive, hence the instantaneous mean μRR (t) is time-varying, which is determined by the time-varying AR coefficients {ai(t)}i=0p. On the other hand, the instantaneous variance of the inverse Gaussian model can be derived as [18]

σRR2(t)=μRR3(t)θ(t),
(3)

In a similar fashion, we can derive the mean and variance of R-R interval for other probability models, such as the Gaussian, lognormal, and gamma models (see Table I for a brief summary).

B. Instantaneous Indices of HR and HRV

Heart rate (HR) is defined as the reciprocal of the R-R interval. For t measured in seconds, a new variable r = c(t − ut)−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(t − ut)−1 is given by

p(r)=|dtdr|p(t),
(4)

and the mean and the standard deviation of HR r can be derived (see Table I). Essentially, the instantaneous indices of HR and HRV are characterized by the mean μHR and standard deviation, σHR respectively. In the case of inverse Gaussian model, we have [2], [3]:

μHR=μ˜1+θ˜1
(5)

σHR=(2μ˜+θ˜)/μ˜θ˜2
(6)

where [mu]=c−1μRR and [theta w/ tilde]=c−1θ.

C. A Bivariate Probabilistic Model for RSA Assessment

In the probability models considered thus far, we have only used the R-R interval time series to estimate the instantaneous mean μRR. Physiology suggests that HR is influenced by other physiological covariates, such as changes in lung volume due to respiratory activity [4], [5]. This fact further motivates us to incorporate respiration as a covariate into the model. Specifically, for the inverse Gaussian model, we may replace the instantaneous mean in (1) by

μtμRR(t)=a0+i=1paiRRti+j=1qbjRPtj
(7)

where RPt−j denotes the previous jth respiration measurement before time t. Now, the instantaneous mean μRR is described by a bivariate AR-type model. It should also be noted that the extended bivariate model shall not be limited to the inverse Gaussian distribution; it is also straightforward to extend it to the Gaussian or lognormal distribution, depending on whichever is more desirable.

Note that in general, the measurements of RR (beat/cycle) and RP have different sampling frequencies. In practice, our framework allows for two ways to tackle this issue. The first approach is to resample the RP measurements to synchronize with the heart beat and obtain the respiration values at the beat time; the second one is to treat them as separate measurements (with different sampling rates) and to conduct frequency analysis with extra caution. The second approach could be useful to avoid aliasing effects in the case of low Nyquist frequency in the presence of long heart beats.

In terms of frequency analysis, the RSA effect is reflected by the fact that the R-R interval has a spectral component modulated by the respiratory variable [27]. With a linear system assumption, RSA can be estimated with transfer function or frequency response analysis using standard signal processing tools [27], [5].

Given the parametric AR model (7), we can evaluate the frequency response for the R-R interval itself

H1(ω)=11i=1pai(k)zi|z=ej2πfs,
(8)

as well as the frequency response for the RSA

H12(ω)=j=1qbj(k)zi|z=ej2πfs1i=1pai(k)zi|z=ej2πfs,
(9)

where fs is the beat rate of the RR, and RP are sampled at the same frequency as the beat. [1Here we assume that RR and RP measurements have the same sampling rate in (7). In practice, the RP time series typically has a higher sampling rate, but it can always be resampled and interpolated to obtain data points at the time of heart beats.] With the estimated time-varying AR coefficients {ai(k)} and {bj(k)} at time t = kΔ, we may evaluate the dynamic power spectrum (parametric autospectrum) or the gain (amplitude) in the frequency domain [7]:

PRR(ω,t)=σRR(t)|H1(ω,t)|,
(10)

RSAGain(ω,t)=|H12(ω,t)|.
(11)

Since two major rhythms in cardiovascular variability analysis are the one occurring at the frequency of the Mayer waves (LF, 0.04–0.15 Hz) and the one triggered by respiration (HF, 0.15-0.5 Hz, ±0.04 Hz around the respiratory rate) [1], we can compute the power or the gain across these frequencies over time for both (10) and (11).

Hence, from (10) and (11) we can estimate the relevant (instantaneous and mean) statistics, such as the LF power, the HF power, and the RSA gain in HF. We may also compute the dynamic LF/HF power ratio. A small (or large) LF/HF ratio indicates relatively predominant vagal (or sympathetic) control [1].

D. Adaptive Point Process Filtering

In practice, we can bin a continuous-time point process with a certain bin size Δ. The bin size has to be small enough to not only contain one event at most inside each bin, but also to characterize the dynamics at the timescale of interest. It is known from point process theory [2], [3], [14] that, the conditional intensity function λ(t) (CIF) is related to the inter-event probability p(t) with a one-to-one relationship:

λ(t)=p(t)1uttp(τ)dτ.
(12)

The estimated CIF can be used to evaluate the goodness-of-fit of the probability model for the heart beat dynamics. The quantity γ(t)Δ yields approximately the probability of observing a beat during the [t,t + Δ] interval.

Let ξ denote a vector that contains all of unknown parameters in any parametric probability model (in the case of inverse Gaussian model, ξ=[{ai}i=0p,{bj}j=1p,θ]T), we can recursively estimate them via adaptive point process filtering [3]:

ξk|k1=ξk1|k1
(13)

Pk|k1=Pk1|k1+W
(14)

ξk|k=ξk|k1+Pk|k1(logλk)[nkλkΔ]
(15)

Pk|k=[Pk|k11+λkλkTΔλk2logλk[nkλkΔ]]1
(16)

where P and W denote the parameter and noise covariance matrices, respectively; Δ =0.005s denotes the time bin size; Δλk=λkξk and Δ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.

Remark: To avoid numerical problems that might occur in the matrix inverse and to increase the numerical stability, one can replace the original posterior covariance update (16) with a Fisher's scoring step [30]:

Pk|k=[Pk|k11+λkλkTΔλk]1.
(17)

In addition, if the prediction covariance is badly conditioned, the prediction covariance is retained as the posterior covariance.

E. Goodness-of-fit Tests

Model goodness-of-fit is assessed based upon the time-rescaling theorem [14]. Given a point process specified by J discrete events:0 < u1 < … < uj < T, define the random variables zj=uj1ujλ(τ)dτ for j=1,2,…,J−1. Then the random variables zj s are independent, unit-mean exponentially distributed. By introducing the variable of transformation vj = 1 − exp(−zj), then vjs are independent, uniformly distributed within the region [0, 1]. Let gj = Φ−1 (vj) (where Φ(•) denotes the cumulative density function (cdf) of the standard Gaussian distribution), then gj s will be independent standard Gaussian random variables. The Kolmogorov-Smirnov (KS) test} is used to compare the cdf of vj against that of the random variables uniformly distributed in [0, 1]. The KS statistic is the maximum deviation of the empirical cdf from the uniform cdf. To compute it, the vj s are sorted from the smallest to the largest value, and plotted against values of the cdf from the uniform density defined as j0.5J. Ideally, if the model is correct, the points should lie on the 45° line, and the 95% confidence interval lines are y=x±1.36(J1)1/2. The KS distance, defined as the maximum distance between the KS plot and the 45° line, is used to measure lack-of-fit between the model and the data.

In addition, we also compute the autocorrelation function of the gj s: ACF(m)=1Jmj=1Jmgjgj+m. If the gjs are independent, they are also uncorrelated; hence, ACF(m) shall be small (around 0 and within the 95% confidence interval 1.96(J1)1/2) for all values of m.

III. Results

In this section, we conduct a probabilistic analysis of heart beat data with the stochastic point process paradigm. The major advantage of casting the heart beat interval within the point process framework is to allow for the possibility to model and evaluate the instantaneous statistics of HR, HRV, and RSA, at arbitrary fine time resolution.

We start by validating the model with synthetic data and illustrating a comparison with estimates from a window-based algorithm to highlight the fast tracking ability of the point process filter. After describing experimental protocol and model initialization, we confirm the better performance of the inverse Gaussian model with respiration as covariate as the model with overall best fit. We then focus on instantaneous RSA estimation from the real data and illustrate several examples where the point process filter estimates novel dynamic signatures of RSA. Finally, we further validate our indices by presenting a concise group study reporting significant differences of RSA with age and changes in posture or drug.

A. Simulated RSA Study with Synthetic Data

To illustrate and validate our proposed algorithm in estimating the RSA, we first conduct a simulation study with synthetic data. Here, the RR and RP time series are generated by a closed-loop bivariate AR(8) model [5]. All RP measurements are assumed to be synchronized at the time of heart beats. The simulated time series last about 23 minutes.

The instantaneous R-R interval is co-modulated by its own history as well as the RP series:

RRk=a0+i=18aiRRki+j=18bjRPkj+εk
(18)

where εk denotes a Gaussian random noise component with zero mean and standard deviation 20. The synthetic RR and RP time series were generated in a way that the RR has one LF at 0.1 Hz and one HF at 0.2 Hz (they are both fixed along time). The RSA gain achieves the maximum gain 1 around 0.2 Hz for the first 10 min and then increases by a factor of 2 in the remaining time. Spectral analysis also confirmed that these two time series have the maximum coherence at 0.2 Hz.

We used the inverse Gaussian distribution as the probability model for the heart beat interval, with its instantaneous mean modeled by (7). The initial bivariate AR coefficients are estimated from the first 100 seconds of the R-R and RP measurements. The noise covariance matrix W in (14) was set as W = diag{3×10−3, 1× 10−7,…, 1×10−6,…, 1×10−2}. Figure 1 shows the estimated time courses of several statistical indices from one simulated heart beat RR and RP time series.

Figure 1
Estimated time courses of μRR, σRR, μRR (instantaneous HR), σHR (instantaneous HRV), and RSA gain statistics for one Monte Carlo realization of the simulated heart beat data. In the first panel, the superimposed curve is ...

As exemplary comparison, we implemented a recursive least-squares (RLS) filter [21] using a moving window-based method (based on RR and RP time series) to estimate the time-varying bivariate AR coefficients. The forgetting factor of the RLS filter was chosen to be 0.98 in order to balance well the bias-variance tradeoff [12]. Specifically, we conducted 40 Monte Carlo experiments for both adaptive point process and RLS filters and compared their averaged RSA estimates (see the top panel of Fig. 2). Table II summarizes the comparison of Monte Carlo mean and standard deviation (std) statistics of these indices in the first and second half of the data. The mean values of the point process algorithm seem to be more accurate, with an averaged variance which is higher by a factor of 3. The lower variance values for the RLS are mainly attributable to the windowing effect. A careful examination of the estimates reveals that the adaptive point process filter has a faster tracking performance (to approach the true RSA gain value, i.e., from 1 to 2 ms/l in our simulation example) than the RLS filter (see the bottom panel of Fig. 2). Our algorithm reaches the 95\% lower threshold towards the expected value at 615s (15s after the abrupt change in gain), whereas the RLS estimate reaches the same threshold only at 720s (120s after the change). This is not very surprising, since the RLS filter uses an intrinsic exponential moving window-based method and updates values only at the beats, in a much greater timescale than the point process filter [6]. These characteristics also yield to a lack of dynamics across time, and a slow tracking performance to reach the steady-state of the RSA gain. Of note, if we pass a smoothing window on our estimates, we could get similar trends as in the RLS case, where such window is implicitly included in the estimation process. In other words, RLS may be considered a window-smoothed, unevenly sampled version of the point process filter.

Figure 2
The comparison between the adaptive point process filter (dark black trace) and RLS filter (red trace) in the estimated RSA curve (averaged over multiple independent Monte Carlo runs). Top panel: the complete trace. Bottom panel: the zoom-in trace from ...
TABLE II
The Monte Carlo RSA gain statistics (mean ± std, time and trial averaging over 40 independent realizations) comparison computed from the adaptive point process and RLS filters in the simulation experiments.

In summary, in light of many experimental observations (including additional simulations not reported here) it appears that not only the point process filter produces a better characterization in the instantaneous estimates of HR and HRV indices (see [3]), but its bivariate extension can simultaneously provide an accurate estimate of the RSA gain (see Monte Carlo comparison). More importantly, we have illustrated how the adaptive point process filter is capable of tracking sharp dramatic changes in RSA gain. As a further validation for the model fit, the KS plot and the autocorrelation plot were generated and shown in Fig. 3. As seen from the figure, the KS fit is within the 95% confidence interval, indicating the model used here is quite satisfactory; in the autocorrelation plot, nearly all of points are also located inside the 95% confidence bounds.

Figure 3
The KS plot and the autocorrelation plot for the simulated heart beat data

B. Experimental Protocol: Autonomic Blockade

Parasympathetic tone usually dominates in healthy, resting individuals. When a resting subject is given atropine (ATR, a muscarinic receptor antagonist that blocks parasympathetic effects), HR usually increases substantially. In contrast, if the subject is given propranolol (PROP, a receptor antagonist that blocks sympathetic effects), HR usually decreases only slightly. If both sympathetic and parasympathetic effects are blocked, HR is called the instrinsic heart rate [9].

The experimental data considered in our study were previously presented in [31], [32], [33]. A schematic diagram is shown in Fig. 4 to describe different stages of the protocol. In each epoch, 5 min segments of continuous ECG and lung volume were recorded. In the drug administered state, either ATR (0.04 mg/kg iv over 5 min, parasympathetic blockade) or PROP (0.2 mg/kg iv over 5 min, sympathetic blockade) was delivered to the subject. In the double blockade (DB), the inputs from both sympathetic and parasympathetic branches of the autonomic nervous system were suppressed. Subjects were randomized to first receive either ATR or PROP, and then the alternate drug; after 10 min all measurements were repeated. All segments were recorded in a steady state condition (or as stationary as possible). R-wave spike trains were detected from 360 Hz ECG recordings. The lung volume data were digitally recorded at 360 Hz, measuring the calibrated outputs corresponding to rib cage and abdominal compartment volume changes associated with respiration. After verifying the Nyquist frequency condition for the R-R intervals, the long volume recordings were resampled at the beat times.

Figure 4
Diagram of the autonomic blockade protocol.

A total of 17 healthy volunteers participated in the study. Fifteen subjects (6 young subjects and 9 old subjects) are included in the present study. Here, for convenience, we have renamed the subjects with numbers. [2Comparing with the original study [31]: 10=NO, 11=AS, 12=CA, 13=DB, 14=GB, 15=GJ, 16=RB, 20=TM, 21=EO, 22=GH, 23=HU, 24=JW, 25=KN, 26=ME, 27=RN.] More specifically, subjects 10–16 belong to the ATR group; and subjects 20–27 belong to the PROP group. In addition, subjects 11, 12, 14, 16, 21, 22, 24, 26, 27 belong to the old age group; and subjects 10, 13, 15, 20, 23, 25 belong to the young age group. Figure 5 shows a snapshot of R-R intervals and respiratory recordings from a representative subject with varying conditions (in this case propranolol is administered first). As expected, note the general RR decrease, both in control and PROP, accompanied by an increase in ventilation when the subject is standing (observed in only some of the subjects). As atropine is administered, the intrinsic HR shows the expected acceleration with markedly reduced dynamics, which are not affected by change in posture. Of note, in some subjects we observed mild increases in HR with change of posture during double blockade, this may be due to waning of the first administered drug effects: in these cases, even a very mild autonomic modulation may overcome any mechanical effect on the SA node [11]. Going back to our representative subject, note the high nonstationarities in the control PROP epoch, with important bradycardic events, supposedly due to exclusive vagal modulation. This example summarizes well some of the main issues involved in HRV studies that motivated our research: (i) presence of fast dynamics, (ii) presence of high nonstationarities even in supposedly steady-state conditions, (iii) dynamic changes in ventilation to be accounted for accurate RSA assessment.

Figure 5
A snapshot of R-R intervals (in msec) and lung volume respiration measures (RP, calibrated and zero mean) under 6 different conditions (subject 20). Top 6 panels: supine posture. Bottom 6 panels: upright posture. Note that all RR (as well as all RP) plots ...

C. Model and Parameter Initialization

The order of the AR model was determined based on the Akaike information criterion (AIC) (by pre-fitting a subset of the data) as well as the KS distance in the post hoc analysis. In all univariate AR cases, the order p = 8 was chosen from {2,4,6,8,10}. In bivariate AR analyses, the order p = q = 8 was used. The initial AR coefficients are estimated by solving the Yule-Walker equation using about 60–80 seconds of the initial recordings [4].

D. Model Comparison

In order to choose a reasonable subset of parametric models we performed a preliminary histogram analysis and probability fit for the recorded R-R interval time series (see Fig. 6 result for the representative subject in Fig. 5). We consequently selected four characteristic probability structures: gamma, Gaussian, lognormal, and inverse Gaussian.

Figure 6
Histogram analysis and probability fit for the control and double blockdade conditions in supine position (subject 20, see Fig. 2). In the probability fit plots, if the data fit the tested probability distribution, the data points will match the straight ...

The comparative results of the four probability models considered for the heart beat data are presented in Table III. The inverse Gaussian model achieves the overall best fit in terms of smaller KS distance, especially during the control and PROP epochs, in both supine and upright positions. The lognormal model achieves better performance during the DB epochs. The gamma model has the worst performance among the four probability models tested here. All models perform rather unsatisfactorily during the ATR epochs.

TABLE III
Performance comparison between four univariate probability models. The number indicates the total number of epochs that specific model achieves the lowest KS distance or the best goodness-of-fit among all epochs of 15 tested subjects. The second last ...

When the bivariate (inverse Gaussian) model was further applied to the heart beat time series together with the lung volume recordings, it was found that the inclusion of the respiration covariate improves the KS fit (second last column, Table III) in the majority of the subjects under all three pharmacological conditions (ATR, PROP, and DB). Figure 7 illustrates a comparative example using (2) and (7) in the “upright+PROP" condition.

Figure 7
Comparison of inverse Gaussian models with the mean as the univariate (top row) and bivariate (bottom row) AR models for subject 20 (upright, PROP). Left panel: estimated time-varying probability density function of the instantaneous μRR(t). ...

In a few subjects that were examined, it was found that the inclusion of respiration measurements does not improve the KS fit. This might be due to the fact that the R-R intervals history itself might be already sufficient to characterize the instantaneous mean statistic for the inverse Gaussian probability distribution (especially in the control cases), and adding covariates makes the on-line adaptation more challenging.

As seen in our experiments, the bivariate model has the best performance in control and PROP conditions, followed by the DB condition, and has the worst fit in the ATR condition. Specifically, for the goodness-of-fit in the control condition, about 85% of the time-rescaled points were inside the 95% confidence interval of the KS plot, which implies that roughly 85% of the beats are correctly predicted by the model (last column of Table III). It is also important to pinpoint that if the same criterion were to be applied to any other window-based model, there would not be a single case where it be satisfied. More specifically, the KS curves would fall far away from the 45° diagonal, with consequent significantly greater KS distance values (see [2]).

E. Instantaneous RSA Analysis

Our point process bivariate framework allowed for estimation of instantaneous measures of RSA. For the protocol analysis, an updating delta interval of 5 ms was chosen. Figure 8 shows data and results from subject 25 for three epochs: Control Supine (CS), Control Upright (CU), and Propranolol Supine (PS). In addition to the expected mean RSA decrease with change in posture, there is an expected mean RSA increase with propranolol administration (when compared to the respective control epoch). The new estimates show novel interesting dynamics. In particular, note the sharp increases in RSA around 160s, 190s, and 230s in CS, the marked oscillatory trend in CU, and the sharp decreases (around 105s, 160s, 185s, and 230s) from a saturated maximum RSA level at 440 ms/l in PS. All these dynamic signatures could not possibly have been inferred by looking at the RR and RP time series, or by applying stationary or window-based analyses. These observations point at our estimates as effective instantaneous measures for a unique novel characterization of vagal modulation at small timescale levels.

Figure 8
The R-R interval, lung volume respiration (RP) measure ((adjusted to zero mean), and estimated dynamic RSA mean gain (0.15-0.5 Hz) in 3 consecutive epochs (subject 25). From control supine to control upright, RR decreases and RP increases signficantly, ...

Figure 9 and Figure 10 plot the instantaneous RSA mean gain (HF: 0.15-0.5 Hz) for all recordings from an old subject from the ATR group and a young subject from the PROP group, respectively. Since the plots are scaled to evidence the dynamics, mean values are reported in the figure to allow for a more direct comparison among epochs. The first observation is that whenever atropine is administered, together with the sharp mean RSA decrease due to the absence of parasympathetic modulation, the RSA dynamics almost disappear to the point that quantization effects due to the bin resolution become evident. On the other hand, interesting trends can be observed also in these two subjects in the presence of parasympathetic modulation. In particular, we observe the saturation effects in Control-Upright for the ATR subject, and in PROP-Upright for the PROP subject. Note also the faster oscillations in the Control-Supine epoch for the young subject (Fig. 10). Again, these trends could not be observed with stationary or window-based analysis.

Figure 9
The estimated instantaneous RSA gain (unit: ms/l) in the HF (0.15-0.5 Hz) range (subject 14 from the old/ATR group). The number in each subplot indicates the mean value of the RSA gain averaged over the entire recording (which is computed using the unnormalized ...
Figure 10
The estimated instantaneous RSA gain (unit: ms/l) in the HF (0.15-0.5 Hz) range (subject 20 from the young/PROP group). The number in each subplot indicates the mean value of the RSA gain averaged over the entire recording (which is computed using the ...

F. Statistical Analysis in Group Study

Upon computing the RSA gain statistics (Table IV), we investigate if there is any statistical difference between posture/pharmacological/age conditions, despite the fact that only a limited number of subjects are available in the present study. Specifically, we applied a nonparametric Mann-Whitney test (also known as the rank sum test) to compare two independent samples. In the case of multiple comparison, we also adjusted the p-value according to the Bonferroni correction. We found that within the young age group, the supine vs. upright posture effect (p<0.01) as well as the control vs. DB drug effect (p <0.01) both show statistical significance. These two effects fail to reach significance within the old age group, possibly due to high inter-subject variability. In the comparison between two age groups in the Control-Supine condition, we also found a statistically significant difference (p <0.01). Our observations are consistent with previous findings pointing at a weaker autonomic control with increasing age [31], [32].

TABLE IV
Mean and standard deviation of RSA gain (unit: ms/l) across subjects in different conditions.

IV. Discussion and Conclusion

We have presented an extended point process paradigm for human heart beat intervals with respiration as covariate to assess autonomic control as quantified by RSA. Our method is validated by Monte Carlo simulations using synthetic data, and the new fast-tracking instantaneous RSA assessment is illustrated as applied to the experimental autonomic blockade recordings. These examples reveal novel interesting dynamic trends reflective of the nonstationary nature of cardiovascular control, whereas simple summary statistics confirm established findings related to RSA measures across different posture/pharmacological/age conditions. Several points below are worth further discussion.

A. Choice of Probability Model

In modeling the heart beat interval during the control epochs, the inverse Gaussian model achieves the best performance, which is in agreement with our earlier claims [2], [3]. The Gaussian model achieves a similar performance since, when the random variable's mean is much greater than its variance, the inverse Gaussian can be well approximated by a Gaussian shape. In modeling the pharmacological autonomic blockade, the inverse Gaussian model is more suited for PROP than ATR---this suggests that the markedly reduced “sympathetic-driven" variability requires more effort for modeling in the absence of parasympathetic modulation. The lognormal model is better fitted for the double blockade---this is partly due to the fact that during DB the lognormal model is more robust in characterizing the significant drop in HRV. With inclusion of the RP covariate into the model, we not only achieve a more accurate physiological model of cardiovascular control (as reflected by a better goodness-of-fit), but we are also able to explicitly monitor the respiratory effects and evaluate the instantaneous RSA gain. In this particular autonomic blockade protocol, the reduced KS fit in the ATR epochs (for most subjects studied in the protocol) still leaves us challenges in choosing appropriate probability models for the heart beat interval.

In fitting experimental data, a perfect KS fit was not always expected, partially because the complex organization of cardiovascular control calls for consideration of several other physiological measurements, such as arterial blood pressure, central venous pressure, and vascular resistance. The higher lack of fit observed in the absence of vagal modulation was also expected, given the markedly reduced variability in these cases. It is important to stress again that the fits for any of our models achieved a far better result when compared to any window-based methods [2]. The possibility of a paradigm that allows for a flexible and adaptive choice of the heart beat interval model will be the subject of future study (e.g., [23]).

B. RSA and Model Identifiability in Adaptive Estimation

The RSA gain is a useful index of vagal control that often correlates with R-R interval modulation. This has been confirmed by our experiments in both simulated and real data. In the latter case the RSA values expectedly decrease in the upright position as compared to supine, and they show significant lower values in DB when vagal activity is absent among healthy young individuals.

The computation of the instantaneous RSA gain depends only on the AR coefficients {aj} and {bj}, which are estimated by the adaptive point process filter. Note that the adaptive point process filter also adapts the variance (or shape) parameter of the probability model at every step; consequently, it changes the model in both mean μRR and variance (or shape) parameters to fit the data.

Since the change rates (within every time bin Δ) of the parameters in vector ξ are determined by their respective random-walk noise variance components (i.e., the diagonal components in matrix W), we might not be able to recover the exact model parameters in ξ in an online manner. Furthermore, even for the bivariate AR model (7) alone (for simplicity assuming the shape parameter of the inverse Gaussian model is fixed), the values of {aj} and {bj} are initialized by a batch least-squares method; however, the model ambiguity problem arises in all adaptive (or online) methods---in other words, there exist many solutions for {aj} and {bj} that can produce an identical mean value μt in (7). For the reasons mentioned above, we can't exclude the possibility of producing a bias in estimating the instantaneous RSA gain in real data.

C. Model Extension

It is noted that thus far the model of μRR, transfer function, and frequency response analysis are all limited by the assumption of a linear system. It is certainly our interest to investigate the nonlinear coupling and nonlinear modulation effects among the cardiovascular/cardiorespiratory systems. Some preliminary work along this direction has been conducted [16]. In the meanwhile, our proposed point process framework is general and similar methodology can be applied to investigate the interaction between heart beat intervals and other cardiovascular covariates, such as the systolic blood pressure [17], for the purpose of studying other cardiovascular functions of interest.

D. Conclusion

To conclude the paper, the probabilistic point process framework is powerful in estimating instantaneous heart beat dynamics involved in autonomic control. In line with the most conventional guidelines [1], [10], our paradigm resolves the old conflict between heart period and heart rate, obviates the need of interpolation with consequent possible aliasing problems, allows for instantaneous measures at virtually any time resolution, and overcomes nonstationarity issues associated with window-based estimation approaches. More importantly, the point process models can be rigorously validated by goodness-of-fit test. Furthermore, our experimental results confirm earlier established findings regarding important physiological mechanisms involved in cardiovascular control, such as RSA, and they also reveal interesting dynamic trends across different posture/pharmacological/age conditions [27], [31], [33]. Our point process approach provides a novel assessment of RSA that we will further validate in a broader range of experimental contexts, and we will use to help answer important remaining questions involving RSA quantification and the role of respiratory activity in cardiovascular control physiology. The dynamic statistical indices (such as HR, HRV, and RSA gain) computed from our point process framework provide the basis for potential real-time indicators for ambulatory monitoring and instantaneous assessment of autonomic control in clinical practice.

Acknowledgment

Preliminary version of this work was presented in Proc. ICASSP'08, Las Vegas, USA. This work was supported by NIH Grants R01-HL084502, R01-DA015644 and DP1-OD003646.

The authors thank Dr. Janice B. Schwartz (University of California, San Francisco) and Dr. Garrett B. Stanley (Georgia Institute of Technology and Emory University) for providing the experimental data used in this study. Recordings were performed during the tenure of a research fellowship from the American Heart Association, California Affiliate (G. B. Stanley), with further support by National Institutes of Health (NIH) Grants GM-26691 and AG-9550. The work was performed at the University of California, San Francisco, General Clinical Research Center supported with funds provided by the Division of Research Grant RR-79. We thank two anonymous reviewers for their valuable comments that help to reshape and improve the presentation of the manuscript.

References

1. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, “Heart rate variability” Circulation. 1996;vol. 93(no 5):1043–1065. [PubMed]
2. 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;vol. 288:424–435. [PubMed]
3. Barbieri R, Brown EN. Analysis of heart beat dynamics by point process adaptive filtering. IEEE Trans. Biomed. Eng. 2006;vol.53(no 1):4–12. [PubMed]
4. 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 bivariate autoregressive techniques. Annals of Noninvasive Electrocardiology. 1996;vol. 3:264–277.
5. Barbieri R, Bianchi AM, Triedman JK, Mainardi LT, Gerutti S, Saul JP. Model dependency of multivariate autoregressive spectral analysis. IEEE Mag. Engr. Med. Biol. 1997;vol. 16(no 5):74–85. [PubMed]
6. Barbieri R, Parati G, Saul JP. Closed- versus open-loop assessment of heart rate baroreflex. IEEE Mag. Engr. Med. Biol. 2001;vol.20(no 2):33–42. [PubMed]
7. Baselli G, Bolis D, Cerutti S, Freschi C. Autoregressive modeling and power spectral estimate of R-R interval time series in arrhythmic patients. Comput Biomed Res. 1985;vol. 18(no 6):510–530. [PubMed]
8. Batzel JJ, Novak V, Kappel F, Olufsen MS, Tran HT. Introduction to the special issues: short-term cardiovascular-respiratory control mechanisms. Cardiovasc. Eng. 2008;vol. 8:1–4. [PMC free article] [PubMed]
9. Berne RM, Levy MN, Loeppen BM, Stanton BA. Physiology. 5th ed. Elsevier; 2004. “Regulation of the heartbeat,” Chap. 17.
10. Bernston GG, Bigger JT, Eckberg DL, Grossman P, Kaufmann PG, Malik M, Nagaraja HN, Porges SW, Saul JP, Stone PH, van der Molen MW. Heart rate variability: origins, methods, & interpretative caveats. Psychophy. 1997;vol. 34:623–648. [PubMed]
11. Bernardi L, Keller F, Sanders M, Reddy PS, Griffith B, Meno F, Pinsky MR. Respiratory sinus arrhythmia in the denervated human heart. J. Appl. Physiol. 1989;vol. 67:1447–1455. [PubMed]
12. Bianchi AM, Mainardi L, Petrucci E, Signorini MG, Mainardi M, Cerutti S. Time-variant power spectrum analysis for the detection of transient episodes in HRV signal. IEEE Trans. Biomed. Eng. 1993;vol. 40(no 2):136–144. [PubMed]
13. Brown EN, Barbieri R, Ventura V, Kass RE, Frank LM. The time-rescaling theorem and its application to neural spike data analysis. Neural Computation. 2002;14:325–346. [PubMed]
14. Brown EN, Barbieri R, Eden UT, Frank LM. Likelihood methods for neural data analysis. In: Feng J, editor. Computational Neuroscience: A Comprehensive Approach. London: CRC; 2003. pp. 253–286.
15. Chen Z, Brown EN, Barbieri R. Proc. ICASSP’2008. Las Vegas, USA: A study of probabilistic models for characterizing human heart beat dynamics in autonomic blockade control; pp. 481–484. [PMC free article] [PubMed]
16. Chen Z, Brown EN, Barbieri R. Characterizing nonlinear heartbeat dynamics within a point process framework. Proc. IEEE 30th Annual Conf. Engineering in Medicine and Biology (EMBC’08); Vancouver, Canada. pp. 2781–2784. [PMC free article] [PubMed]
17. Chen Z, Brown EN, Barbieri R. A point process approach to evaluate baroreflex gain and cardiovascular control. Proc. Computers in Cardiology (CinC’08); Bologna, Italy. 1995. pp. 805–808.
18. Chhikara RS, Folks JL. New York: Marcel Dekker; 1989. The Inverse Gaussian Distribution: Theory, Methodology, and Applications.
19. Christini DJ, Bennett FM, Lutchen KR, Ahmed HM, Hausdofi JM, Oriol N. Application of linear and nonlinear time series modeling to heart rate dynamics analysis. IEEE Trans. Biomed. Eng. 1995;vol. 42:411–415. [PubMed]
20. Eckberg DL. Human sinus arrhythmia as an index of vagal cardiac outflow. J. Appl. Physiol. 1983;vol. 54:961–966. [PubMed]
21. Haykin S. 4th ed. Prentice Hall: 2002. Adaptive Filter Theory.
22. Hirsch JA, Bishop B. Respiratory sinus arrhythmia in humans: How breathing pattern modulates heart rate. Am. J. Physiol. 1981;vol. 241:H620–H629. [PubMed]
23. Iyengar S, Liao Q. Modeling neural activity using the generalized inverse Gaussian distribution. Biol. Cybern. 1997;vol. 77:289–295. [PubMed]
24. Katona PG, Jih F. Respiratory sinus arrhythmia: noninvasive measure of parasympathetic cardiac control. J. Appl. Physiol. 1975;vol. 39:801–805. [PubMed]
25. Malpas SC. Neural influences on cardiovascular variability: possibilities and pitfalls. Am. J. Physiology Heart Cicr. Physiol. 2002;vol. 282:6–20. [PubMed]
26. Motte S, Mathieu M, Brimioulle S, Pensis A, Ray L, Ketelslegers JM, Montano N, Naeije R, van de Borne P, Entee KM. Respiratory-related heart rate variability in progressive experimental heart failure. Am. J. Physiology Heart Cicr. Physiol. 2005;vol. 289(no 4):1729–1735. [PubMed]
27. Saul JP, Berger RD, Chen MH, Cohen RJ. Transfer function analysis of autonomic regulation II. Respiratory sinus arrhythmia. Am. J. Physiology Heart Cicr. Physiol. 1989;vol. 256(no 25):153–161. [PubMed]
28. Saul JP, Cohen RJ. “Respiratory sinus arrhythmia,” Chap. 30 in Vagal Control of the Heart: Experimental Basis and Clinical Implications. In: Levy MN, Schwartz PJ, editors. Armonk, NY: Futura Publishing Inc; 1994.
29. Pinna GD, Maestri R, La Rovere MT, Gobbi E, Fanfulla F. Effect of paced breathing on ventilatory and cardiovascular variability parameters during short-term investigations of autonomic function. Am. J. Physiology Heart Cicr. Physiol. 2006;vol. 290(no 1):424–433. [PubMed]
30. Srinvasan L, Eden UT, Mitter SK, Brown EN. General-purpose filter design for neural prosthetic devices. J. Neurophysiol. 2007;vol. 98:2456–2475. [PubMed]
31. Stanley GB, Verotta D, Craft N, Siegel RA, Schwartz JB. Age and autonomic effects on interrelationships between lung volume and heart rate. Am J. Physiol. Heart Cicr. Physiol. 1996;vol. 270(no 5):1833–1840. [PubMed]
32. Stanley GB, Verotta D, Craft N, Siegel RA, Schwartz JB. Age effects on interrelationships between lung volume and heart rate during standing. Am J. Physiol. Heart Cicr. Physiol. 1997;vol. 273:2128–2134. [PubMed]
33. Stanley GB, Poolla K, Siegel RA. Threshold modeling of autonomic control of heart rate variability. IEEE Trans. Biomed. Engr. 2000;vol. 47(no 6):1147–1153. [PubMed]
34. Taylor JA, Myers CW, Halliwill JR, Seidel H, Eckberg DL. Sympathetic restraint of respiratory sinus arrhythmia: implications for vagal-cardiac tone assessment in humans. Am. J. Physiology Heart Cicr. Physiol. 2001;vol. 280(no 6):2804–2814. [PubMed]
35. White JA, Semple E, Norum RA. Sinus arrhythmia as an estimate of maximal aerobic power in man. J. Physiol. 1990;vol. 432:41.