Acute Cardiorespiratory Coupling Impairment in Worsening Sleep Apnea-Related Intermittent Hypoxemia

Objective: Hypoxic load is one of the main characteristics of obstructive sleep apnea (OSA) contributing to sympathetic overdrive and weakened cardiorespiratory coupling (CRC). Whether this association changes with increasing hypoxic load has remained obscure. Therefore, we aimed to study our hypothesis that increasing hypoxic load acutely decreases the CRC. Methods: We retrospectively analyzed the electrocardiography and nasal pressure signals in 5-min segment pairs (n = 36 926) recorded during clinical polysomnographies of 603 patients with suspected OSA. The segment pairs were pooled into five groups based on the hypoxic load severity described with the the total integrated area under the blood oxygen saturation curve during desaturations. In these severity groups, we determined the frequency-domain heart rate variability (HRV) parameters, the HRV and respiratory high-frequency (HF, 0.15–0.4 Hz) peaks, and the difference between those peaks. We also computed the spectral HF coherence between HRV and respiration in the HF band. Results: The ratio of low-frequency (LF, 0.04–0.15 Hz) to HF power increased from 1.047 to 1.805 (p < 0.001); the difference between the HRV and respiratory HF peaks increased from 0.001 Hz to 0.039 Hz (p < 0.001); and the spectral coherence between HRV and respiration in the HF band decreased from 0.813 to 0.689 (p < 0.001) as the hypoxic load increased. Conclusion and Significance: The vagal modulation decreases and CRC weakens significantly with increasing hypoxic load. Thus, the hypoxic load could be utilized more thoroughly in contemporary OSA diagnostics to better assess the severity of OSA-related cardiac stress.

In addition to sympathovagal modulation, HR is modulated mainly by mechanical coupling and baroreflex due to respiration: HR increases during inspiration and decreases during expiration [9], [10].This phenomenon is called respiratory sinus arrhythmia (RSA) and is the primary expression of cardiorespiratory coupling (CRC) [10].The suggested physiological function is to attain a matching pulmonary blood flow to intermittent airflow [11].In the frequency-domain HRV analysis, the RSA is normally represented as a peak in the high-frequency (HF, 0.15-0.40Hz) band of the power spectral density (PSD) and it is attenuated by hypoxia [9], [10].In addition to vagal dominance, the CRC within the HF band is a biomarker in healthy adults indicating stable sleep and breathing [12], [13], [14].OSA patients, however, exhibit weakened HF coupling [12].
Although RSA and CRC have been studied on OSA patients [12], the effect of the hypoxic load severity on their HRV modulation is currently unknown.In our previous study [15], we observed diminished vagal modulation related to a more severe hypoxic load and almost indistinguishable peaks on the HF band.These results are in line with previous studies [10], [12].Due to the diminished HF band peaks with a more severe hypoxic load, the HF analysis of the HRV alone might not be the optimal method to study the RSA and the effect of hypoxic load severity on the CRC [16].
Thus, our aim was to study the association between HRV and respiration by cross-spectral coherence.The cross-spectral coherence measures the degree of correlation between the spectral components of two signals of interest, e.g., HR and respiration from the nasal pressure signal [17].By utilizing a large clinical population of suspected OSA patients, we aimed to combine spectral coherence analysis with the RSA analysis to obtain a more comprehensive understanding of CRC in OSA patients.We hypothesized that in addition to diminished vagal modulation [15], the increasing hypoxic load acutely weakens the connection between the cardiac and respiratory activities in terms of lower HF coupling.

