Examining the practical importance of nonstationary cardio-respiratory coupling detection in breathing training: a methodological appraisal

Abstract

This study investigates changes in cardiorespiratory coupling during clinic breathing training and its impact on autonomic nervous functioning compared with heart rate variability (HRV). A total of 39 subjects undergoing dynamic electrocardiogram-recorded breathing training were analyzed. Subjects were divided into early- and late-training periods, and further categorized based on changes in HRV indexes. Subtypes were identified using time-frequency cardiorespiratory coupling diagrams. Significant differences were observed in the high-frequency (HF) index between training stages in the subgroup with increasing HF-HRV (p = 0.0335). Both unimodal and bimodal subtypes showed significant high-frequency coupling (HFC) in the mid-training period compared with early and late stages (both p < 0.0001), suggesting improved parasympathetic cardiac regulation or reduced sympathetic control. This study highlights the potential of nonstationary cardiorespiratory coupling analysis alongside traditional HRV in evaluating the therapeutic effect of breathing training on autonomic nervous function. Cardiorespiratory coupling analysis could provide valuable adjunctive information to HRV measures for assessing the impact of breathing training.

Introduction

Reduction of respiratory rate was known to reduce anxiety, e.g., to control physiological arousal under stressful conditions (Wells et al., 2012). Breathing training elicits changes in the autonomic nervous system (ANS) (Paprika et al., 2014), and advances control of autonomic cardiovascular balance (Besnier et al., 2017), which has been reported to reduce respiratory rate and blood pressure, having positive effects on chronic heart failure (Altena et al., 2009; Chaddha, 2015; Drozdz et al., 2016). In brief, breathing training increases parasympathetic cardiac modulation and reduces sympathetic cardiac control (Chang et al., 2015).

Heart rate variability (HRV) (Faust et al., 2022), the variation of the periods between consecutive heartbeats (RR intervals) (Vanderlei et al., 2009), is thought to reflect the adaptability of ANS to varying physiological or pathophysiological states. Evaluations of heart rate variability (HRV) measure cardiac frequency modulations, which reflect the balance between the sympathetic nervous system (SNS) and the parasympathetic nervous system (PNS), with the former associated with acceleration and the latter dependent on the deceleration of heart rate (Adjei et al., 2019). Such a set of HRV measures is accessible with non-invasive monitoring to evaluate autonomic nervous function (Rajendra Acharya et al., 2006; Faust et al., 2022), reduced in conditions such as myocardial infarction (Carney et al., 2001), congestive heart failure (Zhang et al., 2013), sudden cardiac death (Chattipakorn et al., 2007), and hypertension (Pagani & Lucini, 2001). Two main methods for quantifying heart rate variability (HRV) include time-domain indicators and frequency-domain indicators (Camm, 1996). Time domain indicators are based on direct statistical analysis of RR/NN intervals and are particularly sensitive to short or transient changes in heartbeat (Barbieri et al., 2005), while frequency domain indices derived from spectral analysis of RR intervals emphasize the role of autonomic nervous system balance in regulating heart rate dynamics. Three spectral ranges of major interest are high-frequency (HF: 0.15–0.4 Hz), low-frequency (LF: 0.04–0.15 Hz), and very-low-frequency (VLF: <0.04 Hz) bands, corresponding to short-term (5–30 min) to long-term (24 h) electrocardiogram (ECG) recordings. Comparisons and guidance in HRV indices have been well-established (Stein & Kleiger, 1999). However, these HRV measures share indirect correlations to cardiovascular function, which matters to many diseased states (e.g., congestive heart failure), drug usage (e.g., parasympathomimetic drugs), and physiological states (e.g., age) or treatment (e.g., breathing training) (Del Pozo et al., 2004).

