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

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 -type scaling in the power spectrum, also known as the 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 fluctuation means that the observed time series has a non-trivial long-range auto-correlation. 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 [8]). Long-range correlations are hallmarks of complex systems, including living organisms [9-11]. To date, several mathematical models have been proposed to generate the 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 non-invasive 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-19]. By contrast, the relation between 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.

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 at integer time steps. In this process, we assume that the cross-correlation function between and is given by a power-law decay around the time lag where the scaling exponent is defined in the range of , is a constant, 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 . 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 is the frequency, and the asymptotic scaling relation holds for [25]. Although this scaling relation holds only for 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.

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 and its mean cross-increment defined as where is defined as a weighted partial sum with , The weights are assumed to satisfy the self-affinity condition (at least asymptotically) The simplest choice that fulfils this condition is ().

Under the assumption of second-order stationarity, the relations among , and is given by where is defined as and is defined as where denotes the imaginary unit. Moreover, the relation between and is given by This approach provides another scaling analysis as follows: where is the scaling exponent. Note that if and are negatively correlated, then takes a negative value.

When and do not have a long-range correlation with each other, drifts through positive and negative values, and the scaling analysis of (equation (2.10)) does not work. Therefore, to check for uncorrelated cases, we evaluate the multi-scale cross-correlation coefficient [31] as follows: where () is When , , and .

The choice of the weights can improve the performance of the scaling analysis. An optimal choice of is derived from higher-order DMCA. The DMCA procedure is summarized as follows: (1) The observed bivariate time series with zero means and finite variances are integrated separately. (2) The smoothly varying trends included in each of () are approximated by a Savitzky–Golay smoothing filter and removed from the integrated time series. (3) The root-mean cross-increment at each scale is estimated by the detrended profile . 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 . The Savitzky–Golay filter has two parameters: the order of the least-squares polynomial and the sub-interval length , where and 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 , an arbitrary polynomial curve with degree can be traced completely. Thus, DMCA using a higher-order 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 is the unit step function, which equals for and for , and are elements of the inverse matrix of For a practical application of DMCA that analyses an observed bivariate time series , we calculate where (). In this estimation, we employ the temporal sample mean operator defined as where is the temporal sequence number, and is the total number of samples. Finally, the slope of the plot of versus 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 and were shown in [20]. In addition, when , we can derive the scaling relations analytically as and Because the relations listed above require that , an estimation of the time lag 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.

Table 1

Properties of th-order DMCA.

m	detrendable polynomial degree	detectable range of α	corrected time scale s~
0	0	0<α<2	s/1.00
2	2	0<α<4	s/1.93
4	4	0<α<6	s/2.74

To determine the time lag , we estimate the detrended cross-covariance of the observed bivariate time series as where () is a time series filtered by a Savitzky–Golay filter with a window length of , and is a value within the scaling range. The choice of an optimal lag is determined by the maximum value of (equation (2.19)).

In the conventional DMCA (and in DCCA), the reciprocal relation does not always hold. Therefore, a correction must be applied to the scale 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 scale (table 1). For instance, holds approximately when the corrected scale is used for the second-order DMCA.

Analysis methods and data

Analysis methods

Here, we briefly summarize the methodology for practical applications of long-range cross-correlation analysis using DMCA.

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 defined as where is defined as , and the weights are given by equation (2.13). Under the assumption of second-order stationarity, the relation among , the auto-correlation and the power spectrum is given by Note that and are the same functions as those of the DMCA. The DMA scaling exponent is estimated as the slope of a log-log plot of versus . Using equation (3.2), we can analytically derive the scaling relations, when the auto-correlation obeys () and when the power spectral density obeys () [22]. In this analysis, indicates a long-range correlation, whereas represents a long-range anti-correlation. When is the unit impulse function or decays exponentially to zero, the scaling behaviour results asymptotically in .

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)). and were estimated for each pair of time series, where the value of was determined based on the peak position in (equation (2.19)). If no dominant peaks of appeared, was set to be zero.

In this study, the group mean of at each was calculated using the group mean of . When 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 . Using the scaling exponent of , 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 has been expressed by the number of data points. In this study, when we analyse a time series having a sampling frequency , the time scale is expressed by in unit of second instead of the number of data points. As shown in table 1, the relation between and depends on the order of DMA and DMCA. For instance, we use in the case of second-order DMA and DMCA.

Coherence analysis

We estimated the squared coherence [32] as a conventional approach for the cross-correlation analysis. The squared coherence of bivariate time series is defined as where is the cross power spectrum between and , and are the power spectra of and , respectively, and denotes Welch’s overlapped-segment averaging estimator for each spectrum. The squared coherence provides a measure of the degree of cross-correlation between the frequency components of and .

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 cross-correlation with the other time series. We estimated the 95% confidence intervals of each mean based on t-statistics.

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.

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 ( ()). 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 (). 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 (). In the cardiovascular dataset, mean SBP and DBP were and , respectively.

Table 2

Baseline characteristics of HRV parameters of each dataset.

	cardiorespiratory data	cardiovascular data
	(mean ± s.d.)	(mean ± s.d.)
mean NN (ms)	.868±66	1017±151
SDNN (ms)	69.5±26.2	85.8±13.2
ln HF (ln ms^2)	6.25±0.89	6.96±0.83
ln LF (ln ms^2)	7.15±0.89	7.55±0.37
ln VLF (ln ms^2)	7.71±0.59	7.92±0.57
LF/HF	2.77±1.04	2.66±2.32

