Multivariate and Multiscale Complexity of Long-Range Correlated Cardiovascular and Respiratory Variability Series.

Assessing the dynamical complexity of biological time series represents an important topic with potential applications ranging from the characterization of physiological states and pathological conditions to the calculation of diagnostic parameters. In particular, cardiovascular time series exhibit a variability produced by different physiological control mechanisms coupled with each other, which take into account several variables and operate across multiple time scales that result in the coexistence of short term dynamics and long-range correlations. The most widely employed technique to evaluate the dynamical complexity of a time series at different time scales, the so-called multiscale entropy (MSE), has been proven to be unsuitable in the presence of short multivariate time series to be analyzed at long time scales. This work aims at overcoming these issues via the introduction of a new method for the assessment of the multiscale complexity of multivariate time series. The method first exploits vector autoregressive fractionally integrated (VARFI) models to yield a linear parametric representation of vector stochastic processes characterized by short- and long-range correlations. Then, it provides an analytical formulation, within the theory of state-space models, of how the VARFI parameters change when the processes are observed across multiple time scales, which is finally exploited to derive MSE measures relevant to the overall multivariate process or to one constituent scalar process. The proposed approach is applied on cardiovascular and respiratory time series to assess the complexity of the heart period, systolic arterial pressure and respiration variability measured in a group of healthy subjects during conditions of postural and mental stress. Our results document that the proposed methodology can detect physiologically meaningful multiscale patterns of complexity documented previously, but can also capture significant variations in complexity which cannot be observed using standard methods that do not take into account long-range correlations.

1. Introduction

Cardiovascular oscillations are influenced by the combined activity of different physiological regulation processes and, as a consequence, exhibit a rich dynamical structure [1]. Such different processes do not usually work at a single time-scale, but instead operate across different temporal scales, that for example reflect thermoregulatory or neural parasympathetic and sympathetic control. For this reason, different methods have been recently developed to evaluate the ’multiscale complexity’ of cardiovascular oscillations. These methods allow both to characterize the physiological regulatory systems and to extract diagnostic parameters, and thus can have noteworthy clinical implications: in fact, a decrease of dynamical complexity can be related to an impaired capability of the subsystems composing the organism to interact among each other and it has been already proposed as a marker of pathological conditions [2,3]. Among the proposed methods, the one which is likely most popular is the so-called multiscale entropy (MSE), developed by Costa et al. [4]. This method calculates the conditional entropy (CE) of a single time series (usually the heart period, HP) as a function of the time scale of observation; the change of time scale is achieved through averaging consecutive segments of the time series via a procedure that has been lately recognized to correspond to a two-step process consisting of a low-pass filter followed by downsampling [5]. It is worth noting, however, that the initial formulation of multiscale entropy suffered from drawbacks related both to issues due to the filtering and downsampling steps, and to the unsuitability of CE analysis in conditions of data paucity caused by the availability of short time series and by the needs to explore multivariate time series at coarse time scales. Therefore, in the last years, the definition of MSE has been refined to take into account typical requirements of cardiovascular and cardiorespiratory signal analysis, and specifically: (i) to allow the joint calculation of complexity of multiple variables besides HP, for example systolic arterial pressure (SAP) and respiration (RESP) [6]; (ii) to allow the assessment of the complexity of shorter time series, usually few hundred beats long [7,8]. For short-term physiological time series, complexity has been related to the regularity of the temporal patterns observed in the signals, and thus is usually a measure of the unpredictability of the present sample given a limited number of past samples [9]. However, recent studies have recognized the importance of long-range correlations resulting in slowly varying dynamics also for the analysis of short-term complexity [10], and have started to account for these correlations in multiscale entopy-based analysis [8].

In this work, we introduce a novel method to compute multivariate and multiscale complexity of cardiovascular oscillations. This method fits a multivariate time series using a vector autoregressive fractionally integrated (VARFI) model and then provides the multiscale representation of the VARFI parameters using the theory of state space models, allowing to extract from such parameters multiscale and multivariate measures of complexity. This approach presents several advantages if compared to other previous works in the same field [4,5,6,7,8,10], and in particular: (i) the parametric formulation employed permits to work reliably on short time series; (ii) fractional integration allows to take into account not only short-term dynamics, but also the long-range correlations; (iii) the vector formulation permits to compute the overall complexity of multivariate time series, or the complexity of an individual time series when one or more other series are considered. In this work, such approach is applied to HP, SAP and RESP time series measured in healthy subjects monitored in a resting condition and during two types of physiological stress: postural stress provoked by head-up tilt and mental stress induced by mental arithmetics.