A. Polysomnographic Data
The initial clinical population comprised 901 consecutive patients with suspected OSA undergoing type I polysomnography (PSG) in the Sleep Disorders Center of Princess Alexandra Hospital (Brisbane, Australia) between 2015 and 2017.The PSGs were recorded using the Compumedics Grael acquisition system (Compumedics, Abbotsford, Australia).The Institutional Human Research Ethics Committee of Princess Alexandra Hospital has approved the retrospective data collection and re-use (HREC/16/QPAH/021 and LNR/2019/QMS/54313). Experienced sleep specialists manually scored the PSG recordings in compliance with the American Academy of Sleep Medicine (AASM) 2012 scoring criteria [18] using the 3% desaturation rule for hypopneas.After patient exclusion, as described in Fig. 1, we included a total of 603 patients in this study (Table I).

TABLE I CHARACTERISTICS OF THE STUDY POPULATION
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 2. Illustration of used RR interval, nasal pressure, and oxygen saturation signals within one 5-min segment from the reference group, DesSev 0 (DesSev = 0%), and the group with the most severe hypoxic load, DesSev Q4 (DesSev > 2.467%).In DesSevQ4, apneas lead to cyclical variation and losing the high-frequency changes in RR intervals.They also cause desaturations that contribute to hypoxic load.The severity of the hypoxic load is measured with the desaturation severity (DesSev) parameter [21] by dividing the sum of the individual desaturation areas (DesArea) by the duration of the segment (5 min).
frequency of 256 Hz and using the modified lead II [18].The R-peaks were detected from the ECGs using the Kubios HRV Premium 3.4.1 software (Kubios Oy, Kuopio, Finland) with the default settings [19].Before the R peak detection, the software interpolated the R wave to 2000 Hz to improve the detection of the R peak location.The Kubios detects R peaks with an algorithm based on the Pan-Tompkins method [20] using the adaptive amplitude threshold and expected time between adjacent R peaks as the decision rules.In addition, the software automatically corrects the artefacts due to ectopic peaks and missed peak detections by replacing them using cubic spline interpolation.Finally, we divided the resulting R peak time series into non-overlapping 5-min segments during sleep (n = 51 157).
After excluding 14 231 5-min segments (exclusion criteria described in Fig. 1), we divided the remaining 36 926 segments into groups based on the hypoxic load during the segment, similarly as in our previous study [15].We described the severity of hypoxic load with the desaturation severity (DesSev) parameter [21].Due to the lack of strict international scoring rules for desaturations, the desaturation scoring practices vary from sleep center to center.In the dataset applied in the present study, desaturations were scored from the baseline to the baseline.Thus, in contrast to Kulkas et al. [21], the DesSev we used here also includes the recovery area from the nadir back to the baseline.The DesSev describes the severity of the hypoxic load by considering the depth and duration of the individual desaturation events as the integrated area under the blood oxygen saturation curve (Fig. 2).The DesSev is calculated by normalizing the sum of individual desaturation areas in the segment by the duration of the segment (5 min).Each desaturation event was considered to belong to that segment during which it started.The reference group, DesSev 0 , consisted of only segments having zero desaturations (DesSev = 0%, n = 22 139), and the remaining segments were pooled into quartiles DesSev Q1 to DesSev Q4 (0% < DesSev Q1 ≤ 0.301% < DesSev Q2 ≤ 0.898% < DesSev Q3 ≤ 2.467% < DesSev Q4 , Table I).
In this study, HRV was measured using short-term (5 min) frequency-domain analyses from the RR interval series to assess the autonomic nervous system (ANS) activity.After determining the RR interval series for each 5-min R peak time series, the 5-min RR interval segments were resampled by cubic spline interpolation to 4 Hz.We also detrended the segments by the smoothness priors method with λ = 500 [22] to minimize the slow nonstationary trends in the segments.The corresponding power spectral densities (PSDs) for each 5-min RR interval segment were then estimated using Welch's method (eight sections, 50% overlap, Hamming window).Finally, we calculated the frequency-domain HRV parameters from PSDs separately for each 5-min segment: the absolute HF power and the normalized power in the high-frequency band (HF NU = HF/(HF+LF)) to assess the vagal dominance, the ratio of powers in low and high-frequency bands (LF/HF ratio) to assess the sympathetic dominance, and the peak frequency in the HF band that normally matches the respiratory frequency [23], [24].HF and LF represent high-frequency (0.15-0.4 Hz) and low-frequency (0.04-0.15 Hz) bands, respectively.Finally, we calculated the medians of HRV parameter values and the frequencies of the HF peak separately for the DesSev groups.
2) Respiration: Using Welch's method (eight sections, 50% overlap, Hamming window), we also computed the respiratory PSDs from nasal pressure (sampling frequency of 128 Hz) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
for each 5-min segment corresponding to the RR interval segment (Fig. 2) to study the respiratory effect.From the respiratory PSDs, we determined the LF/HF ratio to measure LF modulations unrelated to normal respiration during sleep (e.g., periodic breathing), and the peak frequency in the HF band separately for the DesSev groups.The motivation to compute these frequency-domain parameters also for nasal pressure was to demonstrate the effect of the desaturations on respiration and the interconnectivity of HRV, CRC, and respiration.These parameters were not computed to assess ANS modulation but to aid in the interpretation of the HRV spectrum as respiration affects HRV through RSA.In addition, we determined the mean respiratory frequency (mFr) in each 5-min segment using a peak-conditioned spectral averaging method [25].The respiratory HF-band peak was compared with the corresponding HRV HF-band peak by calculating the frequency difference between the peaks (ΔHF peak = frequency of respiratory HF peakfrequency of HRV HF peak).
Next, the influence of DesSev on cardiorespiratory modulation was quantified as the spectral coherence between the nasal pressure and RR interval segments.Before the coherence analysis, the nasal pressure signal was downsampled to the sampling frequency of the RR interval series, 4 Hz.Welch's method (eight sections, 50% overlap, Hamming window) was also employed for computing the magnitude squared coherence function.The coherence in the HF band was computed only when RR intervals and nasal pressure were linearly related within the HF band: RR intervals and nasal pressure were assumed to be linearly related when the value of spectral coherence exceeded a significant level of coupling.A significant level threshold, ρ = 0.6241, was determined based on a surrogate data analysis setting an α = 1% risk that the RR interval series and nasal pressure signals were coupled when real coupling did not exist [17].The threshold was determined as follows: first, the spectral coherence between two 5-min segments of white noise assumed to be uncorrelated was computed.Second, after 2000 iterations, the maximum coherence in each iteration was obtained and sorted to reduce the risk of setting the threshold too low.Third, the 99th percentile was considered as the threshold value.Finally, we computed Pearson's linear correlation between the spectral coherence and DesSev parameter values within the DesSev quartiles.
3) Statistical Analyses: The DesSev groups could contain numerous 5-min segments from the same patient; thus, we assumed that the samples were not statistically independent.Therefore, we used the Wilcoxon signed-rank test to evaluate the statistical significance of differences between the severity groups and between the frequencies of HRV and respiratory HF peaks.As the Wilcoxon signed-rank test is a pair-wise comparison, we computed a total of 5000 randomly chosen permutations of the parameter pairs.Only one parameter was compared between two severity groups at a time.The median p-value from the 5000 permutations was used to represent the statistical difference.Due to the Bonferroni correction and the large sample sizes, the limit for statistical significance level was set to p < 0.001 for the spectral analyses.The limit for statistical significance in Pearson's correlation analysis was p < 0.05.All spectral, correlation, and statistical analyses were performed with MATLAB R2021b (MathWorks Inc, MA, USA).

