Search tips
Search criteria 


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

Point Process Time-Frequency Analysis of Respiratory Sinus Arrhythmia under Altered Respiration Dynamics


Respiratory sinus arrhythmia (RSA) is largely mediated by the autonomic nervous system through its modulating influence on the heartbeat. We propose an algorithm for quantifying instantaneous RSA as applied to heart beat interval and respiratory recordings under dynamic respiration conditions. The blood volume pressure derived heart beat series (pulse intervals, PI) are modeled as an inverse Gaussian point process, with the instantaneous mean PI modeled as a bivariate regression incorporating both past PI and respiration values observed at the beats. A point process maximum likelihood algorithm is used to estimate the model parameters, and instantaneous RSA is estimated by a frequency domain transfer function approach. The model is statistically validated using Kolmogorov-Smirnov (KS) goodness-of-fit analysis, as well as independence tests. The algorithm is applied to subjects engaged in meditative practice, with distinctive dynamics in the respiration patterns elicited as a result. Experimental results confirm the ability of the algorithm to track important changes in cardiorespiratory interactions elicited during meditation, otherwise not evidenced in control resting states.

I. Introduction

A large number of autonomic and hemodynamic parameters are influenced by respiratory activity. Among them, respiratory sinus arrhythmia (RSA) can be defined as the variations in the heart rate which occur simultaneously with the respiration cycle [1]. At typical resting respiratory frequencies, heart rate increases during inspiration, and decreases during expiration. RSA is usually considered as an indirect measure of vagal cardiac control, although at lower respiratory frequencies sympathetic cardiac control contributes to RSA as well [1]. A key problem in cardiorespiratory engineering is to efficiently and accurately quantify RSA under different physiological conditions, such as variations in respiration and changes in sympathetic-parasympathetic inputs. A solution to this issue could yield critical insights into the mechanisms involved in short-term cardiorespiratory coupling [1] [2].

In previous work [3] [4], RSA was defined using simple time domain measures of beat interval series. Filtering and transfer function approaches were also used in quantifying RSA [5] [6], and a bivariate autoregressive model was further proposed to estimate the time-varying RSA gain [7] [8]. As most of these methods are not able to overcome stationarity issues and estimate fast changes in RSA at arbitrarily small time scales, a point process framework for heart beat dynamics [9] [10] has been proposed to assess RSA within an adaptive point process filtering algorithm [11].

In this paper, we use a novel local maximum likelihood method within a point process framework to allow for instantaneous estimation of RSA. Importantly, as measures based on the traditional subdivision in oscillatory frequency components might not be reliable in the presence of non-stationary respiratory patterns, we further propose a new method for dynamically selecting the RSA gain within the transfer function spectrum, based on a time-frequency characterization of the respiratory cycle and the time-varying coherence between respiration and heart beat series (pulse intervals, PI). Such combined method is capable of computing reliable, instantaneous estimates of RSA by accounting for rapid dynamic changes in both respiration patterns and autonomic inputs. The new algorithm is validated on simulated data, as well as applied to recordings from subjects practicing meditation.

II. Methods

A. Point Process Model of Heart Beat Interval Dynamics

Integrate and fire models are regularly used to simulate heart beats, and such model postulates that the resulting times between two firing events (the PIs) have statistical properties of an inverse Gaussian process [9]. Additionally, autonomic inputs to the SA node are part of the cardiovascular control circuitry, thus the PI variations are dynamic or time-varying [9]. For this reason, we here model the pulse intervals as a history dependent, inverse gaussian point process model with time-varying model parameters. Assume in a given observation interval (0,T], K successive pulses are recorded: 0 < u1 < u2 < … < uKT. Given any pulse-timing uk, the waiting time until the next pulse-timing (ie. Next PI), obeys a history dependent inverse gaussian probability density f(t) given by

equation M1

where t is any time satisfying t > uk, and μPI (t) > 0 is the mean of the distribution, which is an estimation of the instantaneous mean PI. θ(t) > 0 is the shape parameter of the inverse gaussian distribution.

Heart rate is often defined as the reciprocal of PIs [12]. For PIs measured in seconds, r = c/(tuk), where c = 60s/min, gives the estimation of heart rate in beats per minute (bpm). By using a standard change of variables formula [13], mean and standard deviation of the heart rate probability density are given by [9]

equation M2

equation M3

and they respectively define the instantaneous estimates of heart rate and heart rate variability.

B. Bivariate Model for RSA analysis

