Cardiopulmonary Resonance Function and Indices—A Quantitative Measurement for Respiratory Sinus Arrhythmia

Respiratory sinus arrhythmia (RSA) represents a physiological phenomenon of cardiopulmonary interaction. It is known as a measure of efficiency of the circulation system, as well as a biomarker of cardiac vagal and well-being. In this article, RSA is modeled as modulation of heart rate by respiration in an interactive cardiopulmonary system with the most effective system state of resonance. By mathematically modeling of this modulation, we propose a quantitative measurement for RSA referred to as “Cardiopulmonary Resonance Function (CRF) and Cardiopulmonary Resonance Indices (CRI),” which are derived by disentanglement of the RR-intervals series into respiratory-modulation component, R-HRV, and the rest, NR-HRV using spectral G-causality. Evaluation of CRI performance in quantifying RSA has been conducted in the scenarios of paced breathing and in the different sleep stages. The preliminary experimental results have shown superior representation ability of CRF and CRI compared to Heart Rate Variability (HRV) and Cardiopulmonary Coupling index (CPC).

Introduction

There is an urgent need for quantitative assessment of autonomic nervous function. Heart rate variability (HRV) is widely used as a non-invasive method. Particularly, low-frequency (LF, 0.04–0.15 Hz) and high-frequency (HF, 0.15–0.4 Hz) spectral components of HRV are used as the separate metrics of sympathetic and vagal(parasympathetic) functions (Appel et al., 1989). But as a simplified framework, HRV lacks solid physiological foundation, is not able to accommodate varieties of clinical cases (Hayano and Yuda, 2019). For example, HRV measures will change significantly in different physiologic states such as wake and sleep, exercise and rest, circadian rhythms, as well as with pathologic conditions (Task Force of the European Society of Cardiology the North American Society of Pacing Electrophysiology, 1996).

Cardiopulmonary interaction plays important role in the circulation system, and physiologically presents as respiratory sinus arrhythmia (RSA) phenomenon. RSA is regarded as a non-invasive measure of parasympathetic cardiac control (Katona and Jih, 1975; Topcu et al., 2018). The vagal origin of RSA can be found in the vagal synapses, which are faster than the sympathetic ones and are therefore able to translate central respiratory oscillations present in the brainstem to changes of cardiac sinus node discharge rate, which is not capable for slow sympathetic synapses. Tracking the autonomic regulation in RSA using the electrocardiogram and respiratory measurements is a feasible and important approach to gain our knowledge toward autonomic nervous system and its clinical applications.

The quantitative study of RSA has profound significance in physiology and pathology, as well as extensive clinical applications. Some studies have shown that RSA reaches a relatively stable state in deep sleep (Bernston et al., 1997) and a study of the hibernation of 37 polar bears in the University of Minnesota found that RSA reached their peak during the hibernation. At the same time, as a vagal inflammatory reflex was discovered (Tracey, 2002), quantification of the HRV components, which are not directly related to respiration, is important for the analysis of long-range and scaling properties of the cardiac dynamics (Ivanov et al., 1999; Schmitt and Ivanov, 2007). Examples of application of RSA analysis include clinical psychology (Wielgus et al., 2016), treatment of substance use disorder (Price and Crowell, 2016), prediction of the course of depression (Panaite et al., 2016), quantification of cardiac vagal tone and its relation to evolutionary and behavioral functions (Grossman and Taylor, 2007), quantification of vagal activity during stress in infants (Ritz et al., 2012), and even in cancer patients (Moser et al., 2006), to name just a few.

A variety of data analysis techniques quantifying RSA have been proposed in the literature, for a discussion of commonly used metrics and their advantages and drawbacks see, e.g. (Lewis et al., 2012). The techniques quantifying RSA can be divided into two categories, the time domain and the frequency domain. In time domain, continuous wavelet transform (WTC) is used for its advantage to analyze transient and non-linear signals. This method demonstrates the dynamic behavior of respiration sinus arrhythmia through the analysis of the WTC between heart rate and respiration signals (Jan et al., 2019). Phase analysis technique could help to disentangling respiratory sinus in heart rate variability records, but the final HRV index obtained by this technique is complex to calculate in time domain, and its physiological significance is not clear (Topcu et al., 2018). There are works which use respiratory and RR sequences to calculate G-causality and system gain as the measure of RSA. Much further work is needed to make these produced measures useful in clinical research and applications (Fonseca et al., 2013). In the frequency domain, HF of HRV indicators quantifies RSA on specific frequency bands, while Cardiopulmonary Coupling (CPC) measures the correlation between RR interval and respiratory sequence in the frequency domain. Both are empirical, without solid theoretical foundation and systematic design, therefore serious clinical applications are not seen so far (Thomas et al., 2005).

In the rest of this article, our contributions in developing quantitative measures for RSA are described as follows:

In section Cardiopulmonary Resonance Model (CRM), we model RSA as modulation of heart rate by respiration in an interactive cardiopulmonary system with the most effective system state of resonance. Mathematically, it is described by bivariate autoregressive model of respiration series and RR intervals, and quantitatively it is assessed by Granger causality function. The whole model is referred to as Cardiopulmonary Resonance Model (CRM).

In section Cardiopulmonary Resonance Indices (CRI), based on the cardiopulmonary resonance concept, and Granger causality function which is referred to as cardiopulmonary resonance function (CRF) after, a set of quantitative measures for RSA is proposed, and named as Cardiopulmonary Resonance Indices (CRI).

