Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Med Biol Eng Comput. Author manuscript; available in PMC 2013 March 1.
Published in final edited form as:
PMCID: PMC3341131

Point Process Time-Frequency Analysis Of Dynamic Respiratory Patterns During Meditation Practice


Respiratory sinus arrhythmia (RSA) is largely mediated by the autonomic nervous system through its modulating influence on the heart beats. We propose a robust algorithm for quantifying instantaneous RSA as applied to heart beat intervals and respiratory recordings under dynamic breathing patterns. The blood volume pressure derived heart beat series (pulse intervals, PIs) are modeled as an inverse gaussian point process, with the instantaneous mean PI modeled as a bivariate regression incorporating both past PIs 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 via a frequency domain transfer function evaluated at instantaneous respiratory frequency where high coherence between respiration and PIs is observed. 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. The presented analysis confirms the ability of the algorithm to track important changes in cardiorespiratory interactions elicited during meditation, otherwise not evidenced in control resting states, reporting statistically significant increase in RSA gain as measured by our paradigm.

Keywords: Respiratory Sinus Arrhythmia, Heart Rate Variability, Meditation, Point Processes, Time-Frequency Analysis

1 Introduction

A large number of autonomic and hemodynamic parameters are influenced by respiratory activity. Among them, respiratory sinus arrhythmia (RSA) is defined as the variations in heart rate during inspiratory and expiratory phases of the respiratory cycle (43, 13). At typical resting respiratory frequencies, heart rate increases during inspiration, decreases during expiration, and respiratory frequency also influences the phase relationship between respiration and heart rate oscillations (13, 44). RSA serves an important role in providing synchrony between the respiratory and cardiovascular systems, together responsible for maintaining the metabolic equilibrium over a wide range of physical and psychological conditions. RSA magnitude is dependent on both respiratory frequency and tidal volume, even when the autonomic tone remains stable (25). A key problem in cardiorespiratory engineering is to efficiently and accurately quantify and monitor RSA under different physiological and behavioral conditions where frequent variations in respiration and autonomic inputs are present. A solution to this issue could yield critical insights into the mechanisms involved in short-term and long-term cardiorespiratory coupling (13, 43).

Accurate quantification of RSA serves several purposes. First, it is an indirect measure of parasympathetic cardiac control, and it has been shown that pharmacologically induced changes in cardiac vagal tone (eg. atropine administration) can be accurately tracked by RSA measures (15, 20). Thus, RSA mainly reflects cardiac vagal efferent effects on the sinoatrial node (14, 23), though at lower respiratory frequencies sympathetic cardiac control can contribute to RSA as well (43). RSA could be used as a stable measure of parasympathetic cardiac control not only in the controlled environments but also in the ambulatory recordings (19). Second, evidence has been given for using RSA as a predictive marker of risk of physiological morbidity (29). Generally, it has been found that lower RSA is associated with high risk of morbidity. Third, RSA has become a central point in evolution theory of central and vagal control of cardiorespiratory interactions (22).

In early work, RSA was defined using simple time domain measures of beat interval series (28, 21). The algorithm measures the RSA as the difference between the shortest Pulse Interval (PI) during inspiration and the longest PI during expiration within one respiratory cycle, and has been thus referred as peak-valley method. The differences are then averaged across several breaths, thus requiring a predetermined estimation window. Furthermore, when using peak-valley based RSA estimations respiratory parameters need to be controlled carefully, which may not be possible under some behavioral conditions (23). In addition to the standard time domain measures (33), new indices have been proposed for the long-term heart rate variability analysis (41, 34). Several other methods have been devised in the last decades, both in the time and frequency domain, successfully relating measures of heart rate variability to RSA and cardiopulmonary coupling (35, 26, 6, 36, 45, 18, 11, 7, 16, 46, 37). Also, combined linear and nonlinear methods have been devised for heart disease screening (24, 48). In particular, filtering and transfer function approaches were also used in quantifying RSA (44, 49), and a bivariate autoregressive model was further proposed to estimate the time-varying RSA gain (3, 2). As most of these methods are not able to completely overcome non-stationarity issues and not capable of estimating the fast changes in RSA at arbitrarily small time scales, a point process framework for heart beat dynamics (1, 4) has been proposed to assess RSA within an adaptive point process filtering algorithm (9). However, the above approaches do not account for irregular respiration dynamics when assessing RSA.

This paper introduces a maximum likelihood point process framework for instantaneous estimation of RSA. This estimation method has been selected for its robust dynamic identification qualities, and because it is less sensitive to parameter initialization and numerically more stable than other adaptive recursive point process filters. The proposed method allows time increments for parameter update at arbitrary small time intervals, thus achieving instantaneous estimations of HRV and RSA measures. As a consequence, rapid RSA changes at time intervals smaller than pulse or respiratory cycles can be monitored and accounted for, and to get consecutive RSA estimates it is not necessary to average over overlapping time windows, or wait till the next respiratory cycle or the next pulse, as required by other existing methods. 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 estimating 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 the PIs. Such a 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, where respiration patterns are considerably altered.

The structure of this paper is as follows: In Section 2 we present the maximum likelihood point process framework, the methods for evaluating model goodness-of-fit, and the methods for estimating instantaneous RSA gain within the framework. Section 3 reports on results obtained from applying the proposed algorithms on simulated signals with varying respiration and autonomic inputs, and on results of a meditation protocol with non-stationary respiration patterns respectively. Finally, in Section 4 we present discussions and final conclusions.

2 Methods

2.1 Point Process Model of Heart Beat Interval Dynamics

Integrate and fire models are regularly used to simulate heart beats, and such models postulate that the resulting times between two firing events (the PIs) have statistical properties of an inverse Gaussian process (4). Additionally, autonomic inputs to the SA node are part of the cardiovascular control circuitry, thus the PI variations are dynamic, or time-varying (4). For this reason, we here model the pulses 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