Spontaneous resting respiratory frequency is typically around 9-to-24 breaths per minute (bpm) in healthy adults (Song & Lehrer, 2003). Apart from the HRV or respiratory indices (e.g., SNS regulatory capacity) which were reported to share partial linkages to respiratory frequency (Vaschillo, Vaschillo & Lehrer, 2006), respiratory sinus arrhythmia (RSA), which is modulated by the cardiac PNS, occurs in the HF range (0.15–0.4 Hz) (Lehrer, Vaschillo & Vaschillo, 2000; Lehrer, 2007). Numerous research efforts have been dedicated to the intensive exploration of respiratory sinus arrhythmia (RSA) (Zucker et al., 2009). Initially, Thomas et al. (2005) introduced a cardiopulmonary coupling (CPC) method to characterize sleep dynamics using a single-lead ECG. In addition to RR intervals, ECG-derived respiration (EDR) was also estimated as a surrogate respiratory signal. CPC measures the coupling between HRV (RR intervals) and amplitude modulation of ECG driven by respiration (EDR signal). By incorporating RR intervals and EDR signals, a sleep spectrogram was generated to visualize CPC dynamics during sleep (Castaldo et al., 2016). During the past years, CPC studies have been more focused on sleep-state analyses (Seo et al., 2020; Al Ashry, Ni & Thomas, 2021). Compared to the conventional HRV approaches, the CPC indexes were reported to reflect alternating sympathetic and parasympathetic modulations during sleep (Zhu et al., 2020). In recent years, more studies have highlighted the importance of RSA and cardiorespiratory function (Lin et al., 2019), interactions of breathing and heart rate in particular. In 2019, Yeh et al. (2019) verified the importance of RSA in improving cardiovascular function with its linkage to unpredictability and fractal characteristics of HRV. Later, Lin et al. (2020) proposed a multimodal coupling estimate to reveal the dynamics in the parasympathetic nerve and RSA by aging. Recently, Wang, Shi & Yeh (2023) used the synchrosqueezing transform to improve the traditional CPC algorithm, having superior time-frequency resolution, surpassing the baseline method in detecting obstructive sleep apnea (OSA). Yeh et al. (2023) explore the relationship between cross-frequency coupling and neural modulation. To the best of our knowledge, no empirical research exploring the correlations between CPC and ANS under a conscious state.

To this end, the primary objective of this study is to investigate the effects of respiratory training on the autonomic nervous system activity in patients with respiratory disorders. Specifically, we hypothesize that respiratory training will enhance cardiorespiratory coupling (CRC), which serves as a sensitive indicator of the interplay between the sympathetic and parasympathetic nervous systems. CRC reflects the temporal coordination between heart rate variability and respiratory rate, providing a window into the autonomic regulation of cardiovascular function during respiratory challenges. The rationale for selecting CRC as the primary outcome measure lies in its ability to non-invasively assess autonomic balance, which is crucial for understanding the physiological responses to respiratory training and their potential therapeutic implications. By elucidating the mechanisms underlying the effects of respiratory training on CRC, we aim to contribute to the growing body of literature on the use of non-pharmacological interventions to improve autonomic function and overall health outcomes.

Materials & MethodsSubjects and experiments

Thirty-nine patients with respiratory diseases (including chronic obstructive pulmonary disease, pulmonary heart disease, chronic bronchitis, pulmonary fibrosis, etc.) or heart diseases (including heart failure, post-PCI, myocardial infarction, arrhythmia, coronary heart disease, atrial fibrillation, cardiac insufficiency, etc.) were recruited by the Respiratory Department of the General Hospital of the People’s Liberation Army of China. We obtained written consent form and received informed consent from participants of our study. The study protocol was approved by the Ethics Committee on Biomedical Research, Chinese PLA General Hospital (S2023-023-02). The sample size of 39 patients was determined based on a power analysis conducted prior to the study, aiming to achieve a statistical power of at least 80% at a significance level of 0.05, given the expected effect size derived from previous research in the field. This approach ensured that the study had sufficient sensitivity to detect meaningful differences in cardio-respiratory coupling detection during breathing training. All patients enrolled in this study exhibited clinical symptoms indicative of autonomic imbalance, as manifested by irregular breathing patterns, altered heart rate variability, and/or other signs of compromised cardio-respiratory coordination. These symptoms were assessed by experienced respiratory physicians based on medical history and physical examination. Respiratory training was prescribed as a therapeutic intervention aimed at modulating these autonomic imbalances and improving cardio-respiratory coupling. All participants provided informed consent before their involvement, and their confidentiality and privacy were strictly maintained throughout the study. The clinical characteristics of these patients are summarized in Table 1, including age, height, weight, body mass index (BMI), heart rate (HR), and respiratory rate (RR). Data collection was facilitated by a wearable device, SensEcho, which the patients wore comfortably during the study. Each patient was equipped with a bra-like garment, an electrode patch for ECG recording, a blood oxygen saturation tester, and a dynamic ECG recorder. The experimental protocol consisted of two main phases: baseline recording and breathing training. During the baseline recording phase, patients were instructed to breathe autonomously for 1 min in a relaxed state. The patients then underwent breathing training sessions, guided by both auditory and visual cues interactively. Each training session lasted for a predetermined duration of 10 min. Throughout the training, patients were instructed to maintain an inhalation-to-exhalation ratio of 1:2, aiming to optimize respiratory efficiency and promote autonomic modulation. All data were sampled at 200 Hz.

Table 1: Clinical information.