In section Applications Scenarios, to show the effectiveness of CRM and CRI, two application scenarios, paced breathing and sleep stage discrimination, are studied. It has been shown that CRF and CRI provide ideal visual interpretation and numerical measures for cardiopulmonary interactions toward resonance state in paced breathing scenario as the paced breathing rate coming down to 0.1 Hz. The same is true as the sleep stage moves to deep sleep.

Cardiopulmonary Resonance Model (CRM)

We are committed to building a cardiopulmonary resonance model for the purpose of quantitative assessment of RSA with hypothesis that cardiopulmonary interaction is important in circulation system to ensure efficient delivery of oxygen and nutrient, and that the efficiency is optimized at the state of cardiopulmonary resonance. Mathematically, we present a bivariate autoregressive model of respiration series and RR intervals, calculate respiratory and non-respiratory related component on RR intervals in the frequency domain using Granger-causality.

Bivariate Autoregressive Model of Respiration Series and RR Intervals

The cardiopulmonary interaction can be interpreted as functional connectivity analysis such as synchrony (Engel and Singer, 2001) and phase coherence (Nunez et al., 2001) and so on. Our model takes direct central respiratory modulation of the parasympathetic cardiac signal as the main mechanism for RSA. A powerful technique for extracting directed functional connectivity from data is Granger causality (G-causality) (Granger, 1969). According to G-causality, X₂causes X₁ if the inclusion of past observations of X₂ reduces the prediction error of X₁ in a linear regression model of X₁ and X₂, as compared to a model which includes only previous observations of X₁.

The change process of RR can be regarded as a Markov process, ignoring other factors affecting heart rate in short term, we described the RR intervals(X₁(t)) and respiration signal(X₂(t)) (both of length T) by a bivariate auto-regressive model:

$\begin{array}{l}\begin{array}{c}{X}_1(t)={\displaystyle \sum _{j=1}^p}{A}_{11,j}{X}_1(t-j)+{\displaystyle \sum _{j=1}^p}{A}_{12,j}{X}_2(t-j)+{\xi }_1(t)\\ X_2(t)={\displaystyle \sum _{j=1}^p}{A}_{21,j}{X}_1(t-j)+{\displaystyle \sum _{j=1}^p}{A}_{22,j}{X}_2(t-j)+{\xi }_2(t)\end{array}\end{array}$

where p is the maximum number of lagged observations included in the model (the model order, p < T). A contains the coefficients of the model, and ξ₁, ξ₂ are the residuals for each time series.

In order to ensure RR intervals in the normal range and without a mutation, we use interpolation as a substitute for points that do not meet the following conditions:

$\begin{array}{l}|RR{I}_i-\overline{RRI}|<1.5*Std(RRI)\\ 0.7*RR{I}_{i-1}<RR{I}_i<1.3*RR{I}_{i-1}\end{array}$

where RRI is the RR intervals, RRI_i−1and RRI_i are adjacent intervals.

For each record around 120 s, under the assumption of stationary property of signals, and for efficiency of computation, we normalize the RR intervals and respiration series to zero-mean and unit variance. The magnitude of RSA can be measured by the log ratio of the prediction error variances for the restricted (R) and unrestricted (U) models:

$\begin{array}{l}{G}_{2\to 1}=ln\frac{var({\xi }_{1R(12)})}{var({\xi }_{1U})}\end{array}$

where ξ_1R(12) is derived from the model omitting the A_{12, j}(for all j) coefficients in the first equation and ξ_1Uis derived from the full model.

The estimation of the model of each record requires as a parameter the number of time-lags (p) to include, i.e., the model order. A principle means to specify the model order is to minimize a criterion that balances the variance accounted for by the model, against the number of coefficients to be estimated. We chose the Akaike information criterion (Akaike, 1974) for n variables in which the ∑ denotes the noise covariance matrix:

$\begin{array}{l}AIC(p)=ln(det(\sum ))+\frac{2p{n}^2}{T}\end{array}$

Spectral G-causality of Respiration Series and RR Intervals

For the dynamics of the cardiopulmonary system are easier to understand and interpret in the frequency domain, we calculate the Spectral G-causality of respiration series and RR intervals.

The Fourier transform of the bivariate auto-regressive model in time domain gives:

$\begin{array}{l}\left(\begin{array}{cc}{A}_{11}(f)& A_{12}(f)\\ A_{21}(f)& A_{22}(f)\end{array}\right)\left(\begin{array}{c}{X}_1(f)\\ X_2(f)\end{array}\right)=\left(\begin{array}{c}{E}_1(f)\\ E_2(f)\end{array}\right)\end{array}$

in which the components of A are

$\begin{array}{l}{A}_{lm}(f)={\delta }_{lm}-{\displaystyle \sum _{j=1}^p}{A}_{lm}(j)e^{(-i2\pi fj)},\\ {\delta }_{lm}=0(l=m),{\delta }_{lm}=1(l\ne m)\end{array}$

E is the Fourier transform of the residual matrix.

For the sake of calculation, we rewrite it as

$\begin{array}{l}\left(\begin{array}{c}{X}_1(f)\\ X_2(f)\end{array}\right)=\left(\begin{array}{cc}{H}_{11}(f)& H_{12}(f)\\ H_{21}(f)& H_{22}(f)\end{array}\right)\left(\begin{array}{c}{E}_1(f)\\ E_2(f)\end{array}\right)\end{array}$

where H is the transfer matrix. The spectral matrix S can now be derived as