The influence of past autonomic inputs and respiration activity on the PIs are incorporated into the model by defining a bivariate regression on the mean of the point process probability density (1),

equation M4

Note that the respiration signal (RP) is sampled at the pulse timings, so that both respiration and PIs are synchronized.

RSA can then be defined as the transfer function from RP to PI,

equation M5

We propose two methods of estimating the RSA gain from the above transfer function. First, the time-varying respiration spectrum PRP(ω,t) is used to estimate the frequency ωRP(t) where maximum respiration power is concentrated at each time instance, i.e.,

equation M6

Then, RSA gain can be estimated by evaluating (5) at ωRP,

equation M7

Second, we evaluate the RSA gain at the frequency where maximum interaction between PIs and respiration occur. In this regards, we use the time-varying autoregressive coherence spectrum Coh(ω, t) [5] to estimate the frequency ωcoh(t) where coherence is maximum, i.e.,

equation M8

and the RSA gain is evaluated at ωcoh,

equation M9

C. Local Maximum Likelihood Estimation

A local maximum likelihood method [14] was used to estimate the unknown time-varying parameter set ξ = {{ak}pk=0, {bk}qk=1, θ}. In estimating ξ at time t, we take a local likelihood interval (tl, t], where l is the length of the local likelihood observation interval. Within (tl, t], we may observe n pulses, tl < u1 < u2 < … < unt. Then, we consider the local joint probability density of utl:t, where utl:t = { u1,…, un }. Log likelihood of joint probability density is given by

equation M10

where w(tuj) = αtuj, 0 < α < 1, is a weighting function for the local likelihood estimation [14]. The weighting time constant α governs the degree of influence of a previous event observation uj on the local likelihood at time t. Second term of (10) represents the likelihood of partially observed interval since the last observed pulse un (right censoring). To maximize the local log likelihood in (10) we use a Newton- Raphson method, and obtain the local maximum likelihood estimate of ξ. Of note, the time increment Δ for computing the next ξ from t to t+ Δ can be chosen arbitrarily small, thus yielding instantaneous estimates of heart rate and heart rate variability.

Goodness-of-fit of the proposed model was evaluated using a Kolmogorov-Smirnov (KS) test based on the time-rescaling theorem [15]. The test uses the conditional intensity function given by,

equation M11

to transform pulse events into independent observations on the interval [0, 1], and the KS plot allows to test the agreement of the transformed observations and the ideal uniform probability density. The transformed quantiles’ autocorrelation function is further computed to check independence of the transformed intervals.

III. Results

A. Experimental Protocol

Meditation essentially consists of subjects focusing on their breath and passively ignoring everyday thoughts. For this reason, mastering meditation techniques has been often mischaracterized as “simply resting”. To such extent, the present experiment considers two groups of subjects, one group consists of experienced meditators, whereas the other group consists of subjects with no meditation experience. The experiment starts with 6 minutes of control period, followed by 1 minute of fixation, then 24 minutes of meditation, followed by 1 minute of fixation, and finally 6 minutes with the subject silently counting random numbers [16]. The control subjects are asked to simply rest during the “meditation” period. During the experiment, blood volume pressure (BVP) signal, and a belt signal proportional to lung volume changes were recorded at 1 kHz. BVP was used to identify the PIs, and the belt signal was sampled at these pulse timings to obtain the respiration values at the beats.

B. PI and HRV Estimation

First, optimal values for regression orders of the bivariate model p and q, maximum-likelihood interval l, and weighting time constant α were obtained by minimizing the Akaike Information Criterion for maximum likelihood estimation, as well as the KS distance on the KS plot. This empirical optimization yields to p = 4, q = 6, l = 90s, and α = 0.98. The KS plot and autocorrelation function for this final model are shown in Figure 1 for a representative subject.

Fig 1
KS plot (left panel) and Autocorrelation function (right panel) for one experienced meditator. Dashed lines indicate the 95% confidence bounds.

Then, the proposed maximum likelihood point process model was applied to both meditating and control subjects with model parameters as described above. The resulting instantaneous mean PI estimate μPI (t) (4), and HRV index σ2HR(t) (3) for two representative subjects are shown in Figures 2 and and3.3. In particular, sudden variations can be observed in the original PI series of the experienced meditator (Fig. 2). These distinctive contractions are not classifiable as ectopic beats, but could be attributed to fast tachycardic and/or bradycardic events, possibly due to sharp shifts in the sympatho-vagal balance. Sharp increases in the HRV index are consequently observed at corresponding timings. Of note, similar transients were less frequent in the control subject (Fig. 3).