Subjects	Age (yrs)	Height (cm)	Weight (kg)	BMI (kg/m2)	HR (beats/min)	RR (breaths/min)
39	58 ± 14	164 ± 7	65 ± 14	24 ± 5	78 ± 17	18 ± 4

DOI: 10.7717/peerj.18551/table-1

Signal preprocessing

To begin with, 600-s recordings (ECG and respiratory signal) of all subjects before, during, and after respiratory training were included in the following analyses, of which the first and the second half of recordings (300 s of each) were annotated as the early-training and the late-training periods. Firstly, ECG and respiratory signals were both filtered by a 50-Hz notch filter to remove power frequency noise. The consecutive differences of R peaks of an ECG series derived RR intervals. Then, RR intervals were interpolated, and along with the respiratory signal, both were resampled into 4 Hz. Next, the HRV/CRC analyses were performed on the preprocessed data.

HRV analyses

In the HRV analysis, we conducted both the time-domain and the frequency-domain analyses of the RR intervals, with the calculated indexes of the periods before and after training carefully compared. According to the changes of the HRV frequency-domain indicators by training, all subjects were further categorized into subgroups as per indicators including high-frequency power (HF, parasympathetic-related HRV), low-frequency power (LF, sympathetic-related HRV), and LF/HF ratio. The results of the division are summarized in Table 2.

Table 2: Grouping subjects according to the HRV indicators.

Indicators	Ascending	Descending
HF	17	22
LF	24	15
LF/HF	24	15

DOI: 10.7717/peerj.18551/table-2

Nonstationary CRC method

For the nonstationary CRC method, our processing steps are as follows:

1. Uniformly interpolate respiratory and RR signals to 8 Hz.
2. Calculation of the scales for the continuous wavelet is as follows: [Formula omitted. See PDF.] where nv is the number of scales per octave, start_scale is the initial value of the scale, dt is the sampling interval, scale_factor is the scale factor, N is the data length, ncot is the number of octaves, and scales is the scale of the continuous wavelet transforms.
3. Wavelet transform was performed on the respiratory and RR signals to obtain C_R (wavelet representation of the respiratory signal) and C_E (wavelet representation of the RR signal).
4. Calculate the wavelet cross-spectrum Γ_ER of the two signals. The wavelet cross-spectrum is a measure of the distribution of power of the respiratory and RR signals. The cross-wavelet power of these two signals was represented by [Formula omitted. See PDF.] where the superscript * denotes the complex conjugate, and S is a smoothing operator in the time scale, and the smoothing factor can be set based on the data. During the smoothing process, if the size of the average filter is not specified, it defaults to the number of scale factors per octave. This means that the smoothing operation will use the number of scale factors in each octave as the window size for averaging, aiming to achieve a smoothing effect in the scale direction.
5. Calculate the wavelet coherence Λ_ER of the respiratory and RR signals. Wavelet coherence is a measure of the correlation between the two signals in the time-frequency space. The Wavelet coherence of the respiratory and RR signals is defined by [Formula omitted. See PDF.] where the “——” operator represents taking the absolute value or magnitude of a real number.
6. Calculate the cardiorespiratory coupling (CRC) of the two signals. CRC measures the degree of interaction, i.e., [Formula omitted. See PDF.]

Similar to the HRV analyses, all subjects were categorized according to the CRC dynamics in time-frequency visualization. In performing the CRC method, as per the CRC time-frequency visualization, prominent coupling was revealed during the training process, either presenting a single peak or dual peaks over the temporal domain. Following this finding, the early-training and late-training periods required adjustment; precisely, the temporal edges of the aforementioned training periods were tuned according to the peaks of the CRC diagram.

For the condition with a single peak, the interval of 25 s before and after the peak was defined as the mid-training period. The section before the mid-training period was the early-training period, while the section after the mid-training period was the late-training period. For the condition with dual peaks, we considered the segment between the two peaks as the mid-training period, and the periods before and after the mid-training period correspond to the early-training and the late-training period, respectively. Table 3 summarizes the subject numbers as per the number of CRC peaks.

Table 3: Grouping subjects according to the CRC indicators.

Indicator	Single	Dual
Peak #	17	22

DOI: 10.7717/peerj.18551/table-3

Statistical analysis

In this study, comprehensive statistical analyses were employed to systematically summarize the demographic and clinical characteristics of the study participants, as well as to evaluate the respiratory training efficacy based on various cardiorespiratory coupling indicators. Descriptive statistics, including mean and standard deviation, were utilized to summarize both patient characteristics and outcome measures. For comparative analyses, nonparametric Wilcoxon rank-sum tests were used to compare those derived features among groups with different training periods, generating p-values between groups. All statistical tests were two-sided, and statistical significance was set at p < 0.05.