$\begin{array}{l}S(f)=\langle X(f)X^{*}(f)\rangle =\langle H(f)\sum H^{*}(f)\rangle \end{array}$

in which the ∑ denotes the noise covariance matrix.

A split of U into sub-processes X and Y includes a decomposition

$\begin{array}{l}S(f)=\left(\begin{array}{cc}{S}_{xx}(f)& S_{xy}(f)\\ S_{yx}(f)& S_{yy}(f)\end{array}\right)\end{array}$

of the cross-power spectral density and a similar decomposition for the transfer function H(f).

Then S_xx(f) is the spectral density of X, which is given by

$\begin{array}{l}{S}_{xx}(f)=H_{xx}(f)\sum xx{H}_{xx}^{*}(f)+2\mathrm{Re}\{H_{xx}(f)\\ \sum xy{H}_{xy}^{*}(f)\}+H_{\text{xy}}(f)\sum yy{H}_{xy}^{*}(f)\end{array}$

Thus, we can get the Spectral G-causality of respiration series and RR intervals:

$\begin{array}{l}{G}_{Y\to X}(f)=ln\left(\frac{\left|S_{xx}(f)\right|}{|S_{xx}(f)-H_{xy}(f){\sum }_{y|x}{H}_{xy}{(f)}^{*}|}\right)\end{array}$

$\begin{array}{l}{\displaystyle \sum y|x}\equiv {\displaystyle \sum yy}-{\displaystyle \sum yx}{\displaystyle {\sum }_{xx}^{-1}{\displaystyle \sum xy}}\end{array}$

where ∑denotes the residual covariance matrix.

For the non-respiratory components, we get:

$\begin{array}{l}{G}_{N-RESP}(f)=1\text{-}{G}_{Y\to X}(f)\end{array}$

Now we have both the measurement of respiratory and the non-respiratory components effects on RR intervals in the frequency domain, R-HRV and NR-HRV, respectively. Here we focus on the spectral G-causality of respiration series and RR intervals, G_Y→X(f). For convenience, we simply write it as G(f), and rename it as Cardiopulmonary Resonance Function (CRF) in the rest of this article.

Cardiopulmonary Resonance Indices (CRI)

With cardiopulmonary resonance function (CRF), we are now able to establish a quantitative measurement for RSA, referred to as Cardiopulmonary Resonance Indices (CRI), with the hope that it will be able to play a role in quantifying cardiopulmonary system efficiency, and as a biomarker for cardiac vagal tone and well-being, on the basis of CRF and key concept of cardiopulmonary resonance.

Figure 1A shows the power spectral curves of RR interval series and respiration series, as well as the corresponding cardiopulmonary resonance function, G(f). G(f) represent the strength of RSA, the modulation of respiration to heart rate. G(f) is a monotonic function of frequency f with single peak around the main respiration frequency, can be considered as the spectral energy distribution function of cardiopulmonary resonance system. The cardiopulmonary resonance indices (CRI) consists of the following numerical measure:

FIGURE 1

Figure 1. Illustration of Cardiopulmonary Resonance Function (CRF) and Cardiopulmonary Resonance Indices (CRI). (A) the power spectral curves of RR interval series and respiration series, as well as the corresponding cardiopulmonary resonance function, G(f). (B) the schematic diagram of CRA and CRW. CRA is taken from the maximum point of G(f) (CRF); CRW, the bandwidth of CRF as indicated by the bi-directional arrow line.

A) Cardiopulmonary Resonance Amplitude (CRA) is defined as the maximum of Cardiopulmonary resonance function G(f):

$\begin{array}{l}CRA\equiv M\text{ax}CRF\end{array}$

Refer to Figure 1B, in consideration of the main frequency bands of heart rate variability and respiration rate, G(f) is plotted in the frequency range of 0.0033–0.5 Hz. Denote the frequency where the maxima of G(f) appears as cardiopulmonary resonance frequency f_A. CRF is around main respiration rate. In free breathing, respiration rate is around 0.20–0.30 Hz, in the range of HRV high frequency. That is the point of consistence between HRV_HF measure and RSA strength in representing the regal level. As we will see in the next section paced breathing experiments, as paced breathing frequency down to 0.1 Hz, both RSA energy and HRV energy shall move and focus around 0.1 Hz as well. In this case, the hypothesis of HRV_HF representing regal activity may not hold.

B) Cardiopulmonary Resonance bandWidth (CRW). As shown in Figure 1B, CRW is defined as the CRF bandwidth, the degree of RSA energy concentration. CRA and CRW are related. While CRW is narrow, CRA is big.

C) Cardiopulmonary Resonance Quality factor (CRQ). CRQ is defined to measure the merit of the cardiopulmonary resonance system by adopting the quality factor measure for inductor, capacitor, and resistor LCR oscillator where interaction between lung and heart resemble the energy flow between inductor and capacitor, while non-respiration factors are equivalent to resistor, damping the resonance. Mathematically, CRQ is defined as

$\begin{array}{l}CRQ=\frac{f_A}{CRW}\end{array}$

Considering the physiological functions, RSA serves to minimize the energy expenditure of the heart while keeping arterial CO₂ levels at physiological tensions. CRQ measures the energy conversion of the system. The lower the dissipation energy, the higher the quality factor and metabolic efficiency. High CRQ indicates high efficiency of cardiopulmonary metabolic system and relatively healthy physiological and psychological state.

Applications Scenarios

In this section, two application scenarios are presented to demonstrate the descriptive power of CRF and CRI, as well as the application potentials.

Paced Breathing