III. RESULTS
In PSDs derived from the RR interval series, the HF NU band power decreased and the LF/HF ratio increased with increasing hypoxic load (p < 0.001, Table II, Fig. 3(a)).When considering these parameters, the most severe hypoxic load group, DesSev Q4 , differed significantly (p < 0.001) from other severity groups.Moreover, the absolute HF power was significantly increased in DesSev Q4 .However, the frequency of the HF peak did not significantly vary between the severity groups (p > 0.001, Table II, Fig. 3(a)).

TABLE II SPECTRAL CHARACTERISTICS OF HRV AND RESPIRATION
In respiratory PSDs, the HF NU band power decreased and the LF/HF ratio increased from DesSev 0 to DesSev Q4 (p < 0.001, Table II, Fig. 3(b)).All severity groups differed significantly from the reference group, DesSev 0 , and other severity groups (p < 0.001, Table II).The frequency of the HF peak and therefore, the respiratory frequency increased towards more severe groups (p < 0.001, Table II).
The difference between frequencies of HRV and respiratory HF peaks, i.e. the ΔHF peak, increased towards more severe groups (p < 0.001, Table II, Fig. 4).Furthermore, the spectral coherence between RR intervals and nasal pressure in the HF band decreased significantly with more severe desaturations (p < 0.001, Fig. 5(a)).In addition, a strengthening negative correlation was observed between the DesSev parameter values and the spectral coherence towards more severe hypoxic load groups (−0.304, p < 0.05, Fig. 5(b)).

