Assessment of long-range cross-correlations in cardiorespiratory and cardiovascular interactions

We propose higher-order detrending moving-average cross-correlation analysis (DMCA) to assess the long-range cross-correlations in cardiorespiratory and cardiovascular interactions. Although the original (zeroth-order) DMCA employs a simple moving-average detrending filter to remove non-stationary trends embedded in the observed time series, our approach incorporates a Savitzky–Golay filter as a higher-order detrending method. Because the non-stationary trends can adversely affect the long-range correlation assessment, the higher-order detrending serves to improve accuracy. To achieve a more reliable characterization of the long-range cross-correlations, we demonstrate the importance of the following steps: correcting the time scale, confirming the consistency of different order DMCAs, and estimating the time lag between time series. We applied this methodological framework to cardiorespiratory and cardiovascular time series analysis. In the cardiorespiratory interaction, respiratory and heart rate variability (HRV) showed long-range auto-correlations; however, no factor was shared between them. In the cardiovascular interaction, beat-to-beat systolic blood pressure and HRV showed long-range auto-correlations and shared a common long-range, cross-correlated factor. This article is part of the theme issue ‘Advanced computation in cardiovascular physiology: new challenges and opportunities’.

We propose higher-order detrending moving-average cross-correlation analysis (DMCA) to assess the long-range cross-correlations in cardiorespiratory and cardiovascular interactions. Although the original (zeroth-order) DMCA employs a simple movingaverage detrending filter to remove non-stationary trends embedded in the observed time series, our approach incorporates a Savitzky-Golay filter as a higher-order detrending method. Because the nonstationary trends can adversely affect the long-range correlation assessment, the higher-order detrending serves to improve accuracy. To achieve a more reliable characterization of the long-range cross-correlations, we demonstrate the importance of the following steps: correcting the time scale, confirming the consistency of different order DMCAs, and estimating the time lag between time series. We applied this methodological framework to cardiorespiratory and cardiovascular time series analysis. In the cardiorespiratory interaction, respiratory and heart rate variability (HRV) showed long-range autocorrelations; however, no factor was shared between 2021 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Introduction
Heart rate variability (HRV) in a healthy individual displays complex fluctuations even during a resting state [1,2]. One of the typical characteristics of HRV is 1/f β -type scaling in the power spectrum, also known as the 1/f fluctuation [3]. The scaling exponent β during usual daytime activities is controlled to one in healthy young adults [4]. Alterations in β are associated with ageing, disease and mortality [2,5,6]. Characterizing and understanding HRV fluctuations are essential in biomedical applications.
The 1/f fluctuation means that the observed time series has a non-trivial long-range autocorrelation. Specifically, the auto-correlation function of the time series (or its increment) exhibits a slow decay that follows a power law [7] (or logarithmic decay when β = 1 [8]). Long-range correlations are hallmarks of complex systems, including living organisms [9][10][11]. To date, several mathematical models have been proposed to generate the 1/f fluctuation time series [12,13]. However, the underlying mechanisms of fluctuations in biological systems remain unclear.
The biomedical significance of HRV has been well-established [1,2]. For instance, as a noninvasive tool for assessment of autonomic nervous system function, frequency-domain indices, high-frequency (HF: 0.15-0.40 Hz) and low-frequency (LF: 0.04-0.15 Hz) band powers, have been widely used. The HF power in the frequency band around 0.25 Hz is primarily associated with the synchronization between HRV and respiration cycle, mediated by vagal parasympathetic activity. Therefore, HF power has been used as a measure of cardiac parasympathetic activity. The LF power in the frequency band around 0.1 Hz is associated with the synchronization between HRV and arterial blood pressure (BP) oscillation known as the Mayer wave. Although the LF component is believed to be modulated by both sympathetic and parasympathetic nerve activities, its association with sympathetic nerve activity is still under debate [14]. To date, the effects of respiration and BP control on HRV have been evaluated through frequency-specific components such as the HF and LF powers [15][16][17][18][19]. By contrast, the relation between 1/f fluctuations and respiratory-cardiovascular interactions has not been investigated. Therefore, we explore this relation to gain new insights into the complex dynamics of HRV.
Recently, we developed higher-order detrending moving-average cross-correlation analysis (DMCA) to detect the existence of a common long-range correlated factor hidden in two different processes [20]. In our method, non-stationary trends embedded in each time series are removed using a Savitzky-Golay detrending filter [21], because such trends adversely affect the scaling exponent estimation [22]. We have established a rigorous mathematical foundation of the higher-order DMCA [20]. In this paper, we use this method to evaluate cardiorespiratory and cardiovascular interactions. We analyse heart rate (HR) and respiration variability to study cardiorespiratory interactions and HR and beat-to-beat BP variability to study cardiovascular interactions.