Experiment Design

HRV biofeedback has been used for the treatment of depression and other autonomic related problems. HRV biofeedback uses HRV measures, mainly time domain and frequency domain, as feedback cues to guide the subject performing slow paced breathing in order to reach resonance state. The objectives of paced breathing in HRV biofeedback is to gain level of parasympathetic nerves activity and improve the autonomic balance. As such, the measures of current status of the subject play most important role in biofeedback process. So far in the HRV biofeedback HRV measures are used, while HRV measures have problems in representing autonomic regulation status (Vaschillo et al., 2006).

The essential physiological phenomenon of the slow and deep paced breathing in HRV biofeedback is respiration sinus arrhythmia (RSA). The level of RSA should be the natural measure as biofeedback cues. As the quantitative measure of RSA, CRF and CRI provide the best visual cue and numerical cues for biofeedback.

During the paced breathing, Cardiopulmonary Resonance Amplitude (CRA) could help us find the optimal respiratory rate for individuals which is usually around 0.1 Hz. The process of training is the process of making CRA keep approaching 1. As we go from the resting state to paced breathing rate coming down to 0.1 Hz, with the frequency decreases, CRA gets bigger and the bandWidth CRW gets smaller. The frequency of obtaining the maximum value of CRA is the personalized resonance frequency of the subject and also the frequency of biofeedback. As an indicator of cardiopulmonary system metabolism, cardiopulmonary resonance quality factor (CRQ) indicates efficiency of cardiopulmonary metabolic system and relatively healthy physiological and psychological state.

This study was carried out in accordance with the recommendations of guidelines of ethical review of clinical research ethics committee of China-Japanese Friendship Hospital. The number is 2019-GZR-138. The protocol was approved by the clinical research ethics committee of Beijing China-Japanese Friendship Hospital. All participants signed informed consent forms. We collected data from 30 healthy adults in ages of 20–30. The baseline demographic characteristics of 30 participants are shown in Table 1. The subjects' age, height, weight, and mean systolic and diastolic blood pressure were counted and presented as mean ± standard deviation. The paired t-test showed there were no significant differences in age, systolic and diastolic blood pressure between the male and female groups.

TABLE 1

Table 1. Baseline Demographic Characteristics of 30 Participants and p-value between 15 men and 15 women.

The data is collected using one intelligent hardware, worn on the wrist (Figure 2). We collected one-lead ECG and respiratory signals of everyone from resting to biofeedback status. The whole process is recorded. During the process of paced breathing rate down to about 0.1 Hz, we use our method to find the individual resonant frequency for every trainee: Starting from the resting state of the subjects, the breathing rate was gradually reduced at 0.01 Hz intervals guided by voice and image on the computer. Each breathing rate was maintained for at least 1 min. It can be seen that Cardiopulmonary Resonance Amplitude (CRA) gradually increases as the respiratory rate decreases and reaches its maximum value around 0.1 Hz. The frequency of obtaining the maximum value of CRA is the personalized resonance frequency of the subject and also the frequency of biofeedback.

FIGURE 2

Figure 2. The wearable hardware device used to collect ECG and respiration signals.

Statistical Analysis

To demonstrate the advantage of CRI in paced breathing compared to HRV, we calculated the Cardiopulmonary Resonance Indices (CRI) and HRV in different statues. In order to represent CRF and CRA visually, we draw the CRF curves in the frequency domain with HRV in four status of paced breathing from resting status to biofeedback status. The repeated one-way ANOVA, followed by Dunnett's post hoc test was used to represent the significant difference from resting state to biofeedback state of CRI in the breathing training.

CRF and CRI in Paced Breathing

CRF measures the effect of respiration on current heart rate changes in the frequency domain. CRF and corresponding HRV, R-HRV, and NR-HRV in the frequency domain are shown in Figure 3.

FIGURE 3

Figure 3. G(f) (CRF) and corresponding HRV, R-HRV, and NR-HRV in frequency domain.

Respiratory effects in different physiological states have different effects on heart rate. These effects can be directly seen from the power spectrum calculated by spectral G-causality, which is closely related to the current breathing rate of the subjects. CRF expresses the cardiopulmonary interaction at the current time in the frequency domain.

To demonstrate the advantage of CRI in paced breathing compared to HRV, the CRF, HRV, and respiratory power spectral density of one subject of the 30 participants in the experiment from resting to biofeedback status are illustrated in Figure 4. The blue lines show respiratory power spectral density, orange lines show HRV and black lines show CRF. From the top: free breathing, paced breathing at frequency of 0.33, 0.26, and 0.12 Hz. It can be seen that paced breathing increases the strength of RSA, and that as the frequency of paced breathing coming down toward, the cardiopulmonary resonance phenomenon becomes stronger, which is very well-captured by the cardiopulmonary resonance function. We can also see that as the paced breathing frequency approaches 0.1 Hz, around which there is an optimal resonance state for the subject, where HRV HF is small. That is to say that CRI do represent level of parasympathetic nervus activity at various cases, while HRV do not.

FIGURE 4

Figure 4. A typical power spectral density curves of respiratory, HRV and corresponding cardiopulmonary resonance function CRF in the frequency domain of a subject for 4 cases of breathing at resting state. From the top: free breathing, paced breathing at frequency of 0.33, 0.26, and 0.12 Hz. It can be seen that paced breathing increases the strength of RSA, and that as the frequency of paced breathing coming down toward, the cardiopulmonary resonance phenomenon becomes stronger, which is very well-captured by the cardiopulmonary resonance function. It can also seen that as the paced breathing frequency approaches 0.1 Hz, around which there is an optimal resonance state for the subject, where HRV_HF is small. That is to say that CRI do represent level of parasympathetic nervus activity at various cases, while HRV do not.