ResultsHRV analysis

Firstly, we performed the HRV analyses on all subjects in the early and late stages of respiratory training. The time-domain HRV indicators are based on mean NN intervals, including SDNN (mean of the standard deviations of all NN intervals), SDSD (standard deviation of differences between adjacent NN intervals), pNN50 (number of pairs of adjacent NN intervals differing by more than 50 ms divided by the total number of all NN intervals), and pNN20 (number of pairs of adjacent NN intervals differing by more than 20 ms divided by the total number of all NN intervals), while the frequency-domain HRV indexes include the total power, VLF (power in VLF range, ≤0.04 Hz), LF (power in LF range, 0.04–0.15 Hz), HF (power in HF range, 0.15–0.4 Hz), LFn (LF power in normalized units LF/(total power −VLF) ×100), HFn (HF power in normalized units HF/(total power −VLF) ×100), and LF/HF. Table 4 summarizes comparisons of HRV indices with cardio interbeat intervals between the early and late stages of respiratory training. Unfortunately, no indicator presents a significant difference between the early stage and the late stage.

Table 4: Comparison of HRV indices between the early-training and the late-training periods.

	Early	Late	P-value
Time domain
MeanNN (ms)	832.1 ± 19.9	818.2 ± 18.9	0.6134
SDNN (ms)	42.4 ± 2.8	36.9 ± 2.7	0.1591
SDSD (ms)	12.5 ± 1.6	10.3 ± 1.4	0.3073
RMSSD(ms)	12.5 ± 10.2	10.3 ± 8.8	0.3069
pNN50 (%)	1.5 ± 0.6	1.1 ± 0.6	0.5912
pNN20 (%)	7.0 ± 1.7	5.3 ± 1.5	0.4279
Freq. domain
Total Power (lg/ms2)	3.2 ± 0.4	3.2 ± 0.4	0.7590
VLF (lg/ms2)	2.8 ± 0.1	2.6 ± 0.2	0.0947
LF (lg/ms2)	2.8 ± 0.1	2.8 ± 0.2	0.6604
HF (lg/ms2)	2.7 ± 0.1	2.8 ± 0.1	0.5147
LFn (n.u.)	52.9 ± 4.4	51.6 ± 4.4	0.8348
HFn (n.u.)	47.1 ± 4.4	48.4 ± 4.4	0.8348
LF/HF ratio	2.8 ± 0.6	2.9 ± 0.7	0.9193

DOI: 10.7717/peerj.18551/table-4

HRV-subgroup analyses

Given the above-analyzed HRV results, HRV-subgroups were categorized as per frequency-domain indexes of HRV, including changes of HRV-HF, HRV-LF, and HRV-LF/HF (i.e., ascend or descend) by respiratory training. The HRV-LF and HRV-HF correspond to parasympathetic-related and sympathetic-related HRV, respectively. Further comparisons of HRV indices between the early and the late stages of training were performed on the individual HRV-subgroup, and the comparative results are shown in Fig. 1. As shown in Fig. 1, only HF in the HF ascend subgroup showed a significant difference by training, reflecting an increase in parasympathetic function (p = 0.0355). Of note, Table 5 summarizes the significant results for all comparisons of HRV frequency domain features between the early- and late-training periods.

Figure 1: Comparison of low-frequency (HRV-LF) and high-frequency (HRV-HF) heart rate variability (HRV) indices between the early and late stages of training. Participants were grouped based on their HRV changes during training into six categories: HRV-HF ascending, HRV-HF descending, HRV-LF ascending, HRV-LF descending, HRV-LF/HF ratio ascending, and HRV-LF/HF ratio descending. DOI: 10.7717/peerj.18551/fig-1

Table 5: Comparison of HRV indices between the early-training and the late-training periods.

		Early	Late	P-value
HF ascend	LF	615.5 ± 185.8	839.1 ± 268.9	0.5518
HF	326.4 ± 76.2	795.0 ± 169.0	0.0335
HF descend	LF	543.2 ± 126.6	516.9 ± 113.6	0.8779
HF	614.3 ± 135.1	408.4 ± 94.0	0.2176
LF ascend	LF	442.6 ± 88.4	681.2 ± 121.7	0.1389
HF	466.5 ± 107.4	658.4 ± 142.0	0.3106
LF descend	LF	622.8 ± 153.7	439.6 ± 121.9	0.4473
HF	495.9 ± 136.5	386.2 ± 101.6	0.5985
LF/HF ascend	LF	579.3 ± 133.0	685.0 ± 120.3	0.5776
HF	470.8 ± 104.2	395.4 ± 106.6	0.6327
LF/HF descend	LF	440.5 ± 102.9	434.6 ± 123.3	0.9759
HF	490.3 ± 139.8	736.9 ± 145.2	0.3210