Long-range cross-correlation
In this section, we briefly illustrate the long-range cross-correlations in bivariate stochastic processes [23] and its characterization using DMCA [20]. The long-range cross-correlation between two different processes implies that the both processes are affected by a common noise or have a common long-range correlated factor. Therefore, the detection of long-range cross-correlations provides insight into how two processes interact.

(a) Basics of long-range cross-correlation
The basic properties of long-range cross-correlated processes are described under an assumption of second-order stationarity. For a sake of simplicity, let us consider a bivariate stochastic process {(X (1) [i], X (2) [i])} at integer time steps. In this process, we assume that the cross-correlation function between X (1) [i] and X (2) [i] is given by a power-law decay around the time lag κ 0 where the scaling exponent γ is defined in the range of 0 < γ < 1, κ 0 is a constant, E[ · ] denotes the expectation over ensembles, and the tilde symbol denotes proportionality. In this paper, we denote random variables by capital letters and their realizations by the corresponding lowercase letters. The power-law decay of the cross-correlation means the nonzero correlation persists over a long range of time lag k. This persistent cross-correlation can be characterized by γ , while the power-law decay of the correlation has no characteristic time scale. The cross-correlation can also be characterized by the cross-power spectral density, because this spectral density is given by the discrete Fourier transform of the cross-correlation function [24]. Under the assumption of equation (2.1), the corresponding cross-power spectral density is given by where f is the frequency, and the asymptotic scaling relation β = 1 − γ holds for 0 < β < 1 [25]. Although this scaling relation holds only for 0 < γ < 1 due to the stationarity condition, the range of the scaling exponent β is not restricted to this range. Therefore, a scaling analysis using the cross-power spectral density can be applied to diffusive non-stationary processes such as bivariate fractional Brownian motions [26]. However, a spurious scaling exponent can emerge from non-stationary smooth trends embedded in time series [20]. Therefore, detrending-operation-based approaches, such as detrended cross-correlation analysis (DCCA) [27] and DMCA [20,28], have been introduced to achieve reliable estimations of the scaling exponent. DCCA is one of the most widely used methods to detect long-range cross-correlations in real-world time series. However, DMCA has clear advantages over DCCA. For example, DMCA is mathematically simple and has a well-established formulation (see below for details); it also performs detrending with a linear frequency response [29] and a fast implementation algorithm [30]. Consequently, we employ DMCA in our analysis.