The higher degree of cardiopulmonary coupling during paced breathing, the respiration accounts for a higher proportion of HRV. In the resting state, HRV produced by breathing is weak, then NR-HRV can reflect the influence of other physiological activities on heart rate through autonomic nervous activity.

As can be seen from the figure, with the change of respiratory rate, the distribution of CRF and HRV both shift in the corresponding frequency bands. The distribution of the frequency band of HRV is closely related to the change of respiratory frequency, so the degree of biofeedback can be observed from respiration frequency shift (Vaschillo et al., 2006). However, there is no quantitative measure between the peak value of HRV and the respiratory frequency within HRV biofeedback to describe the intensity and depth of cardiopulmonary interaction. Meanwhile, HRV cannot be used as a measure of RSA due to its low repeatability and large individual differences.

On contrary, CRF show a clear trend in paced breathing. Generally, we calculated CRI in four different states from resting to biofeedback status (from 1 to 4) in Table 2. The repeated one-way ANOVA was used to test the differences of CRI in the four states (p < 0.05). In Table 2, the p-values of CRA, CRW and f_Aare smaller than 0.05, and the p-value of CRQ is bigger than 0.05. The results showed that CRA, CRW and f_A have significant differences in the four status. CRQ is defined to measure the merit of the cardiopulmonary resonance system, and did not change significantly at different respiratory rates. In order to confirm the significance and stability of the differences between CRA and CRW in the four states further, we performed Dunnett's post hoc test shown in Table 3. As we can see, CRA increases and CRB decreases during the training. In the most of the pairwise comparisons of state 1, 2, 3, and 4, CRA and CRW show significant differences. These two indicators together represent the intensity of a person's cardiopulmonary interaction and reflect the activity and regulatory capacity of the human vagus nerve with repeatability and stability.

TABLE 2

Table 2. Cardiopulmonary Resonance Indices for the 4 cases of breathing: free breathing and 3 paced breathing at frequency of 0.33, 0.26, and 0.12Hz.

TABLE 3

Table 3. Dunnett's post-hoc test of CRA and CRW for the 4 cases of breathing: free breathing (1), and 3 paced breathing at frequency of 0.33 Hz (2), 0.26 Hz (3), and 0.12 Hz (4).

We can see that CRF and CRI could provide ideal visual interpretation and numerical measures for cardiopulmonary interactions toward resonance state in paced breathing scenario. CRF can be used to analyze the human body in different physiological states, get cardiopulmonary coupling value accurately, and analyze the regulation process of human sympathetic and parasympathetic nerves.

Sleep Stage Discrimination

Experiment Design

Cardiopulmonary Coupling index (CPC) was proposed by Thomas et al. (2005) in 2005, which measures the spectral correlation between heart rate sequence and respiratory signal. Therefore, CPC can be a candidate providing measures for RSA. CPC is defined as the product of the average cross-spectral power divided by the average power of each signal as below.

$\begin{array}{l}CPC({\text{f}}_n)={\langle {\Gamma }_n(R,E)\rangle }^2{\Lambda }_n\end{array}$

$\begin{array}{l}{\Lambda }_n=\frac{{\langle {\Gamma }_n(R,E)\rangle }^2}{\langle {\widehat{R}}_n^2\rangle \langle {Ê}_n^2\rangle }\end{array}$

in which Γ_n(R, E) denotes the cross spectrum of RR intervals and respiratory signals. CPC reflects the degree of sleep and respiratory rhythm disorder through the high frequency, low frequency and very low-frequency parts with a different energy. It overcomes the shortcomings and defects of the HRV method used in the analysis alone. At present, this method has been widely used in the field of evaluating sleep quality and judging sleep and breathing disorders (Yang et al., 2011).

In sleep stage discrimination, we used data from the MIT-BIH database (Ichimaru and Moody, 1999) which has sleep stage labels from polysomnography (PSG). CRI in different sleep stage was calculated and a comparative study was conducted by using Cardiopulmonary Coupling index (CPC) in sleep stage classification. Except for the heart rate, HRV, and respiratory rate, CRA, CRW, and f_A extracted from CRF, meanwhile LF, HF, and LF/HF extracted from CPC are, respectively, used in the classification task to test sleep stage of the subject.

In the experiment, the classifier needs to identify four different sleep stages, including awake, REM, light sleep, and deep sleep. SVM is not able to solve multi-category classification problems directly, but the combination of SVM and decision tree (called DTB-SVM) can be used to solve multi-class classification problems. Based on the structural characteristics of the sleep cycle and the physiological features used in sleep classification, we used three SVM models to classify the sleep stages. As shown in Figure 5, firstly, a classifier is used to separate the awake phase and the sleep phase, and then within the sleep phase, the REM phase and the NREM phase are separated, and finally the light sleep and deep sleep are separated by the last classifier. RBF kernel was used in the model. In order to prevent over-fitting, we selected the optimal parameters in the way of K-fold cross-validation and grid search. For the features of CPC and CRI, classifiers of same structure were used to classify the sleep stages.

FIGURE 5

Figure 5. DTB-SVM model for classifying the sleep stages. First, a classifier is used to separate the awake phase and the sleep phase, and then within the sleep phase, the REM phase and the NREM phase are separated, and finally the light sleep and deep sleep are separated by the last classifier.