IV. DISCUSSION
In this study, we investigated the effect of the hypoxic load severity on HF spectral components of HRV and nasal pressure signal because the effect of the hypoxic load severity on the CRC has remained unrecognized.The findings support our hypothesis: we observed a diminished vagal modulation and a lower level of HF coupling between HRV and respiration with the more severe hypoxic load when analyzing the PSDs of the RR interval series and corresponding nasal pressure signals.In addition, we found that the spectral coherence between HRV and respiration weakened as the hypoxic load increased.As a potential clinical implication of these novel findings, a more severe hypoxic load may, therefore, lead to increased cardiac stress in OSA patients.
In the current study, we observed a diminished vagal modulation and almost indistinguishable, non-varying peaks on the HF band in the severity groups when analyzing HRV (Table II, Fig. 3(a)).Our current and previous results [15] are in line with previous studies [9], [10] showing attenuated RSA with hypoxemia.The diminished vagal modulation (i.e., sympathetic overdrive) due to increased hypoxic load and attenuated RSA are associated with severe health outcomes such as cardiovascular diseases [3], [8], [10], [26], [27], partially explaining the high cardiovascular morbidity of OSA patients [10].However, the severity of the hypoxic load did not significantly affect the HF peak frequency based on the HRV analysis (Table II, Fig. 3(a)) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.despite the significant change in breathing pattern (Table II, Fig. 3(b)).This finding is counterintuitive because the HF band of HRV usually matches the respiratory frequency and is thus well-known as the respiratory band [23], [28].Therefore, factors other than respiration may contribute to HRV in the HF band of HRV and RSA, such as baroreflex or sleep stage transitions [29].The difference between the respiratory and HRV peak frequencies in the HF band and their IQRs increased with the increasing hypoxic load supporting our hypothesis of a weakening CRC (Table II, Fig. 4).Loop gain could be a potential factor explaining some of these differences [30] as cyclical underand overventilation due to unstable respiratory control may also disrupt the coupling between HRV and respiration.The CRC and RSA are traditionally assessed through the peak in the HF band of the PSD in HRV analysis.However, the study by Morales et al. [16] and the differences between our results in HRV and respiratory spectral analyses suggest that HRV analysis within the HF band alone might not be the optimal method for assessing the effect of the hypoxic load severity on CRC.
The methods to assess RSA or CRC are not standardized; therefore, we complemented the traditional RSA analysis with cross-spectral coherence analyses [16], [17].We found that within the HF band, the spectral coherence between HRV and respiration weakened significantly with the increasing hypoxic load (Fig. 5(a)).In addition, the negative Pearson's linear correlation between the spectral coherence and the hypoxic load increased towards the more severe DesSev groups (Fig. 5(b)).Together with previous studies [9], [12], [14], [15], [30], [31], [32], [33], [34], these findings support our hypothesis.Whereas hypoxemia decreases HF coupling and increases LF coupling [9], the effect of the hypoxic load and its severity has remained obscure.Hypoxemia sensitizes the carotid body promoting the high loop gain and changing the breathing pattern [9], [30], [31], [32] thus, disturbing the CRC in terms of HF coupling.Thomas et al. [12], [14], [33] showed that high HF coupling can be used as a biomarker for stable sleep.In patients with no OSA, stable nonrapid eye movement (NREM) sleep is characterized by high HF coupling and stable respiration and oxygenation whereas LF coupling is more prominent during unstable NREM sleep [33].However, patients with OSA experience weakened HF coupling, disordered breathing, and oxygenation abnormalities throughout different sleep stages [12], [14].The weaker HF coupling of OSA patients can partially be explained by arousals which lead to unstable sleep [33], [35].In addition, increasing hypoxic load increases progressively the sympathetic overdrive which in turn gradually decreases the degree of CRC in healthy subjects [15], [34].This topic, however, requires further research, because we did not consider the sleep structure in the current study.Nevertheless, our results and the previous literature emphasize the need to also quantify the hypoxic load and its severity in addition to the conventional OSA parameters.
Our current results add an interesting insight to the discussion of the hypoxic load, the physiological consequences of its severity, and the risk for cardiovascular diseases.It is known that sympathetic overdrive and diminished vagal modulation increase the risk of cardiovascular diseases [3], [5], [6], [7], [8].As the current results demonstrate, the increasing hypoxic load exacerbates the sympathetic overdrive [15] and weakens CRC.The physiological functions of RSA and CRC have previously been proposed to improve gas exchange [36]; however, Ben-Tal et al. [11] suggested that the function would be to reduce the cardiac workload while maintaining a healthy level of blood gases.The cardiac workload may be increased due to a more severe hypoxic load and weakened CRC in terms of HF coupling and RSA: hypoxic load severity is more strongly associated with heart failure and cardiovascular mortality compared to the AHI [26], [37] and heart failure patients have been shown to exhibit weaker CRC compared to the healthy population [38].Altered CRC may also play a role in the development of hypertension, stroke, and possibly other cardiovascular diseases [39], [40].In addition to the physiological phenomenon, genetic factors may connect impaired CRC and cardiovascular diseases [41].The sympathetic overdrive in cardiac control predisposes the OSA patients, for example, to systemic inflammation, oxidative stress, and cardiac remodeling and can further increase their risk for hypertension, heart failure, and arrhythmias [3], [5], [42], [43].Cardiorespiratory instabilities, such as sympathovagal imbalance and impaired CRC, are characterized by diminished cardiac vagal modulation and, therefore, they can indicate an underlying pathology and even an increased risk for sudden cardiac death [44].For these aspects, the hypoxic load, its severity, and CRC can provide valuable information on the cardiovascular stress related to OSA.In addition, the hypoxic load has been shown to have a stronger association with adverse OSA-related health outcomes compared to the AHI [26], [37], [45], [46].Thus, hypoxic load severity could enhance the severity assessment of OSA alongside the conventional OSA parameters.However, assessing the cardiovascular load patient-wise warrants further studies because we performed our analyses in pooled segment groups divided based on their hypoxic load severity.
In our study, the main limitation is that we did not consider the sleep stages or arousals in the present study.HRV, RSA, and HF coupling are modulated by sleep stage transitions and arousals, as all three tend to increase towards deeper NREM sleep and are disturbed by the arousals [29], [33], [47], [48], [49].We decided to exclude these factors from the analyses because assessing the sleep stage modulation or determining a single sleep stage for an individual 5-min segment would have been complicated.Nevertheless, we excluded the segments consisting of ≥ 50% wake to prevent bias caused by the wake wake periods.Restricting the study approach only to linear analyses is also a limitation.In addition to ANS activity, multiple factors modulate HR resulting in complex dynamics in HRV [29].For example, the duration of phase synchronization between HR and respiration decreases with the severity of OSA [50].Therefore, more complex models, e.g., nonlinear metrics show the potential to be utilized to gain complementary, more detailed information about the effect of hypoxic load severity on the intimate coupling of HRV and respiration.However, cardiorespiratory phase synchronization is highly dependent on the sleep stage [29], [51] and thus, we left this approach for future studies in which we consider the effect of sleep stages.
In addition, including the patients with multiple comorbidities and medications is a limitation because such factors can modulate HRV and respiration [52].However, the list of medications and certain comorbidities was incomplete; thus, controlling these factors in the present study was not feasible.HRV and CRC are known to decrease with age [28], although age did not significantly affect the HRV modulation in our previous study with a similar population [15].In addition, OSA seems to be a stronger factor leading to reduced CRC during sleep rather than age [12], [28].Finally, counting the desaturation into the 5-min segment in which it started may be a limitation.The hypoxic load of the segments can be under or overestimated because the final desaturation in some segments can continue to the next segment.Even though we consider this limitation to have only a minor effect on our conclusions, the immediate physiological response of the desaturation may not be reflected in the segment into which it is accounted.More thorough research is warranted on the effect of the hypoxic load severity on HRV and CRC and how the aspects mentioned above affect this association.