The algorithms for multivariate multiscale analysis of physiological time series presented in this work are collected in the MSE-VARFI MATLAB toolbox, which can be freely downloaded from http://www.lucafaes.net/LMSE-MSE_VARFI.html.

2. Methods

2.1. Measures of Complexity in Linear Multivariate Stochastic processes

Considering a dynamical system whose activity is defined by the zero-mean stationary multivariate stochastic process , where each , , denotes a scalar stochastic process , let us define as the dimensional random variable describing the present state of the system, and as the infinite-dimensional vector variable describing its past states. Then, a measure of complexity of the system, typically defined for univariate systems [10] and here extended to the multivariate case, is the entropy rate defined as where is the conditional entropy of the present given the past, with denoting the Shannon Entropy. If the observed process is a jointly distributed Gaussian process, it can be described without loss of generality through a vector linear regression model fed by white and uncorrelated innovations such that, for each , [11]. In such a case, the entropy of the present state of the process and the conditional entropy of the present given the past can be expressed analytically in terms of the covariance of the process, , and the covariance of the innovations, , as [11,12,13]: where stands for matrix determinant. Then, it follows immediately that the complexity of a multivariate linear process with joint Gaussian distribution is given by Equation (2b). In this work we provide an alternative definition of complexity, which includes a normalization of the innovation covariance to the process covariance:Such normalization, which is implicitly implemented in studies assessing the complexity of univariate time series where the series is normalized to unit variance before computing the complexity measure, is formalized here in Equation (3) for multivariate series, and will be fundamental in the definition of multiscale complexity where the process covariance intrinsically changes with the time scale.

The measure of global complexity defined above for a multivariate process can be particularized to the characterization of one of its constituent processes. Specifically, the complexity of a scalar process , , can be defined in an univariate context with respect to its own dynamics only, in a bivariate context with respect to its dynamics and to the dynamics of another scalar process , or in a fully multivariate context with respect to the dynamics of the whole observed vector process [14]. The derivations are based on the knowledge that, in Gaussian processes, a linear parametric representation captures all of the entropy differences that define the conditional entropies [11] and that these entropy differences are related to the partial variances of the present of the target given its past and the past of one or more other processes, intended as variances of the prediction errors resulting from linear regression [13,15]. Specifically, let us denote as and the prediction error of a linear regression of performed respectively on and , and consider that the prediction error of a linear regression of on is the innovation process defined above. Then, denoting the variances of the three prediction errors as , and , the univariate, bivariate and multivariate normalized complexity measures relevant to the process are defined as: where is the diagonal element of .

2.2. Linear Multivariate Stochastic Processes with Long Range Correlations

In this section we present the parametric approach to the description of linear multivariate Gaussian stochastic processes exhibiting both short-term dynamics and long-range correlations. The approach is based on representing the M-dimensional discrete-time, zero-mean and unit variance stochastic process defined in Section 2.1 as a vector autoregressive fractionally integrated (VARFI) process fed by the uncorrelated Gaussian innovations . The VARFI process takes the form [16]: where L is the back-shift operator , ( is the identity matrix of size M) is a vector autoregressive (VAR) polynomial of order p defined by the coefficient matrices , and where , is the fractional differencing operator defined by [17]: with denoting the gamma (generalized factorial) function. The parameter in Equation (5) determines the long-term behavior of each individual process, while the coefficients of the polynomial allow description of the short-term dynamics. Note that the process defined in Equation (5) is a particular case of the broader class of VARFIMA() processes [16], which also contains the class of autoregressive processes VAR(p); here we restrict our analysis to the description of the VARFIMA() process, which we denote as a VARFI() process.