f(t|uk)=[θ(t)2π(tuk)3]12exp {12θ(t)[tukμPI(t)]2μPI2(t)(tuk)},

where t is at 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. The standard deviation of the above PI probability model is given by (10)


Because of the lasting effects of autonomic inputs to the SA node, μPI (t) in the point process probability model should be modeled as dependent on the recent history of the previous PI intervals,


Thus, the mean value of the distribution is modeled as a univariate p-order autoregressive (AR) process. According to the model, the mean PI interval is influenced by the past p PI intervals, thus dependent on the AR coefficients {ak(t)}k=0p, whereas the PI interval variance is determined by{ak(t)}k=0p and θ(t). All the model parameters are time-varying, thus taking into account the non-stationary behavior of heart beat dynamics, and allowing for definition of instantaneous estimates of heart rate variability (HRV). Later we shall see an extension to bivariate AR in Section 2.4 which takes into account non-stationarities in respiration.

2.2 Heart Rate and Heart Rate Variability

Heart rate is often defined as the reciprocal of the PIs (33). 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 (42), we get heart rate probability density f(r) = f(c/(tuk)) as


Mean and standard deviation of the heart rate probability density are given by (4)



Instantaneous estimates of heart rate and heart rate variability are characterized by μ HR (t), and σHR2(t), respectively.

2.3 Heart Beat Spectral Components

Spectral analysis of HRV has been deemed useful in measuring the sympathovagal balance of a subject. The frequency response for the PI interval series is given by


where fs is the beat rate. we can then evaluate the dynamic power spectrum, or the parametric auto-spectrum (5) by


Main spectral components can be identified from the above auto-spectrum by subdividing into three frequency bands as specified by the current HRV standards (33). The low frequency component (LF, 0.04–0.15Hz) indicates the slow autonomic effects (sympathetic and parasympathetic) on the SA node, whereas the high frequency component (HF, 0.15–0.5Hz) indicates only parasympathetic modulation. The very low frequency component (VLF, 0.003–0.04Hz) is dubious and is rarely attributed to physiological processes, especially in short-term recordings. As a consequence, an instantaneous estimation of LF/HF power ratio would provide important information about the dynamic sympathovagal balance of a subject. Furthermore, the HF band has been defined as dynamically centered around the respiratory frequency (37).

2.4 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,


The original respiration signal (RP) is re-sampled at the pulse timings, so that both respiration and PIs are synchronized. Re-sample of RP can be can be performed quite accurately, as current respiration measuring techniques offer high sampling rates.

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


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


Then, RSA gain can be estimated by evaluating Eq. 10 at ωRP,


Alternatively, we evaluate the RSA gain at the frequency where maximum interaction between PIs and respiration occur. In this regards, we use the time-varying coherence spectrum Coh(ω, t) (49)


where PPI−RP (ω, t) is the cross-spectral density between respiration and the PIs, while PRP (ω, t) and PPI (ω, t) are the auto-spectral densities of respiration and PIs, respectively. The time-varying coherence spectrum is then used to estimate the frequency ωcoh (t) where coherence is maximum, i.e.,


and the RSA gain is evaluated at ωcoh,


Although computed differently, both RSA indices are expected to yield similar results, as the PIs are highly correlated with the breathing cycle, and coincident with the main respiratory frequency in spectral terms. In Section 3.5, we test the validity of this assumption, pointing at both RSAgaina(t)andRSAgainb(t) as reliable estimates of instantaneous RSA.

2.5 Local Maximum Likelihood Estimation

A local maximum likelihood method (32, 47) is implemented to estimate the unknown time-varying parameter set ξ={{ak}k=0p,{bk}k=1q,θ}. In estimating ξ at time t, we take a local likelihood interval (tl, t], where l is the length of the local likelihood observation window. 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}. The log-likelihood associated with the joint probability density is given by

log f(utl:t)=j=2nw(tuj) log f(ujuj1)+w(tun) log tunf(ν)dν,

where w(tuj) = αtuj, 0 < α < 1, is a weighting function for the local likelihood estimation (4, 32). The weighting time constant α governs the degree of influence of a previous event observation uj on the local likelihood at time t. The second term of Eq. 16 represents the likelihood of the partially observed interval since the last observed pulse un (right censoring). To maximize the local log likelihood given in Eq. 16 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 as arbitrarily small, thus yielding instantaneous estimates of heart rate and heart rate variability.

2.6 Model Goodness-of-Fit

Goodness-of-fit of the proposed model is evaluated using a Kolmogorov-Smirnov (KS) test based on the time-rescaling theorem (8). The test uses the conditional intensity function λ(t)=f(t)/[1untf(ν)dν] 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. If the model is correct the KS plots should align with the 45 degrees diagonal (within 95% confidence intervals), and the autocorrelation values should be as close as possible to zero (also within 95% confidence intervals). The KS distance is defined as the maximum distance of the KS plot from the diagonal.

3 Results

3.1 Simulations

Algorithms were tested with simulated signals for both constant and dynamic respiratory conditions in order to evaluate the robustness of the proposed methods under altered respiratory patterns. A simple sinusoidal model was used for both PI and RP signals as given below.

PI(t)=μPI+α(t) sin(2πf1(t)t)+β(t) sin(2πf2(t)t)+n1(t)

RP(t)=γ(t) sin(2πf2(t)t)+n2(t)