(b) Detrending moving-average cross-correlation analysis
The mathematical foundation of DMCA is established using the framework of a generalized cross random walk analysis (CRWA) [20]. In this framework, we consider a bivariate time series {(X (1) [i], X (2) [i])} and its mean cross-increment defined as Under the assumption of second-order stationarity, the relations among F 12 (s, κ), C 12 (k) and |S 12 (f )| is given by where L(k, s) is defined as and G s (f ) is defined as where i denotes the imaginary unit. Moreover, the relation between |G s (f )| 2 and L(k, s) is given by This approach provides another scaling analysis as follows: where α is the scaling exponent. Note that if X (1) and X (2) are negatively correlated, then F 2 12 (s, κ) takes a negative value.
When X (1) [i] and X (2) [i] do not have a long-range correlation with each other, F 2 12 (s, κ) drifts through positive and negative values, and the scaling analysis of F 12 (s, κ) (equation (2.10)) does not work. Therefore, to check for uncorrelated cases, we evaluate the multi-scale cross-correlation coefficient [31] as follows: , (2.11) where F l (s) (l = 1, 2) is . The higher-order DMCA employs a Savitzky-Golay smoothing filter to remove non-stationary trends more effectively than the original (zeroth-order) DMCA, which uses a simple moving-average filter. The filtered time series are obtained by the values of the least-squares polynomial at the centre point of each sub-interval with a window length s. The Savitzky-Golay filter has two parameters: the order of the least-squares polynomial m and the sub-interval length s, where m and s are restricted to even and odd numbers, respectively. Note that the zeroth-order Savitzky-Golay filter is equivalent to the simple moving-average filter. Using the Savitzky-Golay smoothing filter with an even polynomial order m, an arbitrary polynomial curve with degree (m + 1) can be traced completely. Thus, DMCA using a higherorder Savitzky-Golay filter can well eliminate non-stationary trends locally approximated by a higher degree polynomial.  In the higher-order DMCA, the weight in equation (2.4) is given by where Θ(x) is the unit step function, which equals 0 for x < 0 and 1 for x ≥ 0, and b i,j (s) are elements of the inverse matrix of (2.14) For a practical application of DMCA that analyses an observed bivariate time series 2). In this estimation, we employ the temporal sample mean operator defined as where i is the temporal sequence number, and N is the total number of samples. Finally, the slope of the plot of log F 12 (s, κ) versus log s provides an estimate of α.
The methodological properties of DMCA are summarized in table 1. The detrendable polynomial degree is the upper limit of the polynomial degree of the trends embedded in the observed time series before integration. Analytical forms of kernels L(k, s) and |G s (f )| 2 were shown in [20]. In addition, when κ = κ 0 , we can derive the scaling relations analytically as Because the relations listed above require that κ = κ 0 , an estimation of the time lag κ 0 is required prior to performing a DMCA. Note that a lag estimation is always required for the time-domain analyses of long-range cross-correlations, such as the cross-correlation function and DCCA.
To determine the time lag κ, we estimate the detrended cross-covariance of the observed bivariate time series {(x (1) is a time series filtered by a Savitzky-Golay filter with a window length of s, and s is a value within the scaling range. The choice of an optimal lag k is determined by the maximum value of |K 12 (k)| (equation (2.19)). In the conventional DMCA (and in DCCA), the reciprocal relation s ≈ 1/f does not always hold. Therefore, a correction must be applied to the scale s to confirm the consistency of the detected characteristic scale in the domains of time and frequency. As we pointed out in [20], the time-frequency relation is distorted because of the band-pass filtering property of the detrending procedure in DMCA. Consequently, we use a corrected time scales (table 1). For instance,s ≈ 1/f holds approximately when the corrected scales = s/1.93 is used for the second-order DMCA.

Analysis methods and data (a) Analysis methods
Here, we briefly summarize the methodology for practical applications of long-range crosscorrelation analysis using DMCA.

(i) Long-range auto-and cross-correlation analysis
We first assessed the long-range auto-correlation of each univariate time series. We employed the second-order detrending moving average analysis (DMA) [22], which is the univariate version of the second-order DMCA. In DMA, we consider an estimator of the root-mean square deviation for the observed time series {x[i]} defined as , and the weights {w s (j)} are given by equation (2.13). Under the assumption of second-order stationarity, the relation among F(s), the auto-correlation C(k) and the power spectrum S(f ) is given by Note that L(k, s) and |G s (f )| 2 are the same functions as those of the DMCA. The DMA scaling exponent α is estimated as the slope of a log-log plot of F(s) versus s. Using equation (3.2), we can analytically derive the scaling relations, α = 1 − γ /2 when the auto-correlation obeys C(k) ∼ k −γ (0 < γ < 1) and α = (β + 1)/2 when the power spectral density obeys S(f ) ∼ f −β (−1 < β < 2m + 3) [22]. In this analysis, 0.5 < α < 1 indicates a long-range correlation, whereas 0 < α < 0.5 represents a long-range anti-correlation. When C(k) is the unit impulse function or decays exponentially to zero, the scaling behaviour results asymptotically in α = 0.5. In our long-range cross-correlation analysis, the second-order DMCA was used because the consistencies of the second-and fourth-order DMCA results were confirmed for our datasets (see §5(a)). F 2 12 (s, κ) and ρ 12 (s, κ) were estimated for each pair of time series, where the value of κ was determined based on the peak position in |K 12 (k)| (equation (2.19)). If no dominant peaks of |K 12 (k)| appeared, κ was set to be zero.
In this study, the group mean of F 12 (s, κ) at each s was calculated using the group mean of F 2 12 (s, κ). When ρ 12 (s, κ) consistently exhibits non-zero values (positive or negative) over the scaling range, it suggests that a common noise source exists. The long-range cross-correlation of the common noise source is characterized by the scale dependence of F 12 (s, κ). Using the scaling exponent α of F 12 (s, κ) ∼ s α , we can characterize the long-range correlated behaviour of the common noise source. The interpretation of α is the same as that of univariate DMA.
In the original DMA and DMCA, the time scale s has been expressed by the number of data points. In this study, when we analyse a time series having a sampling frequency f samp , the time scale is expressed bys/f samp in unit of second instead of the number of data points. As shown in

(ii) Coherence analysis
We estimated the squared coherence [32] as a conventional approach for the cross-correlation analysis. The squared coherence of bivariate time series {(x (1) [i], x (2) [i])} is defined as where S 12 (f ) is the cross power spectrum between x (1) [i] and x (2) [i], S 1 (f ) and S 2 (f ) are the power spectra of x (1) [i] and x (2) [i], respectively, and · denotes Welch's overlapped-segment averaging estimator for each spectrum. The squared coherence provides a measure of the degree of crosscorrelation between the frequency components of x (1) [i] and x (2) [i].

(iii) Surrogate time series
We generated sets of surrogate time series without any cross-correlation and estimated the 95% confidence intervals for an uncorrelated case to determine the statistical significance of the observed cross-correlation between two time series. The surrogate time series were generated using an amplitude adjusted Fourier transform (AAFT) algorithm [33]. We initially generated a stochastic time series with the same power spectrum as the original time series using a Fourier phase randomization. The original time-series values were then sorted based on the order of the phase-randomized time series. We produced 100 surrogate time series with the AAFT algorithm that had the same distribution and a similar power spectrum as the original series, but no crosscorrelation with the other time series. We estimated the 95% confidence intervals of each mean based on t-statistics.

(iv) Conventional HRV parameters
To provide the basic characteristics of HRV,the following HRV parameters were calculated: mean RR intervals (mean NN), standard deviation (s.d.) of RR intervals (SDNN), the variances corresponding to very low frequency (VLF; 0-0.04 Hz), low frequency (LF; 0.04-0.15 Hz), and high frequency (HF; 0.15-0.40 Hz) bands, and LF/HF ratio, all of which were proposed by the Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology [1]. The variances of these frequency components were transformed to natural logarithmic values.

(b) Data
We analyse two datasets obtained from healthy human subjects at rest. The first dataset consists of breathing airflow signals and HRV measurements. Electrocardiography (ECG) and nose-mouth airflow signals were simultaneously recorded in eight male subjects (22.4 ± 0.9 yr (mean ± s.d.)). Each subject was seated in a relaxed position for the three recordings (30 min each). Data were recorded with PowerLab and LabChart Pro software (AD Instruments) at a sampling frequency of 1000 Hz. Breathing airflow measurements were recorded with a thermistor flow sensor (NIHON SANTEKU Co., Ltd) attached to the subject's nose and mouth. RR interval (RRI) time series derived from the ECG recordings were linearly interpolated and resampled at 2 Hz. The airflow signals were also resampled at 2 Hz. In addition, inter-breath intervals (IBI) were extracted by detecting the peaks and valleys in the airflow signals, linearly interpolated, and resampled at 2 Hz. The second dataset consists of beat-to-beat systolic BP (SBP), diastolic BP (DBP), and RRI. ECG and continuous BP were simultaneously recorded at a sampling frequency of 1000 Hz in four male subjects (22.3 ± 1.0 yr). A tonometric device (JENTOW7700, Colin Electronics, Japan) was used to measure continuous BP from the right wrist. To synchronize ECG and continuous BP, both signals were connected to PowerLab and LabChart Pro software (AD Instruments). Data were recorded for 20 min after resting for 5 min in a supine position on a bed. Beat-to-beat SBP and DBP were extracted from the peaks and valleys in the continuous BP signal. The RRI, SBP and DBP time series were then resampled at 2 Hz.
All erroneous peak detections in RRI, IBI, SBP and DBP time series were manually corrected based on the visual inspection of the signals.

Results
The baseline characteristics of conventional HRV parameters in each dataset are shown in table 2. In the cardiorespiratory dataset, mean IBI was 3.8 ± 0.5 s (mean ± SD). In the cardiovascular dataset, mean SBP and DBP were 108.6 ± 7.2 mmHg and 50.9 ± 2.4 mmHg, respectively.
We evaluated long-range correlations between (i) respiration (airflow and IBI) and RRI and

(a) Respiration and heart rate variability
The mean of the univariate DMA scaling exponents for RRI was α RRI = 0.93 ± 0.02 (circles in figure 3c). Note that α RRI = 1 corresponds to the 1/f fluctuation when β = 1. The mean of the univariate DMA scaling exponents for respiration airflow was α airflow = 0.17 ± 0.02 (triangles in figure 3c). Theoretically, if the respiration airflow exhibits a sinusoidal-like oscillation, the scaling exponent approaches zero asymptotically. However, the scaling exponent α airflow appeared to be slightly greater than zero, which indicated a long-range anti-correlation.
We set the lag as κ = 0 in the bivariate DMCA analysis because the primary peaks of K 12 (k) were observed around k = 0 for respiration airflow and not clearly observed for IBI. The relation between the respiration airflow and RRI showed a significant cross-correlation at approximately 4 s in the HF band (figure 3d-f ). The bivariate DMCA result (figure 3c) suggested that the RRI was weakly affected by the fluctuations associated with the respiration airflow in a long-range anti-correlated manner. The common factor between the respiration airflow and RRI behaved like the respiration airflow.
The mean of the univariate DMA scaling exponents for IBI was α IBI = 0.84 ± 0.01 (circles in figure 4c). Although RRI and IBI exhibited long-range correlation, no common factor was observed. RRI and IBI were not cross-correlated over the entire time scale (figure 4d-f ). Therefore, F 12 (s, κ) was not available (not depicted in figure 4a-c).

(b) Blood pressure and heart rate variability
The means of the univariate DMA scaling exponents for RRI and SBP were α RRI = 1.01 ± 0.03 and α SBP = 0.84 ± 0.02, respectively (circles and triangles in figure 5c). Although the bivariate DMCA results showed an increasing trend in the double logarithmic plot of F 12 , a scaling behaviour was not evident. Therefore, the estimated slope values are not shown.
The estimated lags between SBP and RRI were distributed around κ = 6 (SBP is 3 s behind RRI). The relation between RRI and SBP revealed a significant cross-correlation around 10 s in the LF band (figure 5d-f ). The results of each subject consistently suggested that a long-range cross-correlated factor between RRI and SBP existed. However, its correlation was weak, and the scaling law was not evident. The comparison with the surrogate RRI time series reveals that the significance of our observations was marginal (figure 5f ).   The mean of the univariate DMA scaling exponents for DBP is α DBP = 0.57 ± 0.02 (triangles in figure 6c), which was smaller than that of SBP and close to 0.5. Although the bivariate DMCA results exhibited an increasing trend in the double logarithmic plot of F 12 , a scaling behaviour was not evident (like the relation between RRI and SBP). Therefore, the estimated slope values are not shown.
The existence of a long-range cross-correlated factor between RRI and DBP was suggested (as it is between RRI and SBP), although its significance compared with the surrogate RRI time series was less than that of RRI and SBP.

(c) Coherence analysis results
The group means of the estimated squared coherence are shown in figure 7. The significantly cross-correlated fluctuations between respiration airflow and RRI were observed primarily in the HF band (figure 7a). Uncorrelated fluctuations between IBI and RRI were observed across all frequencies (figure 7b) and were consistent with the bivariate DMCA result (figure 4f ). The significantly cross-correlated fluctuations between SBP and RRI (figure 7c) and between DBP and RRI (figure 7d) were observed in the HF and LF bands. In both cases, the peak frequencies of the squared coherence in the LF band were consistent with the bivariate DMCA results (figure 5f and 6f ). In all cases shown in figure 7, no significant cross-correlations were observed in the very low-frequency band (f < 0.04 Hz).   log 10 (s/f samp ) log 10 (s/f samp ) log 10 (s/f samp ) log 10 F 12 (s) log 10 F 2 (s) log 10 F 1 (s) log 10 F 12 (s) log 10 F 2 (s) log 10 F 1 (s) log 10 F 12 (s) log 10 F 2 (s)

Discussion (a) Methodological notes on detrending moving-average cross-correlation analysis
An important scientific contribution of our study is the establishment of a methodological framework for long-range cross-correlation analyses using DMCA. The following procedures are required to achieve a reliable and valid estimate of the long-range cross-correlation quantified by F 12 (s, κ): (i) Time-scale correction in the higher-order DMCAs; (ii) consistency checks for the different order DMCAs; (iii) estimation of the time lag between two time series. These points have not been explicitly considered in previous studies using DMCA.
We illustrate the importance of the points mentioned above in figure 8, which shows bivariate DMCA and DCCA analyses of a single subject's data consisting of RRI and SBP time series. The results obtained from the zeroth-, second-and fourth-order DMCAs (denoted as DMCA0, DMCA2 and DMCA4, respectively) with κ = 6 are shown in figure 8a,  log 10 (s/f samp ) log 10 (s/f samp ) log 10 (s/f samp ) log 10 F 12 (s) DMCA is equivalent to (m + 1)-th order DMCA, where m is an even number. The DMCA2 and DMCA4 results of ρ 12 show the consistency fors/f samp > 0.7. This consistency indicates that the adverse effect induced by the non-stationary trend is sufficiently attenuated by using the secondorder detrending. The corrected time scales makes it possible to confirm the consistency of the peak and valley positions in F 12 (s, κ) and ρ 12 (s, κ). The second-order DMCA results with different lags κ = 6, 0, −6 are shown in figure 8c,d. In this case, |K 12 | (equation (2.19)) showed the highest peak at κ = 6. Therefore, we used κ = 6 in the analysis in the previous section. The scale dependence of F 12 (s, κ) and ρ 12 (s, κ) changes dramatically depending on κ. When κ = 0 and κ = −6, the signs of ρ 12 (s, κ) change depending on the scale. F 12 (s, κ) exhibits sudden drops (vertical dotted lines in figure 8c,d) that correspond to the sign changes of ρ 12 (s, κ). Therefore, the F 12 (s, κ) does not provide any meaningful information without the appropriate lag.
We analysed the same data using DCCA (for methodological details see [27]) to compare DMCA and DCCA (figure 8e,f ). In the original DCCA, the corrected scales and the lag estimation are not introduced. Therefore, the DCCA results are very different from our DMCA results. Moreover, the detrending operation causes the DCCA results to exhibit non-smooth behaviours. In the DCCA, non-stationary trends included in a time series are eliminated using a piecewise least-squares polynomial fitting. This procedure acts as nonlinear filtering and induces high-frequency fluctuations in the detrended time series [34]. Therefore, such high-frequency fluctuations cause non-smoothness in the DCCA results.

(b) Cardio-respiratory and cardiovascular interactions
An oscillatory component observed in the HF range of HRV is modulated by respiration through a process known as the respiratory sinus arrhythmia (RSA) [35]. RSA is modulated by inhibition of cardiac parasympathetic activity during inspiration and return of the parasympathetic activity during expiration. To date, cardiorespiratory interactions have been studied mainly focusing on the oscillation component associated with RSA. Using conventional linear analysis methods such as coherence analysis, the significant interactions between respiration airflow and RRI were observed only in the HF range (figure 7a). By contrast, in our analysis using DMCA, the significant interactions between respiration airflow and RRI were observed in wider time scales including HF and LF bands (figure 3f ). Therefore, our approach could provide a new and important insight into nonlinear interactions in the cardiorespiratory system. Because no correlations were observed between IBI and RRI (figure 4f ), our findings suggest that the cardiorespiratory interactions are mainly modulated by the respiratory amplitude but not by the respiratory rate. To date, complex cardiorespiratory interactions been modelled and analysed using approaches in statistical and nonlinear physics [36,37]. Our approach could facilitate to understand the cardiorespiratory interactions.
The respiration variability described by IBI exhibits long-range auto-correlations. Although long-range correlations in IBI have been reported [38], the interaction between IBI and RRI has not been studied from the perspective of long-range cross-correlation. Our analysis revealed that there was no cross-correlated factor in IBI-RRI interactions. This finding did not help us understand the origin of the long-range auto-correlations of RRI. However, an uncorrelated relation between IBI and RRI suggests that respiration variability may have biomedical significance and applications independent of the well-established features of HRV.
Our results (figures 5 and 6) showed that SBP exhibits long-range auto-correlations and that SBP and RRI share a common long-range correlated factor. Strong negative correlations observed in the LF range (figure 5f ) would be associated with Mayer wave oscillation. However, the common factor observed in the longer time scales (greater than 25 s) corresponding to the very low frequency range (less than 0.04 Hz) is a random fluctuation, because smooth sinusoidal waves are eliminated by the detrending procedure of higher-order DMCA. The estimated lags between RRI and BP indicates that RRI precedes BP. Therefore, the primary source of a long-range crosscorrelation should be RRI. Previous studies have demonstrated that the Starling law and diastolic runoff could explain the causal direction of the interaction from HRV to SBP in the LF band [39]. Therefore, the long-range cross-correlation between HRV and SBP may have the same mechanical origin. Previous studies have proposed that the long-range auto-correlated BP variability has medical significance [40][41][42]. Therefore, the long-range cross-correlation between BP and HRV may represent a new medical bio-marker [43].
Definitely, our findings are still preliminary and require further validation. We tested data only from healthy young male adult subjects under resting condition. Thus, characteristics of the long-range cross-correlation, such as body posture, age and sex dependence, are not fully evaluated. To understand the characteristics and mechanism of cardiorespiratory and cardiovascular interactions, a more extensive and more systematic study is required.

Conclusion
We illustrated the practical usage of higher-order DMCA with a Savitzky-Golay filter and demonstrated its ability to detect long-range cross-correlations. The following three points are essential to achieve a reliable characterization of long-range cross-correlations using DMA and DMCA: (i) time-scale corrections in higher-order DMCAs, (ii) consistency checks of different order DMCAs and (iii) time-lag estimations.
Our approach provided new insights into cardiorespiratory and cardiovascular interactions. IBI and RRI showed long-range auto-correlations in cardiorespiratory interactions but did not share any common factor in the analysed scale. This finding suggests that the cardiorespiratory interactions are mainly modulated by the respiratory amplitude but not by the respiratory rate. In cardiovascular interactions, SBP and RRI showed long-range auto-correlations. In addition, they shared a common factor in the analysed scale, although the scaling law of the cross-correlation was not evident.