DOI: 10.7717/peerj.18551/table-5

Notes:

Bold numbers indicate p-value < 0.05.

CRC-subgroup analyses

In the analysis of time-frequency CRC dynamics, we observed two distinct subtypes among our participants. As illustrated in Fig. 2, the first subtype, referred to as the ’single-peak’ subtype, is characterized by an initial emergence of CRC coupling at the beginning of training, which gradually weakens as training progresses. The left panel of Fig. 2 depicts a typical example of this subtype, where the coupling strength peaks early in the training period and subsequently declines. The middle panel of Fig. 2 presents a variant of this subtype, where the maximum CRC coupling occurs at the midpoint of training, followed by a weakening trend. In contrast, the second subtype, designated as the ‘dual-peaks’ subtype, exhibits an initial increase in CRC coupling, followed by a local maximum, which then weakens and bounces back later in the training process. Based on these distinct temporal patterns in CRC dynamics, we categorized participants into two distinct subtypes (i.e., single-peak subtype and dual-peaks subtype) based on the visual inspection of the CRC spectrogram of each subject. For the condition with a single peak, the mid-training period was defined as the interval of 25 s before and after the peak. The section preceding the mid-training period was referred to as the early-training period, while the section following the mid-training period was known as the late-training period. On the other hand, for the condition with dual peaks, the mid-training period was redefined as the segment between the two peaks. The periods before and after the mid-training period corresponded to the early-training and late-training periods, respectively. This redefinition was done based on the timing of the occurrence of the peaks.

Figure 2: Demonstrations of the two subtypes in the time-frequency CRC dynamics, including the single-peak (left and mid panels) and the dual-peak (right panel) conditions. DOI: 10.7717/peerj.18551/fig-2

Figure 3: Comparison of cardiorespiratory coupling (CRC) indices, specifically low-frequency coupling (LFC) and high-frequency coupling (HFC), across different training stages. Data are presented separately for individuals exhibiting single-peak and dual-peak CRC patterns. DOI: 10.7717/peerj.18551/fig-3

According to the above criteria, either subtype, the nonstationary CRC analysis was applied. Next, both the low-frequency coupling (LFC, 0.1–0.4 Hz) and high-frequency coupling (HFC, 0.1–0.4 Hz) were extracted over the early, mid, and late stages of respiratory training, as compared and summarized in Fig. 3 and Table 6. The coupling strength in the HFC of the unimodal-subtype group during the mid-training period is significantly higher than that of the early and late stages (left panel in Fig. 3, p < 0.0001). Likewise, as for the bimodal-subtype group, the coupling strength in the HFC during the mid-training period is also significantly higher than that of the other two stages (right panel in Fig. 3, p < 0.0001).

Next, in Fig. 4, to compare the CRC indices (LFC and HFC) further between early and late training stages, i.e., to test if the effectiveness of training lasting, subjects were further grouped according to their HRV differences by training, briefly, either ascending (top panels) or descending (bottom panels), either with HRV-HF (left panels), HRV-LF (middle panels), or HRV-LF/HF (right panels). Unfortunately, no statistical significance was shown.

Single-peak CRC-subgroup reaction time analysis

Following the above analysis, we inferred a period of training may be required to achieve stronger CRC or sufficient training effect. Therefore, to examine the transformation of the ANS during the respiratory training process, we further investigated the reaction time (RT) to the maximum HFC (parasympathetic nervous activity) and LFC (sympathetic nervous activity), aiming to quantify changes in parasympathetic and sympathetic dominance during training. We performed the “peak-to-valley” processing on all subjects, i.e., to find the time point when HFC reaches the peak value and the time point when LFC on the trough, which we consider to be the RT for respiratory training to achieve the optimal effect. Following the aforementioned, the varying dominances of the parasympathetic and sympathetic nerves during training were set. As shown in Fig. 5, The reaction time of LFC and HFC is approximately around 300 s.

Table 6: Comparisons of CRC indices include LFC and HFC among training stages for single-peak (unimodal) and dual-peaks (bimodal) subtypes.

		Early	Mid	Late	P-value
Single peak	LFC	5.5 ± 0.5	6.2 ± 0.8	6.3 ± 0.7	0.7094
HFC	5.1 ± 0.4	12.4 ± 0.6	5.3 ± 0.5	<0.0001
Dual peaks	LFC	6.2 ± 0.6	4.8 ± 0.6	5.8 ± 0.6	0.4046
HFC	5.3 ± 0.5	14.3 ± 0.4	5.2 ± 0.5	<0.0001