The LF and HF components of the simulated PIs are modeled by f1(t) and f2(t) with amplitude factors α(t) and β(t) respectively. With the basic assumption that the HF component is coherent with respiration, the RP frequency will be given by f2(t). The amplitude of the RP waveform is modeled by γ(t). n1(t), and n2(t) are additive white gaussian noises. Such formulation is capable of simulating dynamic autonomic inputs via setting the α(t) and β(t) parameters, thus varying the LF and HF powers. The LF/HF ratio, which is proportional to (α(t)/β(t))2, can consequently also be varied. The model is also capable of altering the RSA gain by changing the ratio β (t)/γ(t). The proposed point process algorithm is used to perform two sets of simulations which estimate the dynamic RSA gain at both constant and dynamic respiration by varying β(t), and then β(t) and f2(t), respectively.

3.1.1 Constant Respiration Frequency

The first set of simulations was devised in order to show that the proposed algorithm is capable of accurately estimating the RSA gain while respiratory frequency remains relatively constant, and the results are shown in Figure 1a. Parameters in Eq. 17 and Eq. 18 were set to μPI = 1s, f1(t) = 0.1Hz, f2(t) = 0.3Hz, α(t) = 0.1s, γ(t) = 0.1rpu, β(t) has a step change from 0.2s to 0.1s at 500s, and then a linear increment from 0.1s to 0.2s between 1000s and 1500s. Additive white Gaussian noise was used for n1(t) and n2(t) with signal to noise ratio set at 20dB. A bivariate AR model (order 6) was used to estimate the time-varying respiratory frequency and the coherence between the PIs and respiration at the beats. The regression was defined at the beats, thus neither interpolation nor resampling was necessary. Results demonstrate the ability of the algorithm to track step changes as well as gradual changes in RSA gain while respiratory frequency remains constant. Respiratory frequency was accurately estimated to be 0.3Hz, and a very high coherence was observed (close to 1) between PI and RP at that frequency as seen from Figure 1a. The step change at 500s was tracked accurately, and the point process algorithm was capable of reaching 95% of the lower RSA level within 41s. Gradual changes in RSA from 1000s to 1500s were also followed accurately.