Statistical Analysis

The confusion matrix was used to explain the accuracy of sleep classification results and to compare the performance of CRI and CPC in classification tasks. Each row of the matrix represents the prediction category, and the total number of each row represents the number of data predicted for that category. Each column represents the true category to which the data belongs, and the total number of data in each column represents the number of data instances in that category.

The effectiveness of the features was measured by the ratio of intra-class divergence and inter-class divergence. The ratio between intra-class divergence and inter-class divergence is defined and calculated as follows (Zhou et al., 2010):

Let (X, y)ϵ(Rⁿ × y) be a sample, where Rⁿ is an n-dimensional feature space and y = {1,2,.,s} is the label set. L_i is the number of samples in the ith class, and l is the total number of samples. Let X_ij denote the jth sample in the ith class, m_i the sample mean of the ith class, and m the sample in the ith class, m_i the sample mean of the ith class, and m the sample mean of all class. The within-class scatter matrix S_W, between-class scatter matrix S_B are defined as

$S_w={\displaystyle \sum _{i=1}^s{\displaystyle \sum _{j=1}^{l_i}(X_{ij}-m_i)(X_{ij}-m_i}}{)}^T$

$S_B={\displaystyle \sum _{i=1}^s{l}_i}\left(m_i-m\right){\left(m_i-m\right)}^T$

Large class separability means small within-class scattering and large between-class scattering. A combination of two of them can be used as a measure, |S_W |/|S_B|, where |·| denote the determinant of a matrix. The smaller the ratio, the better the effect of the feature on classification.

In order to demonstrate the good performance of CRI in sleep classification task, we compared the confusion matrix of CRI and CPC classification results. Then, in order to express the role of features further, we conducted the divergence analysis on the features of CRI and CPC. The results showed that CRI was more effective than CPC in sleep classification task, especially in the deep sleep recognition.

CRI in Different Sleep Stages Compared to CPC

To visually compare the difference between CRI and CPC, Figure 6 shows the CPC, HRV, and CRF of one subject. It contains HRV (0.14 Hz as the demarcation line between high frequency and low frequency), respiration power spectral density, CPC index and G(f), which are all discussed in the frequency domain.

FIGURE 6

Figure 6. G(f) (CRF), CPC, power spectral density of HRV and respiration. By definition, cardiopulmonary resonance function CRF, reflects the strength of RSA, with the peak near the mean of the respiration rate, while CPC is the correlation between HRV and respiration with respect to their power spectrum, having multiple peaks. As the name indicated, CRF well capture the resonance nature of cardiopulmonary system.

CPC represents the correlation of RR intervals and respiratory signal. It shows that CPC has multiple peaks in the full frequency band. In low-frequency band, for HRV analysis, people usually think sympathetic nerve and parasympathetic nerve interact together, and CPC also shows a spike which indicates a high correlation between respiration and RR intervals, such as blood pressure, etc. It is difficult to find an exact indicator representing the cardiopulmonary coupling state from CPC. Physiologically, RSA, the strength of respiration modulation of heart rate should appear as CRF, cannot be multiple peaks as CPC.

The indices of one subject in four different sleep stage of one night are shown in Table 4. It shows that our indices could express the cardiopulmonary interaction phenomenon and the degree of cardiopulmonary coupling resonance in different sleep stages. Table 5 shows the performance of CRI indicators and CPC indicators on the whole data set in the classification task. CRW, CRA, f_Aare smaller than LF_CPC and HF_CPC. The results show that the CRI features including CRW, CRA, f_A performed better than CPC features including LF_CPC and HF_CPC in the classification task.

TABLE 4

Table 4. Cardiopulmonary Resonance Indices of one subject in different sleep stages of a whole night.

TABLE 5

Table 5. The divergence analysis of the features of CPC and CRI in the sleep classification task.

The classification results of CRI and CPC were statistically analyzed, respectively, in Table 6. It shows that CRI shows superior in distinguishing deep sleep stage than CPC. The overall accuracy of the classification went up 1.28%. Particularly shown in Table 6, great progress has been made in distinguishing between deep sleep and light sleep, and the recognition rate of deep sleep has been increased by 11.35%. It shows that CRF performances better than CPC, especially in the distinction between light sleep (NREM_1 and NREM _2) and deep sleep (NREM_3 and NREM _4).

TABLE 6

Table 6. Confusion matrix of sleep stage classification using CPC and CRI.

The performance of CPC and CRF in deep sleep and light sleep is shown in Figure 7. In the deep sleep stage, CRF shows more concentrated and indicators we proposed are good indications of this phenomenon. It provides meaningful features for the distinguishing of the two. To illustrate the role of CRI in deep sleep recognition, divergence analysis of the features of CPC and CRI in the deep sleep and light sleep was performed in Table 7. It shows the performance of CRI indicators and CPC indicators for distinguishing the deep sleep stage. CRW, CRA, f_Aare much smaller than LF_CPC and HF_CPC. It shows that CRI features CRW, CRA, f_A performed much better than CPC features LF_CPC and HF_CPC. In the deep sleep, the cardiopulmonary system has the highest metabolic efficiency and the smallest dissipated energy, and the body and mind of the human body can fully rest and recover. This suggests that CRI is a good indicator for different sleep status especially the deep sleep of human body.

FIGURE 7

Figure 7. G(f) (CRF) and CPC in light sleep (top) and deep sleep (bottom).

TABLE 7