The parameters of the VARFI model (5), namely the coefficients of and the variance of the innovations , are typically obtained from process realizations of finite length first estimating the differencing parameters by means of the Whittle semi-parametric local estimator [17] separately for each individual process , then defining the filtered data , and finally estimating the VAR parameters from the filtered data using the ordinary least squares method to solve the VAR model , with model order p assessed through the Bayesian information criterion [18].

2.3. Multiscale Complexity of VARFI Processes

In this section we report how to compute across multiple temporal scales the complexity measures defined in Section 2.1, under the hypothesis that the analyzed multivariate process is properly described by the VARFI representation provided in Section 2.2. The procedure for multiscale complexity analysis, which extends the approach proposed in Reference [19] to multivariate processes incorporating long-range correlations, is presented here reporting the essential steps and described with more mathematical details in the Appendices. There are three appendices to the main text. Appendix A reports the procedure for computing the coefficients of a finite-order VAR process that approximates an assigned VARFI process. Appendix B recalls the derivations relevant to the identification of multiscale state space (SS) processes defined as filtered and downsampled versions of an assigned VAR process; the parameters of the SS process defined at an assigned time scale are exploited to compute the multivariate complexity measure at that time scale. Appendix C describes how to define restricted SS processes and to rearrange them in order to extract the partial variances needed for the computation of the univariate and bivariate multiscale complexity measures. The derivations described in Appendix B and Appendix C are taken from Refs. [7,8,19,20].

Before implementing the rescaling procedure, we approximate the VARFI process (5) with a finite order VAR process by truncating the fractional integration part at a finite lag q, that is, we perform the following approximation: This allows us to rewrite the VARFI() process as a VAR(m) process, where : where the new coefficients are , with as in Equation (5) and as in Equation (7a). The exact procedure to derive the coefficients of the VAR(m) process is explained in Appendix A.

In order to represent a scalar stochastic process at the temporal scale defined by the scale factor , a two step procedure is typically employed which consists first in filtering the process with a lowpass filter with cutoff frequency , and then downsampling the filtered process using a decimation factor [5,19]. Extending this approach to the multivariate process , we first implement the following linear finite impulse response (FIR): where r denotes the filter order and , where the coefficients of the polynomial , , are the same for all scalar processes and are chosen to set up a lowpass FIR configuration with cutoff frequency . The filtering step transforms the VAR(p+q) process (8) into a VARMA(p+q,r) process with moving average (MA) part determined by the filter coefficients:Then, we exploit the connection between VARMA processes and state space (SS) processes [21] to evidence that the VARMA process (10) can be expressed in SS form as: where is a -dimensional state process, is the SS innovation process, and the vectors and and the matrix can be obtained from and (see Appendix B).

The second step of the rescaling procedure is to downsample the filtered process in order to complete the multiscale representation. This is done exploiting theoretical findings [7,22,23] which allow to describe the filtered SS process after downsampling in the form: Equation (12) provides the SS form of the filtered and downsampled version of the original VARFI() process, and its parameters () can be obtained from the SS parameters before downsampling and from the downsampling factor ; moreover, the variance of the downsampled process, , can be computed analytically from the parameters of the SS model (12) by solving a discrete-time Lyapunov equation (these steps are shown in the Appendix B). The parameters relevant to the computation of complexity at scale are the variance of the downsampled process, , and the variance of the corresponding innovations, . These variances can indeed be combined in a similar way to that of Equation (3) to yield the expression of the complexity of the original process when it is observed at scale :Note that this measure tends to the theoretical value when .

Finally, we show how to compute any partial variance appearing in Equation (4) from the parameters of an SS model in the form of (12), so that the three complexity measures relevant to the scalar process can be computed at any assigned time scale . To do this, we exploit the formulations reported in Refs. [22,23], showing that the partial variance , where the subscript a denotes any combination of indexes , can be derived from the SS representation of the innovations of a submodel obtained removing the variables not indexed by a from the observation equation. Specifically, we need to consider the submodel with state Equation (12a) and observation equation: where the additional subscript denotes the selection of the rows with indices a in a vector or a matrix. This restricted SS model can be rearranged to extract the partial variance from the covariance matrix of its innovations (see Appendix C), so that with this procedure the univariate, bivariate and multivariate normalized complexity measures can be computed inserting the partial variances , in Equation (4a) and Equation (4b), and using in Equation (4c), together with from Equation (A9b). These individual measures tend to the theoretical value when .