Fig. 1
fRP (top), coherence (middle), and RSA gain (bottom) as estimated by applying the point process algorithm to simulated PI and RP signals with constant (left) and varying (right) respiratory frequency (dotted curve shows the theoretically expected RSA ...

3.1.2 Dynamic Respiration Frequency

The second set of simulations was done in order to show that the proposed algorithm is robust under sudden and gradual changes in the respiratory frequency, and still capable of accurately estimating the RSA gain. The simulation results for the dynamic respiration is shown in Figure 1b. In this case, respiratory frequency f2(t) was linearly increased from 0.25Hz to 0.35Hz between 0sand 500s, suddenly dropping back to 0.25Hz at 1000s. All the other parameters were set at same values as in Section 3.1.1. Figure 1b shows the simulation results, where we observe a drop in coherence down to 0.7 due to the step change in respiratory frequency at 1000s. More precisely, the coherence drop occurs at 1026s, 26s after the actual frequency drop due the transient effect. The estimated RSA gain during the gradual increase in respiratory frequency (from 0s to 500s) fluctuates just below 2, which demonstrates the robustness of the point process algorithm under altered respiratory conditions. Note that the step change in the respiratory frequency at 1000s does not have an impact on accurate estimation of RSA gain during that period.

3.2 Meditation Protocol

Numerous reports have documented phase-locked decreases in respiratory rate during meditation periods (31, 39), thus this is an ideal protocol for evaluating the robustness of the proposed methods under time-varying respiration activity concurrent with dynamic autonomic inputs. The data were acquired from 13 experienced practitioners of Buddhist Insight (a.k.a. Vipassana) meditation, which is a form of ‘mindfulness’ meditation. Mindfulness is defined as “paying attention in a particular way: on purpose, in the present moment, and non-judgmentally” (27, 40). The study participants were instructed to perform a breath awareness meditation technique, which consists of focusing attention on breathing sensations (flow of air through the nose, or rise and fall of the abdomen) and passively ignoring everyday thoughts. In total 4 females and 9 males, (25 to 49 years old, average 38.4), were included. They had been practicing Insight meditation for between 1 and 20 years (average 8.7), and were required to have been practicing 40 minutes per day at least 5 days per week for at least 1 year and to have attended at least one week-long meditation retreat. None reported taking any medication or having any cardiovascular disease. These subjects are compared to a demographically matched control group of 10 subjects, 4 females and 6 males, (24 to 49 years old, average 35.7), with no previous yoga or meditation experience.

The experiment which was conducted inside a MRI scanner as part of a larger study (30), started with a 6 minute of baseline period, followed by 1 minute of fixation, then 24 minutes of meditation, followed by 1 minute of fixation, and finally 6 minutes of silent random number generation. The control group did not meditate, but rather simply rested throughout the corresponding 24-minute period. For analytical purposes the meditation session is divided into three epochs (early, middle, late) of 8 minutes each, referring to the temporal stage of the meditation. During the experiment, the blood volume pressure (BVP) signal was recorded, and used to identify the PIs. Additionally, two respiration signals, a flow signal proportional to the airflow changes of the breath, and a belt signal proportional to the lung volume changes were recorded. All the raw signals obtained were initially sampled at 1kHz. Respiration signals were resampled at the pulse timings to obtain the respiration values at the beats, after low pass filtering (cutoff 0.5Hz) to avoid aliasing.

A subset of nine meditation subjects further performed an experiment in which they were cued to inhale and exhale in the exact pattern as during a previous meditation period (paced breathing), but without employing any meditation techniques. This was done in order to evaluate the intrinsic effects of meditation as opposed to the effects of ‘meditative’ breathing. External visual stimuli were used to cue the subjects to reproduce the same breathing patterns as in a previous meditation session. Both protocols were approved by the Massachusetts General Hospital Institutional Review Board. Written informed consent was obtained from all study participants and they were compensated for completion of assessment procedures.

3.3 Instantaneous PI and HRV Estimation

PI and RP series from each recording session were used to estimate the optimal model parameters for the proposed point process model. Optimal values for regression orders of the bivariate model p and q, local 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 proposed maximum likelihood point process model was then applied to both meditating and control subjects with these optimal model parameters. The instantaneous mean PI estimate μPI (t), PI variance σPI2(t), mean heart rate μHR(t), and HRV index σHR2(t) are updated every Δ= 5ms.

The estimated instantaneous HR and HRV indices for two representative subjects, one from each group, are shown in Figure 2. In particular, sudden variations can be observed in the estimated mean PI of the experienced meditator. 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 sympathovagal balance. Sharp increases in the HRV index are consequently observed at corresponding timings. Of note, similar transients were less frequent in the mean PI of the control subject.

Fig. 2
Point process instantaneous measures of heart rate and heart rate variability. Top panels show instantaneous mean PI, PI variance, mean heart rate, and heart rate variability indices for (a) an expert meditator and (b) a control subject estimated by the ...

The goodness-of-fit of the model was tested against the experimental data using KS plots and transformed quantiles’ autocorrelation functions for each subject, as shown in Figure 2 for two representative subjects. For all 23 subjects considered, the KS plots approximately follow the uniform cdf (the ideal fit), and mostly stay within the 95% confidence interval. The reportedly small KS distances further indicate that the model fits well with the original data. Additionally, autocorrelation plots were also generated for up to 60 lags (approximately 1 minute), and we observed a low autocorrelation for all the subjects. Almost all the points stay within the 95% confidence interval, which indicates the model has successfully extracted all the dynamic structure from the original data.

3.4 Time-Frequency Analysis of PI and RP

After validating the model, time-frequency analysis of PI, RP, and coherence between the two was performed for each subject in order to identify the dynamic respiration frequency ωRP(t) and the dynamic frequency where maximum coherence occurs ωcoh(t), as a preliminary step of evaluating dynamic RSA. The time-varying bivariate AR model (order 6) was used to estimate the dynamic auto-spectra of PI and RP, as well as the cross-spectrum between the two, thus generating the time-frequency plots shown in Figure 3 for the two representative subjects. From the RP and coherence time-varying spectral distributions, it is evident that the frequency of maximum coherence (ωcoh(t)) closely matches with the respiratory frequency (ωRP (t)), confirming the mentioned assumption that the highest correlation between PIs and RP activity generally occurs at the predominant respiratory frequency.

Fig. 3
Time-frequency distributions of PI power (top), RP power (middle), and Coherence between PI and RP (bottom) for the expert meditator (left panel) and the control subject (right panel).

For the meditation subject, a notable drop in the predominant respiratory frequency is observed during meditation. During baseline and number generation phases the subject breaths at around 0.25Hz, whereas respiratory frequency falls to around 0.2Hz during meditation (even slower in early meditation). Of note, each expert meditator exhibits different trends in respiration, likely due to individual differences in the ability to sustain a deep meditative state. Therefore, a reliable method for autonomic control assessment should be able to account for a wide range of respiratory dynamic changes. As the proposed method relies on estimates of instantaneous respiratory frequency to compute the RSA gain, it is possible to follow such variations in breathing. As we shall see in later analysis, measures based on the standard subdivision of frequency bands are not able to account for such changes.

3.5 Instantaneous RSA Estimation

The bivariate model in Eq. 9, together with the new time-frequency spectral characterization, was used to estimate the instantaneous RSA. As applied to the considered experimental data, the first important finding was that both RSAgaina(t) (respiration based) and RSAgainb(t) (coherence based) gave very similar RSA estimates in all the subjects considered. As previously noted, this finding was somehow expected, as PIs are often highly correlated with respiration activity at the main respiratory frequency, indicating both RSAgaina(t)andRSAgainb(t) as reliable estimates of instantaneous RSA. Results are shown in Figure 4 for one experienced meditator, and for a control subject. 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.

Fig. 4
Meditation subject (left) and control subject (right): From top to bottom: dynamic frequencies where maximum RP power is observed (fRP, blue) and coherence between PIs and RP is maximum (fcoh, red); coherence at fRP (blue) and maximum coherence (red); ...

Second, a significant increase in RSA gain is clearly noticeable during meditation, whereas in the control subject RSA gain seems to fluctuate around baseline levels along the entire experiment1. As the RSA transfer function exhibits low-pass characteristics (44), it could be argued that increases in RSA during meditation are not a result of meditation itself, but are only a result of the drop in breathing frequency. In order to account for such argument, we estimated the RSA gain at a fixed frequency, and Figure 5a shows the comparison between two RSA estimates, one at the dynamic respiratory frequency, and one at a fixed frequency all along the experiment (ffix = 0.223Hz for the particular meditation subject). The fixed frequency ffix was chosen so that, on average, the highest RP power was observed at ffix during the entire experiment. We refer to this fixed frequency based RSA estimation as static evaluation as opposed to the dynamic evaluations, RSAgaina(t) and RSAgainb(t), described in the proposed methods. Higher RSA was found during meditation regardless of whether the static or dynamic evaluation method was used, thus indicating that meditation yields higher RSA independently of changes in respiration. For the RSA estimation based on static frequency, drops in coherence were observed more frequently, especially during the early meditation epoch. Thus, the proposed dynamic RSA estimation based on time-varying respiration relies on higher values of coherence to provide more accurate results on the linear interactions associated with cardiorespiratory coupling. It is important to highlight that the static evaluation requires observing the entire data segments in order to estimate the mean respiratory frequency, and as such it could not be implemented in an on-line fashion for prospective real-time monitoring use. Furthermore, our dynamic evaluations add in flexibility, and overcome at the same time the problems associated with a rigid subdivision in frequency bands (33). A more comprehensive statistical analysis as applied to the two groups considered is presented in Section 3.7.

Fig. 5
(a) Comparison of RSA estimations based on dynamic (blue) and fixed (red) frequency for the expert meditator subject, (b) Comparison of RSA estimations obtained using flow RP (blue) and belt RP (red) signals for the expert meditator subject.

Third, we compare two RSA estimates obtained by using two distinct acquisition methods for the respiration signals. The first method measured the air flow through the subject’s airways (flow RP), whereas the second method measured lung volume from the expansion and contraction of the abdomen using a respiration belt (belt RP). Even though both signals give similar respiratory frequencies, the amplitudes are incomparable, thus they were normalized by the standard deviation of each RP signal in each single experiment. Figure 5b shows a comparison of the RSA instantaneous estimations obtained using flow and belt respiration signals for the previously shown expert meditator subject. The correlation index between RP frequencies estimated for two signals was 0.9593, whereas the corresponding RSA estimates had a correlation of 0.8253. Such high correlation indicates that regardless of the respiration measuring method, both flow and belt RP signals lead to a similar and consistent estimation of RSA. Deviations from a correlation of 1.0 can be reasonably attributed to experimental noise.

3.6 LF-HF Spectral Analysis

We here compare the standard frequency domain indices (33) with the proposed RSA measures. Figure 6 shows the PI power in the LF and HF bands, the LF/HF ratio, and the percentage of RP power in the LF band for an exemplary meditation subject. It is interesting to note that for the considered subject, the average respiratory frequency was found to be 0.146Hz, which is at the borderline between the LF and HF bands. As shown in Figure 6, on average 47.3% of RP power falls in the LF band. This means that the respiratory-driven autonomic control on HRV is now present in both frequency bands, thus causing a bias in the LF/HF ratio as an indicator of the sympathovagal balance.

Fig. 6
From top to bottom: PI power in LF-band 0.04–0.15Hz; PI power in HF-band 0.15–0.5Hz; LF/HF ratio; and percentage RP-power in LF-band 0.04–0.15Hz for the expert meditator subject.

A slight increase in HF is observed during the middle and late meditation epochs possibly due to the effects of meditation, but those changes are also mirrored in LF. Therefore the LF/HF ratio fails to capture the corresponding autonomic dynamics happening during meditation sessions where slow breathing cycles are usually observed. In contrast (Figure 4), the proposed RSA estimates show an increase in RSA during practice epochs as a result of meditation influence, as would be expected from self-reports of increasing relaxation over the course of a typical mediation session.

3.7 Statistical Analysis

In this section we present a more detailed statistical comparison of both meditation and control groups in terms of the proposed RSA indices as well as a set of standard HRV indices. A summary of average pulse intervals, heart rates, their variances, and standard SDNN and pNN50 estimates for the two groups is shown in Table 1. On average, the control group exhibits a lower heart rate when compared to the expert meditator group, although with no statistically significant differences between the two groups. Of note, in Figure 2, the control subject exhibits a higher heart rate compared to the expert meditator, which is in disagreement with the mean statistics, possibly due to large variations in heart rate within each group. For both groups no significant changes in mean heart rate were observed during practice or number generation epochs compared to the baseline. Standard time domain measures SDNN and pNN50 as defined in (33) also fail to show statistical significance between the two groups during practice epochs, although on average both SDNN and pNN50 are somewhat higher for the control group. For both groups a general increase in SDNN is observed during practice epochs compared to baseline or number generation epochs. In summary, standard time-domain measures of HRV are not able to report any significant changes distinguishing meditation practice from very relaxing states.

Table 1
Statistical comparison of point process heart beat indices between meditation and control groups, along with the standard time domain indices (33) (mean ± sd).

Average RSA statistics estimated by the proposed point process algorithm along with respiration and coherence statistics are shown in Table 2. We also compare the proposed RSA estimations with standard frequency domain HRV measures and RSA measures computed within the HF-band. As shown in Table 2 PI power in HF has increased (no statistical significance) in expert meditators during meditation practice, suggesting an improved vagal tone during meditation. We do not observe statistically significant changes in LF/HF as well, possibly due to the considerable transient presence of respiration within the LF-band during meditation practice. Conversely, RSA gain estimated at instantaneous respiratory frequency using the proposed point process model does reveal significant differences between the baseline and meditation condition for the meditators. In fact, all three meditation epochs show a significant increase in RSA (p < 0.05), where increment is between 29% and 44% compared to the baseline. No such difference was found in the control group. In addition, coherence between PIs and RP remains high, between 0.86 and 0.96, at the respiratory frequency for both groups.

Table 2
Statistical comparison of RP, Coherence, and proposed RSA indices between meditation and control groups along with standard frequency domain HRV measures (mean ± sd).

Within the meditation group, meditation epochs result in decrease in respiratory frequency as well as respiratory power, though changes are not significant. Previously it has been reported that the change in respiration rate during meditation compared to baseline is highly correlated with meditation practice (30), suggesting that changes in respiration reflect increasing mastery of the meditation technique. Respiratory rate increases slightly during number generation phase for both groups compared to baseline, also without statistically significant differences.

In order to show that RSA statistics are not biased by lower predominant respiratory frequencies during meditation, we have computed the RSA gains at fixed frequencies where (on average) maximum respiration is observed. We still found mean values comparable to the dynamic assessment (i.e., higher RSA during meditation which is not evident in control subjects), although with higher standard error, lower significance levels, and lower coherence values at constant frequencies. This important outcome indicates that RSA gains estimated accounting for dynamic respiratory frequencies provide a more accurate linear assessment, and consequently yield lower variability of the estimates and less identification uncertainty. During the silent random number generation phase, RSA values are statistically comparable to baseline levels in both groups.

Previous frequency domain methods estimated RSA within the HF-band to disentangle effects from the baroreflex loop (2). The proposed method is independent of such a frequency division, as it is based on time-varying instantaneous respiratory parameters. In order to explore frequency band biases, we have applied our proposed method singularly within the HF-band and LF-band (also in Table 2). Mean respiration power in HF decreases during meditation due to the shift in breathing rate to lower frequencies, then increasing back to baseline levels in the random number generation epoch. A similar trend is noticeable in the control group as well, but it is less significant. As expected, we observe an increase in RSA during meditation, but all the increments are less significant compared to the proposed dynamic RSA estimation where the entire spectrum was accounted for in the analysis. Interestingly, HF-band based estimation fails to show a significant RSA increase (p > 0.1) for the third meditation epoch, compared to the highly significant results given by our dynamic RSA estimation (p = 0.005).

In the LF-band, respiration power increases considerably during meditation. It is interesting to note that high coherence (≥ 0.74) is observed during meditation epochs between respiration and PIs even in this range. This is possibly due to the ‘leakage’ of RP power into LF during meditation, consequently increasing RSA in the LF band within the baroreflex control loop, and biasing the LF/HF ratio as a measure of sympathovagal balance. Importantly, results from the control group do not evidence a similar effect.

Figure 7 shows the mean RSA as estimated by the proposed point process model for the two groups along different epochs in the experiment. Again, ‘practice’ epochs implicate engaging in meditation for the meditator group, while ‘practice’ refers to relaxation for the control group. The expert meditator group has a significantly high RSA gain during meditation, while the control group does not show significant changes compared to the baseline epoch. Expectedly, both dynamic RSA measures proposed in this paper, based on respiration and coherence, give similar RSA estimates statistics.

Fig. 7
Mean and standard deviation of dynamic RSA estimates averaged for both meditator and control groups.

Finally, we further compare the RSA statistics of the meditation group during the paced breathing experiment with the meditation experiment. As expected, paced breathing successfully reproduces respiration characteristics as during meditation, with undistinguished mean respiratory frequencies (fRP = 0.17 ± 0.06Hz paced vs. 0.17 ± 0.04Hz meditation). As we compare results from the respective respiration dynamics within each subject, it is reasonable to assume that the differences between the two experiments solely represent the effects of meditation. On average, at dynamic respiratory frequency, mean coherence was found to be 0.91 ± 0.07 (0.89 ± 0.09 with meditation), and at the fixed frequency it dropped to 0.80 ± 0.08 (0.75 ± 0.14 with meditation). Overall, even in the presence of slower respiratory frequencies, RSA estimation did not alter significantly during paced breathing as compared to the baseline. For the group of subjects considered, sign-rank test between baseline and paced breathing showed no significance (p = 0.10) as opposed to the high significance (p = 0.0005) obtained during meditation. This result further supports the finding that changes in RSA during meditation are not simply due to changes in respiratory rate, but effectively associated to meditation practice.

4 Discussion

We have introduced a point process model for instantaneous estimation of a set of HRV and RSA measures under dynamic respiration characteristics. Within the point process model, the mean pulse interval is modeled as a bivariate regression of both previous PIs and respiration values. We termed the feed-forward dependency of the PIs on respiration as RSA, and a transfer function approach was used to quantify the instantaneous RSA estimates. The model was validated statistically using Kolmogorov-Smirnov (KS) goodness-of-fit tests and independence tests.

Obviating any artificial subdivision in frequency bands, our proposed method is robust to variations in respiration patterns and relies on the instantaneous respiratory spectral power distribution and the dynamic coherence between respiration and PIs. As a further proof that RSA gains are not considerably biased by possible shifts in RP predominant frequency, we have computed similar indices at pre-estimated fixed frequencies where coherence is the highest on average along each experiment, and verified that average statistics using fixed frequency based estimates corroborate the results from the proposed method. We have also demonstrated that considering the entire spectrum instead of designated frequency intervals leads to more reliable and less variable estimates.

Statistical analysis further reports a significant increase of RSA during meditation for the expert meditator group, while no relevant difference and no statistical significance was found for the control group, or during paced breathing. While traditional HRV indices confirm previous findings pointing at an autonomic assessment not necessarily reflective of quiescent cardiac dynamics (39), and possibly biased by altered breathing patterns (31), our statistical RSA appraisal substantiates previous reports of a marked increase of RSA during meditation (12), and further indicates that these changes are not simply due to changes in respiration. Although unlikely, it is possible that mental effort of following the visual cues during the paced breathing condition may have contributed to the observed differences in RSA between this condition and the meditative state. However, the meditative state also requires some mental effort to remain focused, and so differences in cognitive load between the conditions is likely minimal.

Our findings are further substantiated by the absence of an increase in RSA during an experiment in which respiratory dynamics similar to those in meditation were reproduced without effective meditation practice. Comparison between the two protocols allows for isolation of the physiological effects of meditation, thus removing possible additional biases in parameter estimation due to systematic changes in respiration characteristics, and it helps in identifying important changes in cardiorespiratory coupling in the presence of dynamic breathing patterns. Therefore, our overall results support the hypothesis that increased RSA during mindfulness meditation performed by experienced practitioners may reflect adaptive properties of the central-autonomic nervous system to this practice, that are not related to any ‘voluntary intention’ towards the breathing process.

Future research will be required to identify physiological and psychological causes of increased RSA during meditation. As some studies suggest, meditation-based clinical interventions can impact physiological recovery during emotional challenges (17, 38). These RSA measures may provide useful novel applications in dynamic instantaneous impact of meditation and other behavioral interventions on autonomic functioning. Future studies will address whether these methods can also be used to determine an objective measure of meditative states. Such research may lead to finding how meditation practice is capable of positively influencing mood and decreasing arousal with targeted changes to the vagal tone.

Future improvements to the model will consider introducing a vagal filter to the proposed model, thus obtaining more direct estimates of vagal control. Also, devising a multivariate model including arterial blood pressure as further covariate and accounting for directionality and causality will enable to disentangle the closed-loop blood pressure control dynamics from the respiratory-driven effects. Additionally, consideration of nonlinear models to comprise the nonlinear interactions observed in cardiovascular control and cardiorespiratory coupling may also lead to more accurate instantaneous RSA estimates. A precise dynamic RSA assessment could provide a useful complement to current standard HRV measures with the prospect of including RSA quantification into standard clinical or ambulatory monitoring systems to achieve better patient care.

We have proposed a maximum likelihood point process model for instantaneous RSA estimation combined with a time-frequency RSA evaluation based on pulse intervals and respiration spectra 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 or paced breathing experiment, providing interesting insights into the effects of meditation on cardiorespiratory function and 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 real-time 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, Grant K01-AT00694-01.


1Note that for a fair comparison across subjects, the RSA gain was normalized by the standard deviation of the corresponding RP.

Contributor Information

Sandun Kodituwakku, Applied Signal Processing Group, School of Engineering, The Australian National University, Canberra, Australia. Neuroscience Statistics Research Laboratory, Massachusetts General Hospital, Harvard Medical School, Boston, MA 02114 USA, ua.ude.una@ukkawutidok.nudnas, Tel.: + 61-2-612-58689.

Sara W Lazar, Department of Psychiatry, Massachusetts General Hospital, Harvard Medical School, Boston, MA 02114 USA, ude.dravrah.hgm.rmn@razal, Tel.: + 1-617-724-7108.

Premananda Indic, Department of Neurology, University of Massachusetts Medical School, Worcester, MA 01655 USA, ude.demssamu@cidni.adnanamerp, Tel.: + 1-617 287 6050.

Zhe Chen, Neuroscience Statistics Research Laboratory, Massachusetts General Hospital, Harvard Medical School, Boston, MA 02114 USA. Harvard-MIT Division of Health Science and Technology, Massachusetts Institute of Technology, Cambridge, MA 02139 USA, ude.tim@nehcehz, Tel.: + 1-617-324-1882.

Emery N Brown, Neuroscience Statistics Research Laboratory, Massachusetts General Hospital, Harvard Medical School, Boston, MA 02114 USA. Harvard-MIT Division of Health Science and Technology, Massachusetts Institute of Technology, Cambridge, MA 02139 USA. ude.tim.tatsoruen@bne, Tel.: + 1-617-726-7487.

Riccardo Barbieri, Neuroscience Statistics Research Laboratory, Massachusetts General Hospital, Harvard Medical School, Boston, MA 02114 USA. Harvard-MIT Division of Health Science and Technology, Massachusetts Institute of Technology, Cambridge, MA 02139 USA, ude.tim.tatsoruen@ireibrab, Tel.: + 1-617-724-1061.


1. Barbieri R, Brown EN. Analysis of heartbeat dynamics by point process adaptive filtering. IEEE Trans Biomed Eng. 2006;53(1):4–12. [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 bivariate autoregressive techniques. Ann Noninvasive Electrocardiol. 1996;1:264–277.
3. 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]
4. Barbieri R, Matten EC, Alabi ARA, Brown EN. A point-process 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]
5. Baselli G, Bolis D, Cerutti S, Freschi C. Autoregressive modeling and power spectral estimate of RR interval time series in arrhythmic patients. Comput Biomed Res. 1985;18(6):510–530. [PubMed]
6. Berntson GG, Bigger JT, Eckberg DL, Grossman P, et al. Heart rate variability: origins, methods, and interpretive caveats. Psychophysiology. 1997;34(6):623–648. [PubMed]
7. Blain G, Meste O, Bermon S. Influences of breathing patterns on respiratory sinus arrhythmia in humans during exercise. American Journal of Physiology-Heart and Circulatory Physiology. 2005;288(2):H887. [PubMed]
8. 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]
9. 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]
10. Chhikara RS, Folks JL. The inverse gaussian distribution: theory, methodology, and applications. New York, NY, USA: Marcel Dekker, Inc.; 1989.
11. Denver JW, Reed SF, Porges SW. Methodological issues in the quantification of respiratory sinus arrhythmia. Biological Psychology. 2007;74(2):286–294. [PMC free article] [PubMed]
12. 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]
13. Eckberg DL. Human sinus arrhythmia as an index of vagal cardiac outflow. J Appl Physiol. 1983;54:961–966. [PubMed]
14. Eckberg DL. The human respiratory gate. Journal of Physiology. 2003;548(2):339–352. [PubMed]
15. Elstad M, Toska K, Chon KH, Raeder EA, Cohen RJ. Respiratory sinus arrhythmia: opposite effects on systolic and mean arterial pressure in supine humans. Journal of Physiology. 2001;536(1):251–259. [PubMed]
16. 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. Biological cybernetics. 2004;90(6):390–399. [PubMed]
17. Garland EL, Gaylord SA, Boettiger CA, Howard MO. Mindfulness Training Modifies Cognitive, Affective, and Physiological Mechanisms Implicated in Alcohol Dependence: Results of a Randomized Controlled Pilot Trial. J Psychoactive Drugs. 2010;42(2):177–192. [PMC free article] [PubMed]
18. Gilad O, Swenne CA, Davrath LR, Akselrod S. Phase-averaged characterization of respiratory sinus arrhythmia pattern. Am J Physio - Heart and Circul Physio. 2005;288(2):H504. [PubMed]
19. Goedhart AD, Van Der Sluis S, Houtveen JH, Willemsen G, De Geus EJC. Comparison of time and frequency domain measures of RSA in ambulatory recordings. Psychophysiology. 2007;44(2):203–215. [PubMed]
20. Grossman P, Kollai M. Respiratory sinus arrhythmia, cardiac vagal tone, and respiration: Within- and between-individual relations. Psychophysiology. 1993;30(5):486–495. [PubMed]
21. Grossman P, Svebak S. Respiratory sinus arrhythmia as an index of parasympathetic cardiac control during active coping. Psychophysiology. 1987;24(2):228–235. [PubMed]
22. Grossman P, Taylor EW. Toward understanding respiratory sinus arrhythmia: Relations to cardiac vagal tone, evolution and biobehavioral functions. Biological Psychology. 2007;74(2):263–285. [PubMed]
23. Grossman P, Karemaker J, Wieling W. Prediction of tonic parasympathetic cardiac control using respiratory sinus arrhythmia: The need for respiratory control. Psychophysiology. 1991;28(2):201–216. [PubMed]
24. Heitmann A, Huebner T, Schroeder R, Perz S, Voss A. Multivariate short-term heart rate variability: a pre-diagnostic tool for screening heart disease. Medical and Biological Engineering and Computing. 2011:1–10. [PubMed]
25. Hirsch JA, Bishop B. Respiratory sinus arrhythmia in humans: How breathing pattern modulates heart rate. A mer J Physiol. 1981;241:H620–H629. [PubMed]
26. Jo JA, Blasi A, Valladares E, Juarez R, Baydur A, Khoo MCK. Determinants of heart rate variability in obstructive sleep apnea syndrome during wakefulness and sleep. A merican Journal of Physiology-Heart and Circulatory Physiology. 2005;288(3):H1103. [PubMed]
27. Kabat-Zinn J. Full catastrophe living: Using the wisdom of your body and mind to face stress, pain, and illness. Delta. 1990
28. Katona PG, Jih F. Respiratory sinus arrhythmia: noninvasive measure of parasympathetic cardiac control. Journal of Applied Physiology. 1975;39(5):801. [PubMed]
29. Kluge KA, Harper RM, Schechtman VL, Wilson AJ, Hoffman HJ, Southall DP. Spectral analysis assessment of respiratory sinus arrhythmia in normal infants and infants who subsequently died of sudden infant death syndrome. Pediatric research. 1988;24(6):677–682. [PubMed]
30. Lazar SW, Kerr CE, Wasserman RH, Gray JR, Greve DN, Treadway MT, McGarvey M, Quinn BT, Dusek JA, Benson H, L RS, I MC, Fischl B. Meditation experience is associated with increased cortical thickness. Neuroreport. 2005;16(17):1893–1897. [PMC free article] [PubMed]
31. Lehrer P, Sasaki Y, Saito Y. Zazen and cardiac variability. Psychosomatic Medicine. 1999;61(6):812–821. [PubMed]
32. Loader C. Local regression and likelihood. Springer; 1999.
33. 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.
34. Melillo P, Fusco R, Sansone M, Bracale M, Pecchia L. Discrimination power of long-term heart rate variability measures for chronic heart failure detection. Medical and Biological Engineering and Computing. 2011;49:67–74. [PubMed]
35. Mukkamala R, Mathias JM, Mullen TJ, Cohen RJ, Freeman R. System identification of closed-loop cardiovascular control mechanisms: diabetic autonomic neuropathy. American Journal of Physiology-Regulatory, Integrative and Comparative Physiology. 1999;276(3):R905. [PubMed]
36. Mulder LJM. Measurement and analysis methods of heart rate and respiration for use in applied environments. Biological Psychology. 1992;34(2–3):205–236. [PubMed]
37. Orini M, Bailón R, Enk R, Koelsch S, Mainardi L, Laguna P. A method for continuously assessing the autonomic response to music-induced emotions through HRV analysis. Medical and Biological Engineering and Computing. 2010;48(5):423–433. [PubMed]
38. Ortner CNM, Kilner S, Zelazo PD. Mindfulness meditation and reduced emotional interference on a cognitive task. Motivation and Emotion. 2007;31(4):271–283.
39. 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]
40. Peressutti C, Martín-González J, M García-Manso J, Mesa D. Heart rate dynamics in different levels of zen meditation. International Journal of Cardiology. 2010;145(1):142–146. [PubMed]
41. Piskorski J, Guzik P. Asymmetric properties of long-term and total heart rate variability. Medical and Biological Engineering and Computing. 2011:1–9. [PMC free article] [PubMed]
42. Ross SM. Introduction to probability models, 9th edn. Academic Pr; 2007.
43. 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.
44. 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]
45. Sin PYW, Galletly DC, Tzeng YC. Influence of breathing frequency on the pattern of respiratory sinus arrhythmia and blood pressure: old questions revisited. Am J Physio - Heart and Circul Physio. 2010;298(5):H1588. [PubMed]
46. Thomas RJ, Mietus JE, Peng CK, Goldberger AL. An electrocardiogram-based technique to assess cardiopulmonary coupling during sleep. Sleep. 2005;28(9):1151–1161. [PubMed]
47. Tibshirani R, Hastie T. Local likelihood estimation. J Am Stat Assoc. 1987;82(398):559–567.
48. Voss A, Schulz S, Baer KJ. Linear and nonlinear analysis of autonomic regulation of heart rate variability in healthy first-degree relatives of patients with schizophrenia. 2010:5395–5398. [PubMed]
49. Womack BF. The analysis of respiratory sinus arrhythmia using spectral analysis and digital filtering. IEEE Trans Biomed Eng. 1971;18:399–409. [PubMed]