DOI: 10.7717/peerj.18551/table-6

Notes:

Bold numbers indicate p-value < 0.05.

Figure 4: Comparisons of CRC indices include LFC and HFC between the early and late stages of training. Subjects were grouped as per their HRV differences by training, including HRV-HF ascend, HRV-HF descend, HRV-LF ascend, HRV-LF descend, HRV-LF/HF ascend, and HRV-LF/HF descend subtype. DOI: 10.7717/peerj.18551/fig-4

Figure 5: Reaction time (RT) analysis for participants exhibiting single-peak cardiorespiratory coupling patterns. DOI: 10.7717/peerj.18551/fig-5

Discussion

In the present study, we explored the effectiveness of nonstationary CRC and HRV in assessing respiratory training-related changes in the autonomic nervous system. Subtype-grouping techniques include unimodal/bimodal subtypes and HRV ascend/descend were applied. The reaction time of parasympathetic and sympathetic functions during the breathing training was verified.

For the HRV indices, as per Table 4, neither the time-domain indicators nor the frequency-domain indicators showed significant differences. Then, all data were further grouped as per the changes of LF (more to the sympathetic activity), HF (more to the parasympathetic activity), or LF/HF in the HRV, as per Fig. 1 and Table 5, the HRV-HF of the HF-ascend group is significantly higher in the late stage of breathing training than that of the early stage (p = 0.0355), indicating that the parasympathetic nerve dominates in the late stage of training (Chang et al., 2015).

Subsequently, the nonstationary CRC analysis was performed on all data. Followed by all subjects were further divided into the unimodal and bimodal subtypes according to their time-frequency CRC dynamics, as per Fig. 3 and Table 6. Both the unimodal and bimodal subtypes revealed significant differences during the mid-training period compared with the early and late stages. In the unimodal group, the CRC-HFC during the mid-term breathing training was significantly higher than that in the early-training and late-training periods (p <0.001), validating its role of being in an effective state of breathing training and about to enter the exhaustion stage. Moreover, from the early stage to the middle stage of training, CRC-HFC increased significantly, indicating a gradual rising training effect, along with the parasympathetic nerve’s gradual dominance, and vice versa from the middle to the late stages of training. Meanwhile, the response time of the HFC and LFC indicates that about 3 min may be required to achieve optimal respiratory training effect on the ANS function.

While the findings of this study contribute valuable insights into the effects of respiratory training on cardiorespiratory coupling and autonomic nervous system activity, several limitations should be acknowledged. Firstly, since our objective was to evaluate the practical significance of detecting nonstationary cardio-respiratory coupling during breathing training specifically in a patient population with known cardio-respiratory impairments, a healthy control group was not considered in this study. Although this approach allowed us to assess the potential clinical utility of the proposed method in a context where it is most likely to be applied (i.e., in clinical treatment practice), it limits our understanding of the effects of breathing training on normal human physiological mechanisms. Therefore, future studies should consider including a control group to strengthen the causal inference of the findings. Additionally, the study was conducted with a relatively small sample size, which may limit the generalizability of the results to broader populations. Larger, multi-center trials would be necessary to confirm the robustness and reproducibility of our findings.

In terms of novelty, this study represents a unique contribution to the field by employing cardiorespiratory coupling as a novel biomarker to assess the autonomic effects of respiratory training. Previous research has primarily focused on more traditional measures of autonomic function, such as heart rate variability or blood pressure responses. Our analysis offers a more nuanced understanding of the complex interplay between respiratory patterns and autonomic regulation, which has important implications for the development of tailored respiratory training programs for various populations.

Conclusions

In this study, we assessed the impact of breathing training on autonomic nervous and cardiorespiratory dynamics. While HRV analysis showed no significant changes pre- and post-training, except for elevated HF-HRV in the late training period within a subgroup, subjects were categorized based on temporal peaks in HFC from time-frequency CRC into early, mid-, and late-training periods. Coupling intensity in the high-frequency band during the mid-training period significantly surpassed other periods, indicating stronger changes than HF-HRV measures. Our findings underscore the value of nonstationary CRC analysis alongside HRV measures in evaluating breathing training effects, offering more comprehensive insights. Integration of CRC analysis into clinical practice could enhance the assessment and monitoring of patients undergoing breathing training interventions.

Additional Information and Declarations

Competing Interests

The authors declare there are no competing interests.

Author Contributions

Jinfeng Li analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the article, and approved the final draft.