3. Application to Cardiovascular Variability Processes

In this section the proposed method is illustrated using cardiovascular and respiratory series: the heart period (HP), systolic arterial pressure (SAP), and respiration (RESP). Several studies report the interaction between the dynamics of these three time series [24,25,26,27,28], which motivates their use in a multivariate context. The variation of heart period, usually referred as heart rate variability (HRV), reflecting cardiovascular complexity and representing the capability of the organism to react to environmental and psychological stimuli, is the most studied variable and main target variable in cardiovascular and cardiorespiratory spontaneous variability [29,30,31]. For this reason, the conditional measures of a single scalar process will focus on the HP time series, as it has been shown that SAP and RESP have an effect (direct or indirectly) on this process [24,25,26,27,28].

3.1. Experimental Protocol

The HP, SAP and RESP time series were measured in a group of 62 healthy subjects ( years old, 37 females) monitored in the resting supine position (SU), in the upright position (UP) reached through passive head-up tilt, in the recovery supine position (SU) and during mental stress induced by mental arithmetics (MA) [32]. During the measurements, the subjects were free-breathing. For each subject and condition, the multivariate process is defined as . The acquired signals were the surface electrocardiogram (ECG), the finger arterial blood pressure recorded noninvasively by the photoplethysmographic method, and the respiration signal recorded through respiratory inductive plethysmography. For each subject and experimental condition, the values of HP, SAP and RESP were measured on a beat-to-beat basis respectively as the sequences of the temporal distances between consecutive R peaks of the ECG, the maximum values of the arterial pressure waveform taken within the consecutively detected heart periods, and the values of the respiratory signal sampled at the onset of the consecutively detected heart periods. The study was approved by Ethical Committee of the Jessenius Faculty of Medicine, Comenius University and all participants signed a written informed consent. A detailed description of the experimental protocol and signal measurement is reported in Ref. [32].

The analysis was performed on segments of at least 400 consecutive points, free of artifacts and deemed as weak-sense stationary through visual inspection, extracted from the time series for each subject and condition. Missing values and outliers were corrected through linear interpolation and, for HP and when possible, erroneous/missing intervals were substituted by pulse intervals measured as the difference in time between two consecutive SAP measurements (). The three time series were normalized to zero mean and unit variance before multiscale analysis.

3.2. Data Analysis

Two different approaches were followed to compute multiscale complexity: (i) the “eVAR” approach, based on pure VAR model identification, that is, performing the whole procedure described in Section 2 after forcing in Equation (5); (ii) the “eVARFI” approach, based on complete VARFI model identification, that is, following the whole procedure described in Section 2 with , , estimated individually from the original time series and considered in the computations. Pursuing these approaches we compare, respectively, (i) the traditional complexity analysis where long-range correlations are neither removed nor modeled, and (ii) the complexity analysis performed by modeling the long-range correlations and considering them together with the short-term dynamics. Such a comparison is meant to infer the role of long-range correlations in contributing to the multiscale complexity and to its variation between conditions. The VARFI model fitting the time series was identified first estimating the fractional differencing parameter , , individually for each time series using the Whittle estimator, then filtering the time series with the fractional integration polynomial truncated at a lag , and finally estimating the parameters of the polynomial relevant to the short-term dynamics through least squares VAR identification. The value of q has to be selected to approximate the VARFI process, which is theoretically of infinite order, with a finite order VAR process. According to previous studies [8,33] q = 50 is an appropriate value for truncating the VARFI process. By increasing q, we can obtain a more precise approximation of the fractional integration part but with a higher computational cost, while a reduced value (and thus an excessive truncation) causes an underestimation of the complexity and the smoothing of the nonmonotonic trends with the timescale [8]. The VAR model order p was selected as the minimum of the Bayesian information criterion (BIC) figure of merit [34]. Then, multiscale complexity measures were computed implementing a FIR lowpass filter of order , for time scales in the range (1,...,30), which corresponds to lowpass cutoff normalized frequencies . The order of the FIR filter determines its selectivity around the cutoff frequency. In this study, an order r = 48 was selected according to previous settings [8]. The changes in complexity related to the oscillatory rhythms are typically evaluated in cardiovascular variability in low-frequency (LF, from 0.04 Hz to 0.15 Hz) and high-frequency (HF, from 0.15 Hz to 0.4 Hz) spectral bands [29]. To account for the different mean heart rate for each subject and condition, multiscale analysis was performed determining the cutoff frequency of the rescaling filter based on the Nyquist frequency in Hz of the HP time series, which for the time scale is given by , where stands for the mean of the HP process. In this work, considering the physiological LF and HF frequency ranges, four cutoff frequencies were chosen to perform the analysis: 0.4 Hz, 0.15 Hz, 0.1 Hz and 0.04 Hz. To obtain the complexity values at such frequencies, linear interpolation was performed on the profiles .