Table 7. The divergence analysis of the features of CPC and CRI in the deep sleep and light sleep.

The shortcomings of the CPC are obvious. CPC calculates the correlation between RR interval and respiratory signal, with the shape of multiple peaks. Physiologically, RSA, the strength of respiration modulation of heart rate should appear as CRF, cannot be multiple peaks as CPC. Through the study of CRI and Cardiopulmonary Coupling (CPC) in distinguishing deep sleep stage, we got the conclusion that CRI does capture physiologically meaningful characteristics of RSA, therefore, well reflect autonomic status in sleep stages. CRI represents the degree of cardiopulmonary resonance, and reflects parasympathetic nerve activity level well.

Conclusion

Respiratory sinus arrhythmia (RSA) represents a physiological phenomenon of cardiopulmonary interaction. It is known as a measure of efficiency of the circulation system, and a biomarker of cardiac vagal and well-being. In this article, we model RSA as modulation of heart rate by respiration in an interactive cardiopulmonary system with the most effective system state of resonance. Mathematically, it is described by bivariate autoregressive model of respiration series and RR intervals, and quantitatively it is assessed by Granger causality function. The whole model is referred to as Cardiopulmonary Resonance Model (CRM). This method has significant physiological significance in the frequency domain and is convenient for us to explain the experimental results. We suggest using this approach as a universal prepossessing technique which allows a researcher to concentrate on particular properties of the HRV data. Then based on the cardiopulmonary resonance concept, and Granger causality function which is referred to as cardiopulmonary resonance function (CRF) here after, a set of quantitative measures for RSA is proposed, and referred to as Cardiopulmonary Resonance Indices (CRI).

To show the effectiveness of CRM and CRI, two application scenarios, paced breathing and sleep stage discrimination, are studied. It is shown that CRF and CRI provide ideal visual interpretation and numerical measures for cardiopulmonary interactions toward resonance in paced breathing scenario as the paced breathing rate coming down to biofeedback status, and as the sleep stage moves to deep sleep. We draw the conclusion that CRI well represents the degree of cardiopulmonary resonance, and reflects parasympathetic nerve activity level. We think it's a good explanation of the physiological function of RSA and it is also good way to quantify the well-being of human body.

This study has certain limitations. As a measure of RSA under static conditions, CRI was not compared with sympathetic and parasympathetic activity indexes obtained by tilt experiment, nor was it tested under different pathological conditions. We plan to carry out relevant research in the future. In addition, we plan to explore the clinical significance of CRQ and the effects of other physiological activities on heart rate based on NR-HRV data in time and frequency domain. We will gradually accurately analyze the regulatory effects of the autonomic nervous system on various physiological organs and activities through the regulation activities of the autonomic nervous system. It is of great significance for us to understand and monitor the regulation process of the autonomic nervous system in different physiological states.

Data Availability Statement

All datasets generated for this study are included in the article/supplementary material.

Author Contributions

JC, JW, and ZH designed the study and performed the research. JC analyzed the data and wrote the paper. ZH and JW reviewed and edited the manuscript. HJ designed and assisted the experiment. All authors read and approved the manuscript.

Funding

School of Electronic, Electrical and Communication Engineering, University of Chinese Academy of Sciences provided financial support for the experimenters. This research was partially supported by National Natural Science Foundation of China, contract 61431017. CAS Institute of Healthcare Technologies provided the equipment of this research.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Thanks are due to Digital Health Research Institute of Zhong Ke (Nanjing) for provision of laboratory equipment and to Xiang'ao Meng and Zhongdi Liu for the assistance with the experiments.

References

Akaike, H. (1974). A new look at the statistical model identification. IEEE Trans. Autom. Control 19, 716–723. doi: 10.1007/978-1-4612-1694-0_16

CrossRef Full Text | Google Scholar

Appel, M. L., Berger, R. D., Saul, J. P., Smith, J. M., and Cohen, R. J. (1989). Beat to beat variability in cardiovascular variables: noise or music? J. Am. Coll. Cardiol. 14, 1139–1148. doi: 10.1016/0735-1097(89)90408-7

CrossRef Full Text | Google Scholar

Bernston, G. G., Bigger, J. T., Eckberg, D. L., Grossman, P., Kaufmann, P. G., Malik, M., et al. (1997). Heart rate variability: origins, methods, and interpretative caveats. Psychopathy 34, 623–648. doi: 10.1111/j.1469-8986.1997.tb02140.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Engel, A., and Singer, W. (2001). Temporal binding and the neural correlates of sensory awareness. Trends Cogn. Sci. 5, 16–25. doi: 10.1016/S1364-6613(00)01568-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Fonseca, D. S., Beda, A., Miranda de, S. Á., Antonio, M. F. L., and Simpson, D. M. (2013). Gain and coherence estimates between respiration and heart-rate: differences between inspiration and expiration. Auton. Neurosci. 178, 89–95. doi: 10.1016/j.autneu.2013.03.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Granger, C. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica 37, 424–438. doi: 10.2307/1912791

CrossRef Full Text | Google Scholar

Grossman, P., and Taylor, E. W. (2007). Toward understanding respiratory sinus arrhythmia: relations to cardiac vagal tone, evolution and biobehavioral functions. Biol. Psychol. 74, 263–285. doi: 10.1016/j.biopsycho.2005.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayano, J., and Yuda, E. (2019). Pitfalls of assessment of autonomic function by heart rate variability. J. Physiol. Anthropol. 38:3. doi: 10.1186/s40101-019-0193-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Ichimaru, Y., and Moody, G. B. (1999). Development of the polysomnographic database on CD-ROM. Psychiatry Clin. Neurosci. 53, 175–177. doi: 10.1046/j.1440-1819.1999.00527.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C., Amaral, L. A. N., Goldberger, A. L., Havlin, S., Rosenblum, M. G., Struzik, Z., et al. (1999). Multifractality in human heartbeat dynamics. Nature 399, 461–465. doi: 10.1038/20924