Yong Fan conceived and designed the experiments, performed the experiments, analyzed the data, authored or reviewed drafts of the article, and approved the final draft.

Wenbin Shi conceived and designed the experiments, performed the experiments, analyzed the data, authored or reviewed drafts of the article, and approved the final draft.

Mengwei Li analyzed the data, authored or reviewed drafts of the article, and approved the final draft.

Lixuan Li analyzed the data, authored or reviewed drafts of the article, and approved the final draft.

Wei Yan analyzed the data, authored or reviewed drafts of the article, and approved the final draft.

Muyang Yan analyzed the data, authored or reviewed drafts of the article, and approved the final draft.

Zhengbo Zhang analyzed the data, authored or reviewed drafts of the article, and approved the final draft.

Chien-Hung Yeh conceived and designed the experiments, performed the experiments, analyzed the data, authored or reviewed drafts of the article, and approved the final draft.

Human Ethics

The following information was supplied relating to ethical approvals (i.e., approving body and any reference numbers):

The Ethics Committee on Biomedical Research, Chinese PLA General Hospital

Data Availability

The following information was supplied regarding data availability:

The raw ECG signals for all the use cases are available in the Supplemental File.

Funding

This work was supported by tthe National Natural Science Foundation of China (grant no. 62171028, no. 62471043, no. 62171471), the Beijing Natural Science Foundation (grant no. L232139), the Open Project of Key Laboratory of Medical Electronics and Digital Health of Zhejiang Province (no. MEDC202303), and the National High-level Fellow Research Fund Program (grant no. 3050012222022). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

Adjei T, von Rosenberg W, Nakamura T, Chanwimalueang T, Mandic DP. 2019. The classA framework: HRV based assessment of SNS and PNS dynamics without LF-HF controversies. Frontiers in Physiology 10:505

Al Ashry HS, Ni Y, Thomas RJ. 2021. Cardiopulmonary sleep spectrograms open a novel window into sleep biology—implications for health and disease. Frontiers in Neuroscience 15:755464

Altena MR, Kleefstra N, Logtenberg SJ, Groenier KH, Houweling ST, Bilo HJ. 2009. Effect of device-guided breathing exercises on blood pressure in patients with hypertension: a randomized controlled trial. Blood Pressure 18:273-279

Barbieri R, Matten EC, Alabi AA, Brown EN. 2005. A point-process model of human heartbeat intervals: new definitions of heart rate and heart rate variability. American Journal of Physiology-Heart and Circulatory Physiology 288:H424-H435

Besnier F, Labrunee M, Pathak A, Traon APavy-Le, Gales C, Senard J-M, Guiraud T. 2017. Exercise training-induced modification in autonomic nervous system: an update for cardiac patients. Annals of Physical and Rehabilitation Medicine 60:27-35

Camm J. 1996. Task force of the european society of cardiology and the north american society of pacing and electrophysiology. Heart rate variability: standarts of measurement, physiological interpretation and clinical use. Circulation 93:1043-1065

Carney RM, Blumenthal JA, Stein PK, Watkins L, Catellier D, Berkman LF, Czajkowski SM, O’Connor C, Stone PH, Freedland KE. 2001. Depression, heart rate variability, and acute myocardial infarction. Circulation 104:2024-2028

Castaldo R, Melillo P, Izzo R, De Luca N, Pecchia L. 2016. Fall prediction in hypertensive patients via short-term HRV analysis. IEEE Journal of Biomedical and Health Informatics 21(2):399-406

Chaddha A. 2015. Slow breathing and cardiovascular disease. International Journal of Yoga 8:142-143

Chang Q, Liu R, Li C, Shen Z. 2015. Effects of slow breathing rate on blood pressure and heart rate variabilities in essential hypertension. International Journal of Cardiology 185:52-54

Chattipakorn N, Incharoen T, Kanlop N, Chattipakorn S. 2007. Heart rate variability in myocardial infarction and heart failure. International Journal of Cardiology 120:289-296

Del Pozo JM, Gevirtz RN, Scher B, Guarneri E. 2004. Biofeedback treatment increases heart rate variability in patients with known coronary artery disease. American Heart Journal 147:545

Drozdz T, Bilo G, Debicka-Dabrowska D, Klocek M, Malfatto G, Kielbasa G, Styczkiewicz K, Bednarek A, Czarnecka D, Parati G. 2016. Blood pressure changes in patients with chronic heart failure undergoing slow breathing training. Blood Pressure 25:4-10

Faust O, Hong W, Loh HW, Xu S, Tan R-S, Chakraborty S, Barua PD, Molinari F, Acharya UR. 2022. Heart rate variability for medical decision support systems: a review. Computers in Biology and Medicine 145:105407