We evaluated long-range correlations between (i) respiration (airflow and IBI) and RRI and (ii) BP (beat-to-beat SBP and DBP) and HRV (figure 1). Representative time series, including the surrogate time series of RRI, are shown in figure 2. The results of the auto- and cross-correlation analyses using DMA and DMCA, respectively, are shown in figures 3 and 4 for the respiration-HRV interactions and in figures 5 and 6 for the BP-HRV interactions. In each figure, the panels in the left two columns show representative results for individual subjects, and the panels in the rightmost column show the group means for all subjects. The individual subject ID corresponds to the file name provided as electronic supplementary material.

Figure 1

Cardio-respiratory (a) and cardiovascular and (b) datasets analysed in this study.

Figure 2

(Left) Representative cardiorespiratory time series: (a) RR interval (RRI), (b) surrogate RRI, (c) respiration airflow,(d) inter-breath interval (IBI). (Right) Representative cardiovascular time series: (e) RRI; (f) surrogate RRI, (g) Beat-to-beat systolic blood pressure (SBP), (h) beat-to-beat diastolic blood pressure (DBP). (Online version in colour.)

Figure 3

Long-range cross-correlation analysis for RRI and respiration airflow. (a–c) Univariate DMA results of RRI (circle) and respiration airflow (triangle) and bivariate DMCA results between RRI and respiration airflow (cross symbol). (d–f) Multi-scale cross-correlation coefficients between RRI and respiration airflow (cross symbol). The left and centre panels each show representative results (means of three recordings) for an individual subject. The rightmost panels show the groups means for all subjects. In the bivariate analysis, we set the lag as . Scaling exponents were estimated in the range of (corresponding to to 270 s). The estimated slope values are indicated in the figure. (Online version in colour.)

Figure 4

Long-range cross-correlation analysis for RRI (circle) and IBI (triangle). Other details are the same as figure 3. (Online version in colour.)

Figure 5

Long-range cross-correlation analysis for RRI and SBP. (a–c) Univariate DMA results of RRI (circle) and SBP (triangle) and bivariate DMCA results between RRI and SBP (cross symbol). (d–f) Multi-scale cross-correlation coefficients between RRI and SBP (cross symbol). The left and centre panels each show representative results (means of three recordings) for an individual subject. The rightmost panels show means for all subjects. In the bivariate analysis, we set the lag as (SBP was 2.8 0.4 s behind RRI). Scaling exponents were estimated in the range of (corresponding to to 270 s). The estimated slopes are indicated in the figure. (Online version in colour.)

Figure 6

Long-range cross-correlation analysis for RRI (circle) and DBP (triangle). In the bivariate analysis, we set the lag as (DBP was behind RRI). Other details are the same as figure 5. (Online version in colour.)

Respiration and heart rate variability

The mean of the univariate DMA scaling exponents for RRI was (circles in figure 3c). Note that corresponds to the fluctuation when . The mean of the univariate DMA scaling exponents for respiration airflow was (triangles in figure 3c). Theoretically, if the respiration airflow exhibits a sinusoidal-like oscillation, the scaling exponent approaches zero asymptotically. However, the scaling exponent appeared to be slightly greater than zero, which indicated a long-range anti-correlation.

We set the lag as in the bivariate DMCA analysis because the primary peaks of were observed around for respiration airflow and not clearly observed for IBI. The relation between the respiration airflow and RRI showed a significant cross-correlation at approximately 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 (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, was not available (not depicted in figure 4a–c).

Blood pressure and heart rate variability

The means of the univariate DMA scaling exponents for RRI and SBP were and , respectively (circles and triangles in figure 5c). Although the bivariate DMCA results showed an increasing trend in the double logarithmic plot of , a scaling behaviour was not evident. Therefore, the estimated slope values are not shown.

The estimated lags between SBP and RRI were distributed around (SBP is 3 s behind RRI). The relation between RRI and SBP revealed a significant cross-correlation around 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 (triangles in figure 6c), which was smaller than that of SBP and close to . Although the bivariate DMCA results exhibited an increasing trend in the double logarithmic plot of , 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.

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 ( Hz).

Figure 7

Coherence analysis results. Mean coherence with RRI over all subjects is described by red solid lines. (a) Respiration airflow (), (b) IBI (), (c) SBP (), (d) DBP (). The grey shaded areas show the 95% confidence intervals for the uncorrelated surrogate time series. (Online version in colour.)

Discussion

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 : (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 are shown in figure 8a,b. Note that -th order DMCA is equivalent to -th order DMCA, where is an even number. The DMCA2 and DMCA4 results of show the consistency for . This consistency indicates that the adverse effect induced by the non-stationary trend is sufficiently attenuated by using the second-order detrending. The corrected time scale makes it possible to confirm the consistency of the peak and valley positions in and .

Figure 8

Comparison of bivariate analyses using different methods and parameters. RRI and SBP time series of subject L shown in figure 5b,e were analysed. (a,b) Results of zeroth-, second- and fourth-order DMCAs (denoted as DMCA0, DMCA2 and DMCA4, respectively) with . (c,d) Second-order DMCA results with , , . (e,f) Results of first-, second- and third-order DCCAs (denoted as DCCA1, DCCA2 and DCCA3, respectively) with . (Online version in colour.)

The second-order DMCA results with different lags are shown in figure 8c,d. In this case, (equation (2.19)) showed the highest peak at . Therefore, we used in the analysis in the previous section. The scale dependence of and changes dramatically depending on . When and , the signs of change depending on the scale. exhibits sudden drops (vertical dotted lines in figure 8c,d) that correspond to the sign changes of . Therefore, the 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.

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 ) corresponding to the very low frequency range (less than ) 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 cross-correlation 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-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.

Long-range correlations are ubiquitous in real-world time series [44]. Hence, Our approach is well-suited to applications in multivariate time-series analysis.