Fig 2
Meditation subject: Instantaneous mean PI (top), and Heart rate variance (bottom) estimated by the point process model. Dashed lines (from left to right) indicate the start times of control, meditation, numbers epochs respectively, and dotted lines indicate ...
Fig 3
Control subject: Instantaneous mean PI (top), and Heart rate variance (bottom) estimated by the point process model.

Average statistics for 6 subjects from each group were performed for mean PI, HRV, and the LF/HF index [12]. Results show no significant changes in all parameters in the control group, whereas a general increase in HRV was observed in the meditating group. Importantly, we observed that the respiratory cycle often tends to slow down during both meditation practice and relaxation, and some of the respiratory variability may fall well within the low frequency (LF) band of the standard HRV frequency band division [12]. As a consequence, the RSA-elicited PI variability would mirror the respiratory spectral distribution, with vagal modulated cardiac dynamics distributed between both LF and HF bands, thus potentially biasing the LF/HF index.

C. Instantaneous RSA Estimation

Instantaneous time-frequency RSA estimates are computed using the bivariate AR point process and transfer function methods described in Section II-B. The algorithms were first tested with artificial series simulating both sudden and gradual shifts in RSA gain, as well as shifts in the main frequency where the two series are coupled. Results (not shown) demonstrate that the proposed algorithm is very robust and stable even under sudden changes in the foremost oscillatory (respiratory) frequency, while still being capable of accurately estimating sudden and gradual changes in the RSA gain.

As applied to the considered experimental data, the first important finding was that both RSAagain(t) and RSAbgain(t) gave very similar estimates in all subjects considered. This finding was somehow expected, as PIs are often highly correlated with respiration activity at the main respiratory frequency, and indicates both RSAagain(t) and RSAbgain(t) as reliable estimates of instantaneous RSA. Results are shown for one experienced meditator in Figure 4, and for a control subject in Figure 5. For both subjects we generally observed a high coherence (> 0.8), with only one evident exception in the experienced meditator where coherence drops below 0.8 at the timings corresponding to some of the fast transients in the heart beat dynamics mentioned before. Here, the significant increase in RSA gain during meditation is clearly noticeable, whereas in the control example RSA gain seems to fluctuate around baseline levels along the entire experiment1.

Fig 4
Meditation subject: (a) Frequency where coherence between PI and RP is maximum (fcoh), (b) Maximum coherence, (c) RSA gain estimated at fcoh by the point process model (nrpu = normalised respiratory units).
Fig 5
Control subject: (a) Frequency where coherence between PI and RP is maximum (fcoh), (b) Maximum coherence, (c) RSA gain estimated at fcoh by the point process model.

Average RSA statistics for 6 subjects from each group are shown in Table I. Results point at a major increase (58%) in RSA during meditation, with high statistical significance (p = 0.031) compared to the baseline epoch. On the other hand, the control group does not show any significant change in RSA during relaxation (p > 0.2). During the silent random number generation phase, RSA values are statistically comparable to baseline levels for both groups.

Table 1
Average RSA gain of 6 subjects each in the meditation and control groups. The percentage change, and Wilcoxon signed rank test results compared to the baseline epoch are also reported.

Differently from the LF/HF assessment, our proposed method is robust to variations in respiration patterns as we rely on the respiratory spectral power distribution and the dynamic coherence between respiration and PIs, without any artificial subdivision in frequency bands. As a further proof that RSA gains are not considerably biased by possible shifts in respiration frequency, we have computed similar indices at a fixed frequency, where coherence is the highest on average along each experiment, and verified that average statistics using these quantities confirm the presented results.

In summary, while traditional HRV indices confirm previous findings pointing at an autonomic assessment not necessarily reflective of quiescent cardiac dynamics [17], and possibly biased by altered breathing patterns [18], our statistical RSA appraisal substantiates previous reports of a marked increase of RSA during meditation [19].

IV. Conclusions

We have proposed a maximum likelihood point process model for instantaneous RSA estimation combined with a time-frequency RSA evaluation based on respiration spectrum and coherence. Our novel framework allows for robust tracking of RSA changes at any time resolution, as well as overcoming potential problems associated with traditional subdivision in standard frequency bands in the presence of distinctive respiratory dynamics. The new algorithm was applied to subjects either practicing meditation, or just asked to relax under equal conditions. Results show a significant increase in RSA under meditation practice which is not evident in the control group, encouraging further investigation into the effects of meditation techniques on cardiovascular control and into the potential benefits of meditation on cardiovascular health. Overall, the dynamic statistical measures computed from our point process framework provide the basis for potential realtime indicators for ambulatory monitoring and instantaneous assessment of autonomic control in clinical practice.