Lehrer PM. 2007. Biofeedback training to increase heart rate variability. In: Lehrer PM, Woolfolk RL, eds. Principles and practice of stress management (4th edition). New York: The Guilford Press. 264-302

Lehrer PM, Vaschillo E, Vaschillo B. 2000. Resonant frequency biofeedback training to increase cardiac variability: rationale and manual for training. Applied Psychophysiology and Biofeedback 25:177-191

Lin C, Lin P-F, Wang C-H, Juan C-H, Tran T-T, Pham V-T, Nien C-T, Lin Y-J, Wang C-Y, Yeh C-H, Lo M-T+1 more. 2020. Probing age-related changes in cardio-respiratory dynamics by multimodal coupling assessment. Chaos: An Interdisciplinary Journal of Nonlinear Science 30:033118

Lin C, Yeh H-M, Lo M-T, Yeh C-H, Wang C-Y, Shi W, Bmaf Serafico, Wang C-H, Juan C-H, Young H-WVincent, Lin Y-J+1 more. 2019. Robust fetal heart beat detection via R-peak intervals distribution. IEEE Transactions on Biomedical Engineering 66:3310-3319

Pagani M, Lucini D. 2001. Autonomic dysregulation in essential hypertension: insight from heart rate and arterial pressure variability. Autonomic Neuroscience 90:76-82

Paprika D, Gingl Z, Rudas L, Zöllei E. 2014. Hemodynamic effects of slow breathing: does the pattern matter beyond the rate? Acta Physiologica Hungarica 101:273-281

Rajendra Acharya U, Paul Joseph K, Kannathal N, Lim CM, Suri JS. 2006. Heart rate variability: a review. Medical and Biological Engineering and Computing 44:1031-1051

Seo MY, Hwang SJ, Nam KJ, Lee SH. 2020. Significance of sleep stability using cardiopulmonary coupling in sleep disordered breathing. The Laryngoscope 130:2069-2075

Song H-S, Lehrer PM. 2003. The effects of specific respiratory rates on heart rate and heart rate variability. Applied Psychophysiology and Biofeedback 28:13-23

Stein P, Kleiger RMD. 1999. Insights from the study of heart rate variability. Annual Review of Medicine 50:249-261

Thomas RJ, Mietus JE, Peng CK, Goldberger AL. 2005. An electrocardiogram-based technique to assess cardiopulmonary coupling during sleep. Sleep 28(9):1151-1161

Vanderlei LCM, Pastre CM, Hoshi RA, Carvalho TD de, Godoy MF de. 2009. Basic notions of heart rate variability and its clinical applicability. Brazilian Journal of Cardiovascular Surgery 24:205-217

Vaschillo EG, Vaschillo B, Lehrer PM. 2006. Characteristics of resonance in heart rate variability stimulated by biofeedback. Applied Psychophysiology and Biofeedback 31:129-142

Wang Y, Shi W, Yeh C-H. 2023. A novel measure of cardiopulmonary coupling during sleep based on the synchrosqueezing transform algorithm. IEEE Journal of Biomedical and Health Informatics 27(4):1790-1800

Wells R, Outhred T, Heathers JA, Quintana DS, Kemp AH. 2012. Matter over mind: a randomised-controlled trial of single-session biofeedback training on performance anxiety and heart rate variability in musicians.

Yeh C-H, Juan C-H, Yeh H-M, Wang C-Y, Young H-WV, Lin J-L, Lin C, Lin L-Y, Lo M-T. 2019. The critical role of respiratory sinus arrhythmia on temporal cardiac dynamics. Journal of Applied Physiology 127:1733-1741

Yeh C-H, Zhang C, Shi W, Lo M-T, Tinkhauser G, Oswal A. 2023. Cross-frequency coupling and intelligent neuromodulation. Cyborg and Bionic Systems 4:0034

Zhang Y, De Peuter OR, Kamphuisen PW, Karemaker JM. 2013. Search for HRV-parameters that detect a sympathetic shift in heart failure patients on β-blocker treatment. Frontiers in Physiology 4:81

Zhu T, Zhou J, Zhou J, Feng L, Yang J, Wang G. 2020. High-frequency cardiopulmonary coupling during sleep correlates with memory in depressed patients: A pilot study. Journal of Affective Disorders 270:118-123

Zucker TL, Samuelson KW, Muench F, Greenberg MA, Gevirtz RN. 2009. The effects of respiratory sinus arrhythmia biofeedback on heart rate variability and posttraumatic stress disorder symptoms: a pilot study. Applied Psychophysiology and Biofeedback 34:135-143