V. CONCLUSION
Our results show that a more severe hypoxic load is associated with decreased HF power in HRV analyses, increased LF/HF ratio, and decreased HF coupling between HRV and respiration, implying diminished vagal modulation and weakened CRC.These phenomena are indicators of impaired cardiac health and recognized risk factors for several cardiovascular diseases and sudden cardiac death.A more severe hypoxic load leads to stronger physiological stress; therefore, a more thorough desaturation analysis could be beneficial when assessing the severity of OSA.In addition, although HRV is an effective and non-invasive measure of the ANS and ECG is routinely recorded during PSG, neither are commonly utilized in contemporary OSA diagnostics.Hypoxic load and CRC analyses could easily be implemented in the OSA diagnostics, and therefore, they could provide valuable information on the OSA-related cardiovascular risk.

Fig. 3 .
Fig. 3. Median power spectral densities (PSDs) derived from (a) RR interval series (HRV PSDs) and (b) nasal pressure signal (respiratory PSDs) in the groups pooled based on the desaturation severities (Des-Sev) within the 5-min segments.Circles present the median frequency of the high-frequency peak.Vertical dashed lines separate the frequency bands used in heart rate variability analyses: the low-frequency band (0.04-0.15 Hz) and the high-frequency band (0.15-0.40 Hz).