The differencing parameters were estimated individually for each time series in the interval since the VARFI model is stationary for , while it is nonstationary but mean reverting for . As such, the subjects with an estimation of in at least one condition were removed given that the series is possibly non mean reverting; three subjects were removed, so that a total of 59 subjects was considered for the statistical analysis.

3.3. Statistical Analysis

Significant changes in complexity across the pairs of experimental conditions SU vs. UP and SU vs. MA were assessed via a linear mixed-effects model, that is, a linear regression model that incorporates both fixed and random effects [35]. In our case, the fixed-effects (or factors) were condition and scale, while the random-effect was the subject-dependent intercept that allows for the random variation between subjects. Additionally, the interaction between the factors is also considered. The significance of both the effects (fixed and random) and the interaction between fixed effects was assessed by the significance of the corresponding estimated coefficients at the level of . Residuals were checked for whiteness.

To evaluate the changes of interest, estimated marginal means (EMM) [36] were computed for each difference, or contrast, UP-SU and MA-SU, at each frequency of interest (0.4 Hz, 0.15 Hz, 0.1 Hz and 0.04 Hz). A Z-test was applied to check the significance of these differences with a significance level of and Tukey correction for multiple comparisons. The packages lme4 [37] and emmeans [38] of the R software [39] were used to build the model and to compute EMM, respectively.

4. Results

This section presents the results of multiscale analysis performed for the multivariate complexity as defined in Equation (3), as well as for the complexity computed for the scalar process in four different ways: the univariate complexity of HP with respect to its own dynamics, (Equation (4a)), the bivariate complexity of HP with respect to its own dynamics and to the dynamics of SAP, , or RESP, (Equation (4b)), and the multivariate complexity of HP with respect to the dynamics of the whole trivariate process, (Equation (4c)).

Figure 1 presents the median and quartiles across subjects of the multivariate complexity computed for eVAR (first row) and eVARFI (second row) as a function of the time scale , for SU vs. UP (left column) and SU vs. MA (right column). The measure always tends to its theoretical value when computed for eVAR at high values of . In fact, the complexity value is in general lower for eVARFI than for eVAR, and presents more variability across subjects. From a visual inspection of the multiscale patterns one can infer that for eVAR there is a decrease in complexity moving from SU to UP, evident at short time scales (), and that the complexity seems not to change while moving from SU to MA. For eVARFI, the decrease from SU to UP is visible across the whole range of time scales and there seems to be a decrease from SU to MA for scales larger than 4.

Figure 1

Distribution across subjects of the multivariate complexity measure as a function of the time scale for eVAR (first row) and eVARFI (second row), for SU vs. UP (left column) and SU vs. MA (right column).

Figure 2 presents the same as the previous figure but for the univariate (Figure 2a), bivariate (Figure 2b with SAP and Figure 2c with RESP) and multivariate (Figure 2d) complexity measures relevant to the HP time series. Again, eVARFI estimation presents more variability than eVAR, visible across all measures, and does not reach the theoretical value at long time scales. Although the four measures are similar with each other, slight changes occur when the complexity of HP is computed accounting for the dynamics of the other time series. When compared to the univariate measure , the median complexity value decreases for the bivariate measures at lower time scales, being lower with RESP () than with SAP (); the multivariate measure presents even lower median complexity values for lower time scales. For both eVAR and eVARFI, the complexity of HP decreases at short time scales for all four measures. The differences are more subtle during MA, where only a slight decrease in the median is noticeable at very short time scale for eVARFI.