PubMed Abstract | CrossRef Full Text | Google Scholar

Jan, H. Y., Chen, M. F., Fu, T. C., Lin, W. C., Tsai, C. L., and Lin, K. P. (2019). Evaluation of coherence between ECG and PPG derived parameters on heart rate variability and respiration in healthy Volunteers With/Without Controlled Breathing. J. Med. Biol. Eng. 39, 783–795. doi: 10.1007/s40846-019-00468-9

CrossRef Full Text | Google Scholar

Katona, P. G., and Jih, F. (1975). Respiratory sinus arrhythmia: a noninvasive measure of parasympathetic cardiac control. J. Appl. Physiol. 39, 801–805. doi: 10.1152/jappl.1975.39.5.801

PubMed Abstract | CrossRef Full Text | Google Scholar

Lewis, G. F., Furman, S. A., McCool, M. F., and Porges, S. W. (2012). Statistical strategies to quantify respiratory sinus arrhythmia: are commonly used metrics equivalent? Biol. Psychol. 89:349364. doi: 10.1016/j.biopsycho.2011.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Moser, M., Fruhwirth, M., Peter, R., and Winker, R. (2006). Why life oscillates - from a topographical towards functional chronobiology. Cancer Causes Cont. 17, 591–599. doi: 10.1007/s10552-006-0015-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunez, P. L., Wingeier, B. M., and Silberstein, R. B. (2001). Spatial-temporal structures of human alpha rhythms: theory, microcurrent sources, multiscale measurements, and global binding of local networks. Hum. Brain Mapp. 13, 125–164. doi: 10.1002/hbm.1030

PubMed Abstract | CrossRef Full Text | Google Scholar

Panaite, V., Hindash, A. C., Bylsma, L. M., Small, B. J., Salomon, K., and Rottenberg, J. (2016). Respiratory sinus arrhythmia reactivity to a sad film predicts depression symptom improvement and symptomatic trajectory. Int. J. Psychophysiol. 99, 108–113. doi: 10.1016/j.ijpsycho.2015.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, C. J., and Crowell, S. E. (2016). Respiratory sinus arrhythmia as a potential measure in substance use treatment-outcome studies. Addiction 111, 615–625. doi: 10.1111/add.13232

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritz, T., Bosquet Enlow, M., Schulz, S. M., Kitts, R., Staudenmaier, J., and Wright, R. J. (2012). Respiratory sinus arrhythmia as an index of vagal activity during stress in infants: respiratory influences and their control. PLoS ONE 7:e52729. doi: 10.1371/journal.pone.0052729

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmitt, D. T., and Ivanov, P. C. (2007). Fractal scale-invariant and nonlinear properties of cardiac dynamics remain stable with advanced age: a new mechanistic picture of cardiac control in healthy elderly. Am. J. Physi. Regul. Integr. Comp. Physiol. 293, R1923–R1937. doi: 10.1152/ajpregu.00372.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. (1996). Heart rate variability. Circulation 93, 1043–1065.

Google Scholar

Thomas, R. J., Mietus, J. E., Peng, C. K., and Goldberger, K. L. (2005). An electrocardiogram-based technique to assess cardiopulmonary coupling during sleep. Sleep 28, 1151–1161. doi: 10.1093/sleep/28.9.1151

PubMed Abstract | CrossRef Full Text | Google Scholar

Topcu, C., Fruehwirth, M., Moser, M., Rosenblum, M., and Pikovsky, A. S. (2018). Disentangling respiratory sinus arrhythmia in heart rate variability records. Physiol. Meas. 39:054002. doi: 10.1088/1361-6579/aabea4

PubMed Abstract | CrossRef Full Text | Google Scholar

Tracey, K. J. (2002). The inflammatory reflex. Nature 420, 853–859. doi: 10.1038/nature01321

PubMed Abstract | CrossRef Full Text | Google Scholar

Vaschillo, E. G., Vaschillo, B., and Lehrer, P. M. (2006). Characteristics of resonance in heart rate variability stimulated by biofeedback. Appl. Psychophysiol. Biofeedback 31, 129–142. doi: 10.1007/s10484-006-9009-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Wielgus, M. D., Aldrich, J. T., Mezulis, A. H., and Crowell, S. E. (2016). Respiratory sinus arrhythmia as a predictor of self-injurious thoughts and behaviors among adolescents. Int. J Psychophysiol. 106, 127–134. doi: 10.1016/j.ijpsycho.2016.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, A. C., Yang, C. H., Hong, C. J., Tsai, J., Kuo, C. H., Peng, C. K., et al. (2011). Sleep state instabilities in major depressive disorder: detection and quantification with electrocardiogram-based cardiopulmonary coupling analysis. Psychophysiolog 48, 285–291. doi: 10.1111/j.1469-8986.2010.01060.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, L., Wang, L., and Shen, C. (2010). Feature selection with redundancy-constrained class separability. IEEE Transac. Neural Netw. 21, 853–858. doi: 10.1109/TNN.2010.2044189

PubMed Abstract | CrossRef Full Text | Google Scholar