Fig. 4 .
Fig. 4. Median difference between the respiratory and HRV peak frequencies in the HF band (red line) with their interquartile ranges (blue boxes) in the groups pooled based on the desaturation severities (Des-Sev) within the 5-min segments.Whiskers correspond to approximately ±2.7 standard deviation and 99% data coverage.Each severity group DesSev Q1-4 differed statistically significantly (p < 0.001) from the reference group DesSev 0 and each other (p < 0.001).Abbreviations: ΔHF Peak = Difference between the respiratory and HRV peak frequencies in the HF band (frequency of respiratory HF peak -frequency of HRV HF peak), HRV = heart rate variability.

Fig. 5 .
Fig. 5. (a) Spectral coherence between RR intervals and nasal pressure at high-frequency band (0.15-0.40 Hz) in the groups pooled based on the desaturation severity (DesSev) in the 5-min segments.Red lines present the median spectral coherences and blue boxes their interquartile ranges.Whiskers correspond to approximately ±2.7 standard deviation and 99% data coverage.Each severity group DesSev Q1-4 differed statistically significantly (p < 0.001) from the reference group DesSev 0 and each other (p < 0.001).(b) The correlation between spectral coherence and DesSev parameter values within the severity groups.* denotes a statistically significant (p < 0.001) correlation coefficient.