Figure 2

Distribution across subjects of the univariate (a), bivariate (b) with SAP and (c) with RESP, and multivariate (d) complexity measures as a function of the time scale τ for eVAR (first row) and eVARFI (second row), for SU1 vs. UP (left column) and SU2 vs. MA (right column).

Figure 3 reports the distribution of the multivariate complexity measure (blue dots) (values and boxplot) computed at the four predetermined cutoff frequencies of the multiscale filter (0.4 Hz, 0.15 Hz, 0.1 Hz and 0.04 Hz) for each experimental condition. Statistically significant changes () in complexity across the pairs SU vs UP or SU vs MA, as assessed by the estimated marginal means based on the linear mixed-effects model, are marked with *. Comparing SU and UP, the multivariate complexity decreases significantly both for eVAR and eVARFI at 0.4 Hz, while no significant differences are observed at lower cutoff frequencies. Moving from SU to MA, the measure increases significantly for eVAR at all frequencies except 0.04 Hz, but decreases for eVARFI at frequencies 0.1 Hz and 0.04 Hz.

Figure 3

Distribution of the multivariate complexity measure , depicted as boxplot (mean and confidence intervals, yellow filled box; standard deviation, black vertical line) and original values (dots) for selected frequencies ( 0.4 Hz; 0.15 Hz; 0.1 Hz; 0.04 Hz), computed for the four experimental conditions using eVAR (first row) and eVARFI (second row) identification methods. Statistically significant differences between pairs of conditions are marked with an asterisk.

Figure 4 presents the same information as the previous figure but considering the HP only as target process and computing the corresponding univariate (Figure 4a), bivariate (Figure 4b with SAP and Figure 4c with RESP) and multivariate (Figure 4d) complexity measures. We find that, moving from SU to UP, the univariate measure evaluated at a time scale corresponding to the cutoff frequency of 0.4 Hz decreases significantly for both eVAR and eVARFI; this decreased complexity of HP in the UP condition is observed identically when the dynamics of SAP, of RESP, and of both SAP and RESP are included in the conditional entropy measure. Moreover, the inclusion of one or both the other time series makes such a decrease statistically significant also with a cutoff of 0.15 Hz when eVAR analysis is performed. Moving from SU to MA, we observe that only the eVARFI analysis detects statistically significant differences. Specifically, using eVARFI the univariate complexity of HP decreases during MA at time scales compatible with the cutoff frequency of 0.4 Hz, and also with cutoff 0.15 Hz when SAP dynamics (but not RESP dynamics) are considered. For both protocols and using both estimation methods, the analysis at longer time scales (cutoff 0.1 Hz and 0.04 Hz) does not lead to any significant variation in any of the HP complexity measures.

Figure 4

Distribution of the univariate (a), bivariate (b) with systolic arterial pressure (SAP) and (c) with respiration (RESP), and multivariate (d) complexity measures depicted as boxplot (mean and confidence intervals, yellow filled box; standard deviation, black vertical line) and original values (dots) for selected frequencies (ftHz = 0.4 Hz; 0.15 Hz; 0.1 Hz; 0.04 Hz), computed for the four experimental conditions using eVAR (first row) and eVARFI (second row) identification methods. Statistically significant differences between pairs of conditions are marked with an asterisk.

Figure 5 depicts the distribution across subjects of the differencing parameter computed using the Whittle semiparametric estimator for the three time series and for each condition. The differencing parameter is positive for HP and SAP, while it is lower and centered around zero for RESP. The statistical significance of the variations of this parameter moving from SU to UP, or from SU to MA, was tested with the same linear mixed-effects model used for the complexity measures. Results indicate that there are statistically significant differences in the values of d only for SAP, documenting that conditions of stress evoked by UP or MA determine an increase of the differencing parameter estimated for this time series.

Figure 5

Distribution of the long-rang parameter for each of the time series considered (heart period (HP), SAP and RESP) and for the four conditions.

5. Discussion