This work was supported by the National Institutes of Health (NIH) under Grant R01-HL084502, Grant R01-DA015644, Grant DP1-OD003646.


1Note that for a fair comparison across subjects, the RSA gain (9) was normalized by the standard deviation of the corresponding respiration signal.


1. Saul JP, Cohen RJ. Vagal control of the heart: Experimental basis and clinical implications. Futura Publishing Co. Inc; 1994. Respiratory sinus arrhythmia; pp. 511–536.
2. Eckberg DL. Human sinus arrhythmia as an index of vagal cardiac outflow. J Appl Physiol. 1983;54:961–966. [PubMed]
3. Hirsch JA, Bishop B. Respiratory sinus arrhythmia in humans: How breathing pattern modulates heart rate. Amer J Physiol. 1981;241:H620–H629. [PubMed]
4. White JA, Semple E, Norum RA. Sinus arrhythmia as an estimate of maximal aerobic power in man. J Physiol. 1990;432:41.
5. Womack BF. The analysis of respiratory sinus arrhythmia using spectral analysis and digital filtering. IEEE Trans Biomed Eng. 1971;18:399–409. [PubMed]
6. Saul JP, Berger RD, Chen MH, Cohen RJ. Transfer function analysis of autonomic regulation. II. Respiratory sinus arrhythmia. Am J Physiol - Heart and Circulatory Physiol. 1989;256(1):H153–H161. [PubMed]
7. 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. Ann Noninvasive Electrocardiol. 1996;1:264–277.
8. Barbieri R, Bianchi AM, Triedman JK, Mainardi LT, Cerutti S, Saul JP. Model dependency of multivariate autoregressive spectral analysis. IEEE Eng Med Bio Mag. 1997;16(5):74–85. [PubMed]
9. Barbieri R, Matten EC, Alabi ARA, Brown EN. A pointprocess model of human heartbeat intervals: new definitions of heart rate and heart rate variability. Am J Physio- Heart and Circulatory Physiology. 2005;288(1):H424–H435. [PubMed]
10. Barbieri R, Brown EN. Analysis of heartbeat dynamics by point process adaptive filtering. IEEE Trans Biomed Eng. 2006;53(1):4–12. [PubMed]
11. Chen Z, Brown E, Barbieri R. Assessment of Autonomic Control and Respiratory Sinus Arrhythmia Using Point Process Models of Human Heart Beat Dynamics. IEEE Trans Biomed Eng. 2009;56(7):1791–1802. [PMC free article] [PubMed]
12. Malik M, Bigger JT, Camm AJ, Kleiger RE, Malliani A, Moss AJ, Schwartz PJ. Heart Rate Variability Standards of Measurement, Physiological Interpretation, and Clinical Use. Am Heart Assoc Circulation. 1996;93(5):1043–1065.
13. Ross SM. Introduction to probability models. 9. Academic Pr; 2007.
14. Loader C. Local regression and likelihood. Springer; 1999.
15. Brown EN, Barbieri R, Ventura V, Kass RE, Frank LM. The time-rescaling theorem and its application to neural spike train data analysis. Neural Comput. 2002;14(2):325–346. [PubMed]
16. Lazar SW, Kerr CE, Wasserman RH, Gray JR, Greve DN, Treadway MT, McGarvey M, Quinn BT, Dusek JA, Benson H, RSL, MCI, Fischl B. Meditation experience is associated with increased cortical thickness. Neuroreport. 2005;16(17):1893–1897. [PMC free article] [PubMed]
17. Peng CK, Henry IC, Mietus JE, Hausdorff JM, Khalsa G, Benson H, Goldberger AL. Heart rate dynamics during three forms of meditation. Int J Cardio. 2004;95(1):19–27. [PubMed]
18. Lehrer P, Sasaki Y, Saito Y. Zazen and cardiac variability. Psychosomatic Medicine. 1999;61(6):812–821. [PubMed]
19. Ditto B, Eclache M, Goldman N. Short-term autonomic and cardiovascular effects of mindfulness body scan meditation. Annals of Behavioral Medicine. 2006;32(3):227–234. [PubMed]