Building on recent developments regarding linear multiscale time series analysis [7,8,19,20], the present study introduces a novel approach to assess the dynamical complexity of multivariate processes accounting for the presence of short term dynamics and long-range correlations. The adopted linear parametric framework retains the advantage of previous formulations [7,19,20] related to the computational reliability of complexity estimated even at coarse time scales and over relatively short time series (few hundred points), and is here integrated with the incorporation of long-range dynamics (fundamental to a proper evaluation of complexity at coarse time scales [8]) and extended for the first time towards a fully multivariate implementation.

The usefulness of the proposed multivariate and multiscale measure of complexity was demonstrated in practice by the analysis of cardiovascular and cardiorespiratory interactions. In particular, the multivariate measure of complexity allowed us to evaluate the overall impact that the physiological mechanisms related to postural and mental stress have on the joint dynamics of HP, SAP and RESP. Simultaneously, the multiscale assessment of complexity led to elicit the contribution of mechanisms operating across multiple time scales (evaluation at fine time scales, with cutoff frequency of 0.4 Hz encompassing VLF, LF and HF oscillations) and that of mechanisms confined to slower oscillations (evaluation at coarser time scales, with cutoff at 0.15 Hz excluding HF oscillations and cutoffs at 0.1 Hz and 0.04 Hz excluding progressively also LF oscillations). In addition the evaluation of the complexity of HP, either considering its own dynamics only or the dynamics of SAP and/or RESP, was exploited to relate the overall complexity trends reflected in the multivariate measure to variations specifically related to physiologically well-known mechanisms of cardiovascular and cardiorespiratory interactions. To aid interpretation of the results and the discussion of the related physiological mechanisms, we report in Table 1 the statistically significant increase or decrease in complexity observed without (eVAR method) or with (eVARFI method) modeling long-range correlations during head-up tilt (comparison SU vs. UP) or during mental arithmetics (comparison SU vs. MA) at the time scales encompassing the whole spectral content (0.4 Hz), filtering the HF oscillations (0.15 Hz), and filtering in part (0.1 Hz) or completely (0.04 Hz) the LF oscillations. Complexity is linearly dependent on tilt inclination: a steeper inclination of tilt table produces higher degrees of sympathetic activation, with corresponding lower values of entropy-based complexity measures [40]. Studying the trends of the multivariate complexity measure, its decrease with tilt (Figure 3) documents a simplification of the overall dynamics of HP, SAP and RESP. This result can be mainly driven by the less complex dynamics of heart rate variability in the upright position, that we document calculating the complexity of HP (Figure 4a). Indeed it is well known that the dynamics of HP tend to become less complex as a consequence of the rise of LF oscillations and the weakening of HF oscillations related to the shift of the sympathovagal balance towards sympathetic activation and vagal deactivation during head-up tilt [15,40]. The fact that the decrease of multivariate complexity is not statistically significant for cutoff 0.15 Hz and lower suggests that it is related mainly to the weakening of HF oscillations; indeed, HF oscillations of both HP and SAP are known to be blunted with tilt [15,41,42]. In addition, the finding that the decrease in complexity is observed in the same way for VAR and VARFI analysis suggests that the impact of long-range correlations does not change substantially from rest to tilt. Overall, these results suggest that cardiovascular and respiratory oscillations in the LF and VLF band considered together do not alter significantly their complexity in the transition from the supine to the upright body position.

Table 1

Significant differences (p-value < 0.05) between pair of conditions for each measure and frequency. The arrows indicate if the measure increases or decreases from rest to stress.

Measure	Approach	SU_1 → UP				SU_2 → MA
		0.4	0.15	0.1	0.04	0.4	0.15	0.1	0.04
C¯_X	eVAR	↘				↗	↗	↗
	eVARFI	↘						↘	↘
C¯_X_HP|X_HP	eVAR	↘
	eVARFI	↘				↘
C¯_X_HP|X_SAP,X_HP	eVAR	↘	↘
	eVARFI	↘				↘	↘
C¯_X_HP|X_RESP,X_HP	eVAR	↘	↘
	eVARFI	↘
C¯_X_HP|X	eVAR	↘	↘
	eVARFI	↘					↘

On the other hand, the increase of the global complexity with MA observed for eVAR modeling across multiple time scales (cutoffs 0.4, 0.15, 0.1) is in agreement with previous findings relevant to the HP time series [7]. The fact that such increase is not observed for eVARFI indicates that long range correlations have a different impact on the cardiovascular and respiratory dynamics during MA compared with the relaxed condition SU. In particular, eVARFI estimation reports values of the multivariate complexity which are unchanged at low time scales (cutoff 0.4 Hz), tend to decrease at cutoff 0.15 Hz, and decrease significantly at cutoffs 0.1 Hz and 0.04 Hz. We ascribe this lower complexity to a larger contribution of long-range correlations acting in LF and especially VLF bands, which is supported by the fact that the differencing parameter increases significantly, for HP and especially for SAP, during MA. This confirms the regularizing role of long range correlations on physiological dynamics [8,10].

The complexity of HP assessed at lower time scales (cutoff 0.4 Hz) always decreases in the UP position, for all measures (univariate, both bivariate, and multivariate), as shown in Figure 4. This documents the well-known simplification of heart rate variability induced by head-up tilt, which is known to evoke sympathetic activation and vagal withdrawal making the cardiac dynamics more regular [31,40,43,44]. The fact that this result is observed identically for the eVAR and eVARFI approaches indicates that long-range correlations do not impact significantly the evaluation of complexity performed at short time scales. On the other hand, focusing on intermediate time scales for which HF oscillations are cut off ( = 0.15 Hz), we find that the decrease in complexity is lost when considering the HP dynamics individually, but is maintained when SAP and RESP, taken individually or together, are used to reduce the complexity of HP; moreover, this holds for eVAR but not for eVARFI. Taken together, these results suggest that slow oscillations of SAP and/or RESP reduce the complexity of HP in the transition from rest to tilt, and such regularizing action is related to long-range correlations. This finding is compatible with the known increase of cardiovascular interactions during tilt, documented in previous studies using a number of causality measures including information-theoretic ones [13,24,26,27], and also with the existence of slowly varying respiratory patterns that enhance their effect on the variability of HP during postural stress [24]. Here, the presence of cardiorespiratory and cardiovascular interactions is documented indirectly observing the lower values during tilt of the bivariate and multivariate complexity measures compared with the univariate one, through an approach similar to that proposed in Ref. [14].

The complexity of HP tends to decrease also with the mental stress induced by MA. Contrary to the postural stress induced by HUT, the decrease is observed only using the VARFI approach, and not using the VAR. The absence of significant changes in complexity using methods that do not model long-range correlations is in agreement with previous findings [7,14]. Our results document that the decrease in complexity is due to the different impact of long-range correlations, detected only modeling them through the VARFI approach, and again are confirmed by the higher values of the differencing parameter estimated for HP and for SAP (but indeed not for RESP) during MA compared with the relaxed condition (Figure 5). A differencing parameter d = 0 indicates that there are no long-range correlations. The rather small values of d observed for the RESP time series likely reflect the fact that, even if not controlled, the respiratory activity is confined in a narrow band of the frequency spectrum and thus it does not exhibit the slow trends typical of long-range correlated dynamics. Moreover, the changes of the differencing parameter d observed in the upright position for the systolic pressure (and only to a smaller extend for the heart period) document that the sympathetic activation produced by the tilt table inclination may modulate the impact that long range correlations have on the cardiovascular variability. We note also that, when SAP is considered, the significant changes extend to the second cutoff frequency (0.15 Hz), indicating the increasing contribution of SAP long-memory dynamics in reducing the complexity of HP. Therefore, we state the importance of accounting for long-range correlations in the assessment of the changes in the complexity of heart rate variability induced by mental stress. This may have relevance for the practical applications focused on the detection of mental workload or stress [45,46,47,48]).

Future extensions of the present work can be targeted first at investigating methodological aspects, such as the possibility of describing interactions between processes at the level of long-range correlations (possibly with the modeling of cointegration [49]) or the extension of the framework to the direct evaluation of Granger-causal interactions [20] in a multivariate context where long-range correlations are modeled [8]. Regarding applicative contexts, the analysis of multivariate multiscale dynamics is of particular interest in econometrics [50] and in the neurosciences [51,52], where dynamics spanning several temporal scales are commonly observed and multichannel data acquisition is ubiquitous.