Multiscale partial information decomposition of dynamic processes with short and long-range correlations: theory and application to cardiovascular control

In an information-theoretic framework, the directed transfer of information between scalar sub-processes is assessed by the Transfer Entropy (TE). Specifically, the TE quantifies the amount of information that the past of the source provides about the present of the target process over and above the information already provided by the past of the target itself (Schreiber 2000). Transfer entropy between the source process component i and the target component j is defined as

$$\begin{eqnarray}&&{T}_{i\to j}=I\left({X}_{j,n};{X}_{i,n}^{-}| {X}_{j,n}^{-}\right),\end{eqnarray} \tag{ 1 }$$

where I( · ; · ∣ · ) denotes conditional Mutual Information (MI) (Cover and Thomas 2005). Considering two sources X_i and X_k and a target X_j , the information transferred towards X_j from the sources X_i and X_k taken together is quantified by the joint TE (JTE)

$$\begin{eqnarray}&&{T}_{{ik}\to j}=I\left({X}_{j,n};{X}_{i,n}^{-},{X}_{k,n}^{-}| {X}_{j,n}^{-}\right).\end{eqnarray} \tag{ 2 }$$

Generally, the joint TE differs from the sum of the two individual TEs, since the source processes (i, k) typically interact with each other while they transfer information to the target process, j. The joint TE (2) can be decomposed under an Interaction Information Decomposition (IID) framework as (Faes et al 2017)

$$\begin{eqnarray}&&{T}_{{ik}\to j}={T}_{i\to j}+{T}_{k\to j}+{I}_{{ik}\to j},\end{eqnarray} \tag{ 3 }$$

where I_ik→j is denoted as Interaction Transfer Entropy (ITE) since it is equivalent to the interaction information (McGill 1954) computed between the present of the target and the past of the two sources, conditioned on the past of the target

$$\begin{eqnarray}&&{I}_{{ik}\to j}=I\left({X}_{j,n};{X}_{i,n}^{-};{X}_{k,n}^{-}| {X}_{j,n}^{-}\right).\end{eqnarray} \tag{ 4 }$$

The ITE quantifies the modification of the information transferred from the source processes X_i and X_k to the target X_j . The ITE can take positive and negative values. Positive values of I_ik→j denote synergy, where the joint TE is greater than the sum of the two individual TEs $\left({T}_{{ik}\to j}\gt {T}_{i\to j}+{T}_{k\to j}\right)$. In contrast, negative values of I_ik→j refer to redundancy, which occurs when the information transferred from the sources to the target overlapped, meaning that the sum of individual TEs is larger than the joint TE $\left({T}_{i\to j}+{T}_{k\to j}\gt {T}_{{ik}\to j}\right)$.

The main drawback of IID is that the interaction TE is quantified using only one measure, and thus makes redundancy and synergy mutually exclusive. This disadvantage can be overcome by the Partial Information Decomposition (PID) (Williams and Beer 2010) encompassing four distinct positive quantities (Faes et al 2017)

$$\begin{eqnarray}&&{T}_{{ik}\to j}={U}_{i\to j}+{U}_{k\to j}+{R}_{{ik}\to j}+{S}_{{ik}\to j},\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{T}_{i\to j}={U}_{i\to j}+{R}_{{ik}\to j},\end{eqnarray} \tag{ 5b }$$

$$\begin{eqnarray}&&{T}_{k\to j}={U}_{k\to j}+{R}_{{ik}\to j}.\end{eqnarray} \tag{ 5c }$$

The terms U_i→j and U_k→j quantify the parts of the information transferred to the target process X_j , which are unique to the source processes X_i and X_k , respectively, thus reflecting contributions to the predictability of the target that can be obtained from one of the sources alone. Then, the terms R_ik→j and S_ik→j quantify the redundant and synergistic interaction between the two sources and the target, respectively.

When compared to IID (3), the PID (5) has the advantage that it provides distinct non-negative measures of redundancy and synergy, therefore allowing the simultaneous presence of redundancy and synergy as distinct elements of information modification. Moreover, the IID and PID are related to each other as

$$\begin{eqnarray}&&{I}_{{ik}\to j}={S}_{{ik}\to j}-{R}_{{ik}\to j},\end{eqnarray} \tag{ 6 }$$

thus showing that the interaction TE is actually a measure of the 'net' synergy manifested in the transfer of information from the two sources to the target (Krohova et al 2019).

The main issue with the PID (5) is that its constituent measures cannot be obtained from the classic information theory simply subtracting conditional MI terms, as done for the IID. Therefore, to complete the PID an additional ingredient to the theory is needed to get a fourth defining equation to be added to (5) for providing an unambiguous definition of U_i→j , U_k→j , R_ik→i and S_ik→j . Several PID definitions have been proposed arising from different conceptual definitions of redundancy and synergy (Harder et al 2013, Bertschinger et al 2014, Griffith et al 2014). Here we make reference to the so-called Minimum Mutual Information PID (MMI-PID) (Barrett 2015). In this approach, redundancy is defined as the minimum of the information provided by each individual source to the target. This leads to the following definition of the redundant TE

$$\begin{eqnarray}&&{R}_{{ik}\to j}=\min \left\{{T}_{i\to j},{T}_{k\to j}\right\}.\end{eqnarray} \tag{ 7 }$$

This definition satisfies the desirable property that the redundant TE is independent of the correlation between the source processes. Furthermore, if the observed processes have a joint Gaussian distribution, all previously-proposed PID formulations reduce to the MMI PID (Barrett 2015).

The classic parametric approach to describe linear Gaussian stochastic processes exhibiting both short-term dynamics and long-range correlations is based on representing an M-dimensional discrete-time, zero-mean and unitary variance stochastic process X_n as a Vector Autoregressive Fractionally Integrated (VARFI) process fed by uncorrelated Gaussian innovations E_n . The VARFI(p, d) process is expressed as (Tsay 2010)

$$\begin{eqnarray}&&{\bf{A}}(L)\mathrm{diag}\left({{\rm{\nabla }}}^{{\bf{d}}}\right){{\bf{X}}}_{n}={{\bf{E}}}_{n},\end{eqnarray} \tag{ 8 }$$

where L is the back-shift operator $\left({L}^{i}{{\bf{X}}}_{n}={{\bf{X}}}_{n-i}\right)$, ${\bf{A}}(L)={{\bf{I}}}_{M}-{\sum }_{i\,=\,1}^{p}{{\bf{A}}}_{i}{L}^{i}$ (I_M is the identity matrix of size M) is a vector autoregressive (VAR) polynomial of order p defined by the M × M coefficient matrices A₁,...,A_p , and

$$\begin{eqnarray*}\mathrm{diag}\left({{\rm{\nabla }}}^{{\bf{d}}}\right)=\left[\begin{array}{cccc}{\left(1-L\right)}^{{d}_{1}} & 0 & \ldots & 0\\ 0 & {\left(1-L\right)}^{{d}_{2}} & \ldots & 0\\ \vdots & \vdots & & \vdots \\ 0 & 0 & \ldots & {\left(1-L\right)}^{{d}_{M}}\end{array}\right],\end{eqnarray*}$$

and ${\left(1-L\right)}^{{d}_{i}},i\,=\,1,\ldots ,M$, is the fractional differencing operator defined by:

$$\begin{eqnarray}&&{\left(1-L\right)}^{{d}_{i}}=\sum _{k=0}^{\infty }{G}_{k}^{(i)}{L}^{k},\quad {G}_{k}^{(i)}=\displaystyle \frac{{\rm{\Gamma }}\left(k-{d}_{i}\right)}{{\rm{\Gamma }}\left(-{d}_{i}\right){\rm{\Gamma }}(k+1)},\end{eqnarray} \tag{ 9 }$$

with Γ( · ) denoting the Gamma (generalized factorial) function. The VARFI model is stationary when all the roots of $\det [{\boldsymbol{A}}(L)]$ are outside the unit circle and − 0.5 < d_i < 0.5 for i = 1...M, while it is nonstationary but mean reverting for 0.5 ≤ d_i < 1 (Baillie 1996, Velasco 1999). The coefficients of the polynomial A(L) allow the description of the short term dynamics, while the parameter ${\bf{d}}=\left({d}_{1},\ldots ,{d}_{M}\right)$ in equation (8) determines the long-term behavior of each individual process.

The parameters of the VARFI(p, d) model (8), namely the coefficients of A(L) and the variance of the innovations ${{\rm{\Sigma }}}_{{\bf{E}}}={\mathbb{E}}\left[{{\bf{E}}}_{n}^{T}{{\bf{E}}}_{n}\right]$, are generally obtained from process realizations of finite length first estimating the differencing parameters d_i using the Whittle semi-parametric local estimator (Beran et al 2016) individually for each process X_i ; then defining the filtered data ${X}_{i,n}^{(f)}={\left(1-L\right)}^{{d}_{i}}{X}_{i,n};$ and finally estimating the VAR parameters from the filtered data ${{\bf{X}}}_{n}^{(f)}$ using the ordinary least squares method to solve the VAR model ${\bf{A}}(L){{\bf{X}}}_{n}^{(f)}={{\bf{E}}}_{n}$, with model order p assessed through the Bayesian information criterion (Faes et al 2012).

Here, we approximate the VARFI process (8) with a finite order VAR process by truncating the fractional integration at a finite lag q, as follows

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{diag}\left({{\rm{\nabla }}}^{{\rm{d}}}\right)\approx {\bf{G}}(L)=\mathrm{diag}\left[\displaystyle \sum _{k=0}^{q}{G}_{k}^{(1)}{L}^{k}\quad \ldots \quad \displaystyle \sum _{k=0}^{q}{G}_{k}^{(M)}{L}^{k}\right]\\ \quad =\displaystyle \sum _{k=0}^{q}\mathrm{diag}\left[{G}_{k}^{(1)},\ldots ,{G}_{k}^{(M)}\right]{L}^{k}=\displaystyle \sum _{k=0}^{q}{{\bf{G}}}_{k}{L}^{k}.\end{array}\end{eqnarray} \tag{ 10 }$$

This allows us to express the VARFI(p, d) process as a VAR(m) process, with m = p + q

$$\begin{eqnarray}&&{\bf{B}}(L){{\bf{X}}}_{n}={{\bf{E}}}_{n},\end{eqnarray} \tag{ 11 }$$

with the coefficients in B(L) given by

$$\begin{eqnarray}&&{\bf{B}}(L)={\bf{A}}(L){\bf{G}}(L)=\left({{\bf{I}}}_{M}-\sum _{i=1}^{p}{{\bf{A}}}_{i}{L}^{i}\right)\left(\sum _{k=0}^{q}{{\bf{G}}}_{k}{L}^{k}\right)={{\bf{I}}}_{M}-\sum _{k=0}^{p+q}{{\bf{B}}}_{k}{L}^{k},\end{eqnarray} \tag{ 12 }$$

yielding, for q ≥ p,

$$\begin{eqnarray}&&\begin{array}{l}{B}_{0}={I}_{M},\\ {B}_{k}=\left\{\begin{array}{l}-{G}_{k}+\displaystyle \sum _{i=1}^{k}{A}_{i}{G}_{k-i},\quad k=1,\ldots ,p\\ -{G}_{k}+\displaystyle \sum _{i=1}^{p}{A}_{i}{G}_{k-i},\quad k=p+1,\ldots ,q\\ \displaystyle \sum _{i=0}^{p+q-k}{A}_{i+k-q}{G}_{q-i},\quad k=q+1,\ldots ,q+p\end{array}\right..\end{array}\end{eqnarray} \tag{ 13 }$$

In this section, we describe how to compute the information measures defined in section 2.1 across multiple temporal scales, under the hypothesis that the analyzed multivariate process is appropriately modeled by the VARFI representation provided in section 2.2. The procedure for multiscale analysis extends the rescaling approach proposed in (Faes et al 2017). Here only the fundamental steps are presented, the mathematical details are provided in (Faes et al 2017) and in the Appendixes of (Martins et al 2020).

Typically, to represent a scalar stochastic process at the temporal scale defined by the scale factor τ, a two-step procedure is employed which consists in first filtering the process with a low pass filter with cutoff frequency f_τ = 1/(2τ), and then downsampling the filtered process using a decimation factor τ (taking one every τ samples) (Porta et al 2006, Faes et al 2017). Extending this approach to the multivariate case, we first implement the following linear finite impulse response (FIR) filter

$$\begin{eqnarray}&&{{\bf{X}}}_{n}^{(r)}={\bf{D}}(L){{\bf{X}}}_{n},\end{eqnarray} \tag{ 14 }$$

where r denotes the filter order ${\bf{D}}(L)={\sum }_{k\,=\,0}^{r}{{\bf{I}}}_{M}{D}_{k}{L}^{k}$, and the coefficients of the polynomial D_k , k = 1,...,r, are the same for all scalar processes X_j ∈ X and are chosen to set up a low pass FIR configuration with cutoff frequency 1/(2τ).This step transforms the VAR(p + q) process (11) into a VARMA(p + q, r) process with moving average (MA) part determined by the FIR filter coefficients

$$\begin{eqnarray}&&{\bf{B}}(L){{\bf{X}}}_{n}^{(r)}={\bf{D}}(L){\bf{B}}(L){{\bf{X}}}_{n}={\bf{D}}(L){{\bf{E}}}_{n}.\end{eqnarray} \tag{ 15 }$$

Then, we exploit the connection between VARMA processes and state space (SS) processes (Aoki and Havenner 1991) to evidence that the VARMA process (10) can be expressed in SS form as

$$\begin{eqnarray}&&{{\bf{Z}}}_{n+1}^{(r)}={{\bf{B}}}^{(r)}{{\bf{Z}}}_{n}^{(r)}+{{\bf{K}}}^{(r)}{{\bf{E}}}_{n}^{(r)},\end{eqnarray} \tag{ 16a }$$

$$\begin{eqnarray}&&{{\bf{X}}}_{n}^{(r)}={{\bf{C}}}^{(r)}{{\bf{Z}}}_{n}^{(r)}+{{\bf{E}}}_{n}^{(r)},\end{eqnarray} \tag{ 16b }$$

where ${{\bf{Z}}}_{n}^{(r)}={\left[{{\bf{X}}}_{n-1}^{(r)}\cdots {{\bf{X}}}_{n-m}^{(r)}{{\bf{E}}}_{n-1}\cdots {{\bf{E}}}_{n-r}\right]}^{T}$ is a (m + r) - dimensional state process, ${{\bf{E}}}_{n}^{(r)}={{\bf{D}}}_{0}{{\bf{E}}}_{n}$ is the SS innovation process, and the vectors K^(r) and C^(r) and the matrix B^(r) can be obtained from B(L) and D(L). Further details can be found in appendix B of (Martins et al 2020).

In the second step of the rescaling procedure, the filtered process is downsampled in order to complete the multiscale representation. This is achieved by applying the results in, (Barnett and Seth 2015, Solo 2016, Faes et al 2017) which allow describing the filtered SS process after downsampling in the form

$$\begin{eqnarray}&&{{\bf{Z}}}_{n+1}^{(\tau )}={{\bf{B}}}^{(\tau )}{{\bf{Z}}}_{n}^{(\tau )}+{{\bf{K}}}^{(\tau )}{{\bf{E}}}_{n}^{(\tau )},\end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray}&&{{\bf{X}}}_{n}^{(\tau )}={{\bf{C}}}^{(\tau )}{{\bf{Z}}}_{n}^{(\tau )}+{{\bf{E}}}_{n}^{(\tau )}.\end{eqnarray} \tag{ 17b }$$

equations (17) provide the SS form of the filtered and downsampled version of the original VARFI(p, d) process, and parameters $\left({{\bf{B}}}^{(\tau )},{{\bf{C}}}^{(\tau )},{{\bf{K}}}^{(\tau )},{{\rm{\Sigma }}}_{{{\bf{E}}}^{(\tau )}}\right)$ can be obtained from the SS parameters before downsampling and from the downsampling factor τ.

In this section, we show how to compute analytically the information decomposition of a jointly Gaussian multivariate stochastic process starting from its associated SS model (17).

The derivations are based on the knowledge that the linear parametric representation of Gaussian processes captures all the entropy differences that define the various information measures (Barrett et al 2010). These entropy differences are related to the partial variances of the present of the target conditioned to its past and the past of one or more sources. The partial variances can be formulated as variances of the prediction errors resulting from linear regression (Faes et al 2015, 2017). Specifically, let us denote as ${E}_{j| j,n}={X}_{j,n}-{\mathbb{E}}\left[{X}_{j,n}| {X}_{j,n}^{-}\right]$ and ${E}_{j| {ij},n}={X}_{j,n}-{\mathbb{E}}\left[{X}_{j,n}| {X}_{i,n}^{-},{X}_{j,n}^{-}\right]$ the prediction errors of a linear regression of X_j,n on ${X}_{j,n}^{-}$ and $\left({X}_{j,n}^{-},{X}_{i,n}^{-}\right),$ respectively. Then, the TE from X_i to X_j can be expressed as

$$\begin{eqnarray}&&{T}_{i\to j}=\displaystyle \frac{1}{2}\mathrm{ln}\displaystyle \frac{{{\boldsymbol{\Sigma }}}_{{E}_{j| j}}}{{{\boldsymbol{\Sigma }}}_{{E}_{j| {ij}}}}.\end{eqnarray} \tag{ 18 }$$

Similarly, the joint TE from $\left({X}_{i},{X}_{k}\right)$ to X_j can be defined as

$$\begin{eqnarray}&&{T}_{{ik}\to j}=\displaystyle \frac{1}{2}\mathrm{ln}\displaystyle \frac{{{\boldsymbol{\Sigma }}}_{{E}_{j| j}}}{{{\boldsymbol{\Sigma }}}_{{E}_{j}}},\end{eqnarray} \tag{ 19 }$$

where ${{\boldsymbol{\Sigma }}}_{{E}_{j}}={\mathbb{E}}\left[{E}_{j,n}^{2}\right]$ is the variance of the prediction error of a linear regression of X_j,n on ${{\bf{X}}}_{n}^{-},{E}_{j,n}={X}_{j,n}\,-{\mathbb{E}}\left[{X}_{j,n}| {{\bf{X}}}_{n}^{-}\right]$. Based on these derivations, one can easily complete the IID decomposition of TE by computing T_k→j as in (18) and deriving the interaction TE from (3) and the PID decomposition, as well by deriving the redundant TE from (7), the synergistic TE from (6) and the unique TEs from (5).

Finally, we show how to compute any partial variance from the parameters of an SS model in the form (17) at any assigned time scale τ (Barnett and Seth 2015, Solo 2016). The partial variance ${{\rm{\Sigma }}}_{{E}_{j| a}^{(\tau )}}$, where the subscript a denotes any combination of indexes ∈ {1,...,M}, 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 (17a ) and observation equation

$$\begin{eqnarray}&&{{\bf{X}}}_{a,n}^{(\tau )}={{\bf{C}}}_{a}^{(\tau )}{{\bf{Z}}}_{n}^{(\tau )}+{{\bf{E}}}_{a,n}^{(\tau )},\end{eqnarray} \tag{ 20 }$$

where the additional subscript _a denotes the selection of the rows with indices a in a vector or a matrix. These submodels can be converted to the SS form as in (17), with innovation covariance ${{\rm{\Sigma }}}_{{{\bf{E}}}_{a}^{(\tau )}}$, so that the partial variance ${{\rm{\Sigma }}}_{{E}_{j| a}^{(\tau )}}$ is derived as the diagonal element of ${{\rm{\Sigma }}}_{{{\bf{E}}}_{a}^{(\tau )}}$ corresponding to the position of the target X_j,n (Faes et al 2017, Martins et al 2020).

The cardiovascular and respiratory time series were measured in a group of 61 healthy subjects (19.53.3 years old, 24 males) 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 arithmetic test (MA). The head-up tilt protocol consisted in tilting the motor-driven bed table where the subject was lying to 45 degrees in order to evoke mild orthostatic stress. During the MA task the subjects were instructed to sum up three-digit numbers displayed on the ceiling of the examination room by a data projector, and to decide whether the final one-digit number was odd or even pushing the button projected on the ceiling. During the whole protocol, the volunteers were asked to avoid disturbing movements and speaking (Javorka et al 2017). During all measurements, the subjects were free-breathing. The study was approved by Ethical Committee of the Jessenius Faculty of Medicine, Comenius University (Slovakia) and all participants signed a written informed consent. More details about the experimental protocol are reported in (Javorka et al 2017).

The acquired signals were the surface electrocardiogram (ECG, horizontal bipolar thoracic lead; CardioFax ECG-9620, NihonKohden, Japan), the finger arterial blood pressure (Finometer Pro, FMS, Netherlands) recorded noninvasively by the volume-clamp photoplethysmographic method, and the respiration signal recorded through respiratory inductive plethysmography (RespiTrace, NIMS, USA). All measured signals were digitised at 1000 Hz. For each subject and experimental condition, the values of H, S and R 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 three time series H, S and R were interpreted as realizations of the stochastic processes descriptive of the cardiac, vascular, and respiratory dynamics. For each subject and condition, the analyzed multivariate process is defined as ${\bf{X}}=\left[{X}_{{\rm{H}}},{X}_{{\rm{S}}},{X}_{{\rm{R}}}\right]$.

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 H and when possible, erroneous/missing intervals were substituted by pulse intervals measured as the difference in time between two consecutive S measurements $\left({{\rm{\Delta }}}_{{t}_{{\rm{S}}}}(n)={t}_{{\rm{S}}}(n+1)-{t}_{{\rm{S}}}(n)\right)$. The three time series were normalized to zero mean and unit variance before multiscale analysis.

To compute the IID and PID measures, the approach based on complete VARFI model identification defined in section 2.3 was applied. The VARFI model was identified first estimating the fractional differencing parameter d_i , i = 1,...,3, individually for each time series using the Whittle estimator, then filtering the time series with the fractional integration polynomial truncated at a lag q, and finally estimating the parameters of the polynomial relevant to the short-term dynamics via least squares VAR identification. Theoretically, the VARFI is of infinite order, hence the value of q has to be selected to approximate the VARFI process with a finite order VAR process. Several previous studies (Bardet et al 2003, Faes et al 2019) defined q = 50 as 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 TE measures and the smoothing of the non-monotonic trends with the time scale (Faes et al 2019). The order p of the VAR model was assessed individually for each subject and in each phase of the experimental protocol, by the Bayesian information criterion (BIC) (Martins et al 2020). Then, multiscale TE measures were computed implementing FIR low pass filter of order r = 48, for time scales τ in the range (1,...,12), which corresponds to low pass cutoff normalized frequencies f_τ = (0.5,...,0.04). The value r = 48 was set according to previous settings (Faes et al 2019).

The differencing parameters d_i were estimated individually for each time series in the interval [−0.5, 1]. For 3 individuals the estimated d_i parameters were near to 1 which indicate that the estimated VARFI models were nonstationary and thus only 59 subjects were considered for further statistical analysis.

Significant changes in the information transfer and modification measures across the pairs of experimental conditions SU₁ versus UP and SU₂ versus MA are evaluated via a linear mixed-effects model, incorporating both fixed and random effects (Pinheiro and Bates 2006). 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. Furthermore, the interaction between the factors is also considered. To assess the changes of interest, estimated marginal means (EMM) (Searle et al 1980) are obtained for each difference, UP − SU₁ and MA − SU₂, at each time scale, τ = 1,...,12. A Z-test is applied to check the significance of these differences at a level p < 0.05 with the Tukey correction for multiple comparisons.

Additionally, we used a similar approach to evaluate significant differences between measures, particularly, between the individuals TEs from S to H and from R to H, and between the redundant and synergistic TEs, at each time scale τ. This allows to ascertain which of the source processes, S or R, is prevalent in driving H at a given time scale, or whether redundant effects are prevalent over synergistic ones. The comparison between redundancy and synergy is also very important as it provides clues about the type of interactions between source processes, making it possible to corroborate some physiological assumptions. In this case, the fixed-effects (or factors) are measure and scale, and the random-effect is the subject-dependent intercept as in the previous case. For both models, residuals were checked for whiteness. The packages lme4 (Bates et al 2015) and emmeans (Lenth 2019) of the R software (R Core Team 2016) were used to build the models and to compute EMM, respectively.

The identification of the VAR model returned rather stable values of the model order: the median value was p = 4 for each phase of the protocol, while the interquartile interval was [4, 5] for SU₁, [3, 5] for UP, [3, 4] for SU₂, and [3, 4] for MA. The results of multiscale analysis performed for the IID terms as defined in equations (1)–(4) as well for the redundant and synergistic TEs of the PID equations (5) $\left({R}_{{ik}\to j}\right.$ and S_ik→j ) are provided next.

First, we analyze the multiscale interactions when the heart period time series was considered as the target, and the systolic pressure and respiration were considered as the sources. Figure 3 presents the median and quartiles across subjects of the six information measures computed as a function of the time scale τ = 1,...,30, for SU₁ versus UP, panel (a) and SU₂ versus MA, panel (b). Statistically significant changes (p < 0.05) in TE measures at each time scale across the pairs SU₁ versus UP or SU₂ versus MA are marked with ∗. From a visual inspection of the multiscale patterns one can infer a markedly higher T_S→H at lower scales up to τ ≈ 10 moving from SU₁ to UP, panel (a). The values of T_R→H are lower in UP for scale 1. In contrast, for scales 3—7 this measure is higher in UP phase. We can observe a similar behaviour in the T_S,R→H where the values of this measure are significantly higher from τ = 2 to τ = 9. The I_S,R→H decreases significantly with tilt, mainly in the first three time scales. Then, significant differences between SU₁ and UP can be observed in the mid-range time scales τ = 5, 6, 7. Observing R_S,R→H and S_S,R→H as two distinct information modification measures we see that both exhibit significant differences up to time scale 7 with higher values in the upright phase (UP). Note that in the first few time scales the redundancy is greater than the synergy, hence the negative values of ITE in these time scales.

Zoom In Zoom Out Reset image size

Moving from SU₂ to MA, we find significantly higher values of T_S→H in the mid-range of time scales. The T_R→H exhibits statistically significant differences (a decrease) only for the first scales. Similar behaviour is observed for the joint TE, with significant higher values in the MA phase for τ = 1, 4. Regarding I_S,R→H , R_S,R→H and S_S,R→H the model was not able to find significant differences. The profile of I_S,R→H , R_S,R→H and S_S,R→H in recovery supine position (SU₂) and in the mental stress phase (MA) are quite similar. For this reason, no significant difference were detected. Figure 4 presents the estimated marginal means for the differences between T_S→H − T_R→H and R_S,R→H − S_S,R→H and the respective 95% confidence intervals for the 4 positions of the experimental protocol. Statistically significant differences (p < 0.05) are marked with ∗.

Zoom In Zoom Out Reset image size

In the resting supine position (SU₁, figure 4(a)) we find prevalence of the redundancy over synergy at τ = 1. For the other time scales, no significant differences are observed for these two measures of information modification. Comparing T_S→H and T_R→H , the T_R→H overcomes the T_S→H in the first three time scales. This behavior is inverted at τ = 5, 6 where the value of T_S→H is superior. In the remaining time scales no significant differences were found.

Moving to the upright position (UP, figure 4(b)) a prevalence of the redundancy is again observed in the first time scales τ = 1, 2. For the other time scales no significant differences between R_S,R→H and S_S,R→H were detected. Regarding the individuals TEs T_S→H and T_R→H we observed a similar behaviour to that observed in the SU₁, that is, the oscillations observed in the heart period are essentially of respiratory origin. For the mid-range scales τ = 5, 6, 7, 8, an inversion is observed, where oscillations of vascular origin prevail over those of respiratory origin. No significant differences were found for the remaining time scales.

The results of the estimated EMM for the differences in recovery supine position (SU₂, figure 4(c)) are similar to those observed previously. We find only one significant positive difference for R_S,R→H − S_S,R→H at τ = 1 denoting a preponderance of redundancy. Regarding T_S→H − T_R→H , we note prevalence of the T_R→H at τ = 1, 2. On the other hand, at τ = 4, 5, 6 the individual TE T_S→H plays a more dominant role.

Finally, in the mental stress phase (MA, figure 4(d)) a prevalence of the redundancy is observed at the first time scale and no significant difference R_S,R→H − S_S,R→H was found in the remaining scales. The T_S→H − T_R→H are significant up to time scale τ = 9, however in the first two temporal scales τ = 1, 2 the oscillations observed in the heart period are predominantly of respiratory origin. In the remaining time scales, i.e., τ = 3...8, oscillations of vasomotor origin have a leading role. Figure 5 presents the median and quartiles across subjects of the six information measures computed as a function of the time scale τ, for the four conditions analyzed when S is considered the target process. In this case, moving from SU₁ to UP, a markedly lower value of the TE from heart period to systolic pressure, T_H→S , is observed for scales up to τ = 3, while the trend is inverted for scales from τ = 4 to τ = 9 where the TE increased significantly. Both the TE from respiration to systolic pressure and the joint TE directed to systolic pressure display significantly lower values across a wide range of scales (τ = 1,...,5 for T_R→S , and τ = 2,...,8 for T_H,R→S ). In the first two time scales there is a significant increase of the ITE (figure 5(a4)) going from SU₁ to UP. This trend is also observed in the S_H,R→S profile (figure 5(a6)). Interestingly, the ITE I_H,R→S and the redundant transfer R_H,R→S exhibit a different response to tilt at different time scales: when τ = 1 the ITE increased as a consequence of the reduction in the redundancy, while for τ ≥ 3 the ITE decreased as a consequence of the increased redundancy (figures 5(a4), (a5)).

Zoom In Zoom Out Reset image size

A very different behaviour of the information measures is reported moving from SU₂ to MA (figure 5(b)). The individual and joint information transfer to S displayed no statistically significant variations moving from relaxation to mental stress (figure 5(b1), (b2) and (b3)). In the first two time scales the R_H,R→S is higher than the S_H,R→S (figure 5(b5) and (b6)), determining negative values of ITE (figure 5(b4)). The mental stress condition induces a decrease of redundancy and has no influences on slower time scales (τ > 2). Figure 6 presents the median and quartiles across subjects of the six information measures computed as a function of the time scale τ, for the four condition analised when R is considered the target process. The postural stress determined a decrease of the information transfer from H to R at small time scales (τ = 1, 2, figure 6(a1)), and an increase of the information transfer from S to R at intermediate time scales (τ = 2, 3, 4, figure 6(a2)). The joint information transfer from H and S to R decreased moving from SU₁ to UP for τ = 1, and increased markedly for τ = 2, 3, 4 (figure 6(a3)). The ITE directed to R increased during the postural stress for τ = 1, 2, 3; this trend is the result of a significant increase of synergy at small time scales (τ = 1 − 4, figure 6(a6)), partly compensated for τ = 1 by an increase of redundancy (figure 6(a5)).

Zoom In Zoom Out Reset image size

Again, the behaviour of the information measures is very different during mental stress compared with postural stress (figure 6(b)). As observed before when S was the target, also with target R only a few significant changes were detected moving from SU₂ to MA, consisting in increased values of T_S→R , T_H,S→R , R_H,S→R , and S_H,S→R at the faster time scales.

The increase of transfer entropy from S to H at lower time scales (figure 3 a1) when going from rest to tilt is in agreement with previous studies (Westerhof et al 2006, Faes et al 2011a, Krohova et al 2019) indicating indicating the dominance of baroreflex-mediated interactions during orthostatic stress. The statistically significant differences evidenced for mid-range time scales also indicate that S has the most relevant part of its dynamics in the LF and VLF bands, as already noted in (Krohova et al 2019), and that postural stress induces changes to such dynamics noticeable in a wider range of time scales. For the UP position, at τ = 1 the information transfer from R to H (figure 3 a2) is significantly lower than in rest. This is in accordance with previous findings evidencing that the respiratory sinus arrhythmia (RSA) decreases with tilt (Porta et al 2012, Javorka et al 2018, Krohova et al 2019), probably due to the decreased parasympathetically mediated RSA. The trend is different for longer time scales, at which the information transfer from R to H becomes significantly higher than in rest, indicating the prevalence of slower oscillations in the information transfer from R to H, i.e. slowly varying respiration influences (mostly related to spontaneous changes of the respiratory pattern) are transferred more to slower HRV oscillations during postural stress (Krohova et al 2019). The two previous effects produce a higher joint information transfer (figure 3 a3) during UP for time scales from 2 to 9, denoting stronger redundancy, as seen in figure 3 a4 from the negative I and from figure 3 a5-a6 (higher redundancy than synergy) for time scales up to 7. Physiologically, such results confirm previous findings indicating that postural stress produces a strong activation of baroreflex-mediated RSA, path R → S → H, especially for time scales longer than τ > 1. This confirms that baroreflex is a fundamental mechanism for slower heart rate oscillations, which should be studied in the LF band (Krohova et al 2019, Pernice et al 2021).

Regarding the comparison between mental stress and SU₂, the transfer entropy T_S→H is higher during MA if compared to rest at mid-range time scales (3 ≤ τ ≤ 6) (figure 3 b1), reflecting the activation of vasomotion associated with enhancement of slower blood pressure oscillations. Conversely, there is a significant decrease of information transfer from R to H at lower time scales (τ = 1, 2) during mental stress (figure 3 b2), indicating an overall weakening of the influence of respiration on heart rate, which is in agreement with vagal inhibition, but also the lower involvement of baroreflex-mediated RSA, provoked by stress challenges already demonstrated in previous works (Faes et al 2011a, Javorka et al 2018, Krohova et al 2019), and it is also in accordance with reduced cardiorespiratory interactions and synchronization occurring during mental task (Zhang et al 2010, Widjaja et al 2013, Pernice et al 2021). The combination of these results suggests that mental stress not only produces the weakening of RSA due to vagal inhibition, but also the non-activation of baroreflex-mediated RSA (R → S → H), differently from postural stress; these observations are supported by previous findings in the literature (Krohova et al 2019). The results of joint information transfer (figure 3 b3) evidence the progressive increase occurring during the mental stress from time scale 1 to time scale 4: for τ = 1 there is a significant decrease during MA driven by the prevalence of respiratory dynamics, then there is the activation of baroreflex effects for mid-range time scales that produces the prevalence of information transfer from S dynamics with regard to slower oscillations (the increase becomes statistically significant in MA just for τ = 4). Such effect is not statistically significant for τ = 5 and higher scales, differently from postural stress (figure 3 a3), given the lower influence of respiratory dynamics already starting from τ = 3. This may be due to the the combined effect of (i) vasomotor reactions elicited by MA through cortical mechanisms reflected in SAP changes and then transferred to HRV through the baroreflex and (ii) the reduced RSA due to the significant increase in the breathing rate during MA if compared to SU₂; again, such interpretations agree with (Krohova et al 2019). This difference between the stress conditions and the corresponding resting state across the various time scales emphasizes the importance of employing a multiscale approach when studying cardiovascular and cardiorespiratory interactions, suggesting also that such interactions include complex multiscale patterns which respond flexibly to stress challenges (Widjaja et al 2013, Krohova et al 2019). No statistically significant differences between MA and SU₂ are detected, at any time scale, analyzing the interaction transfer entropy (Figure 3 b4), and neither with regard to redundancy (Figure 3 b5) or synergy (figure 3 b6), conversely to postural stress, evidencing that the mechanisms underlying postural and mental stress are different, not only when 'raw' signals (i.e. at τ = 1) are considered, but also when going through a full multiscale analysis. This also supports the importance of employing a multiscale approach to shed more light on such mechanisms, and its potential usefulness to differentiate between stress conditions.

The analysis of marginal means (figure 4) indicates the prevalence of redundancy only for τ = 1 during SU₁, and also for τ = 2 during tilt, and this can be put in relation to the activation of baroreflex due to the postural stress. Moreover, the difference between T_S→H and T_R→H indicates for τ < 4 the prevalence of information transfer from respiration to heart rate, evidencing the dominance of RSA on short time scales. The opposite is instead observed for mid-range scales (τ from 5 to 9), thus confirming once again that S has the most relevant part of its dynamics in the LF and VLF bands, and that transfers more information to H than respiration when assessed only for slower oscillations (Krohova et al 2019). Similar trends are reported comparing the two different stress typologies (MA and UP), with the only difference being evidenced for τ = 3, with a prevalence of redundancy for S at τ = 2 during postural stress. Conversely, we notice an increase of the information transfer from S (instead than R) to H during mental stress,ggesting that mental stress strengthens respiration-unrelated baroreflex effects. Overall, the obtained results highlight that head-up tilt induces scale-dependent variations in the transfer entropy of arterial pressure, higher in the mid-scales associated with slow oscillations, and lower associated to the effects of respiration. This result is similar to what observed in (Faes et al 2019) with regard to complexity of arterial pressure time series.

It is well known that the pathway from H to S is influenced by the Frank-Starling law according to which heart period variations affect the end diastolic volume and thus the strength of the systolic contraction. The decrease of TE from H to S during postural stress (figure 5(a1)) is in agreement with the decrease of the causal interactions from H to diastolic arterial pressure (DBP) and from DBP to S reported in (Westerhof et al 2009, Javorka et al 2017). Physiologically, this may be due to the changes in the peripheral vascular resistance associated with sympathetic activation during baroreceptors unloading, or alternatively to the a lower magnitude of heart rate oscillations associated with orthostasis (Cooke et al 1999). Our results expand such findings evidencing that these trends are then inverted for longer time scales, for which higher TE from H to S is observed; such increased transfer is associated with a tilt-induced shift from increasing ITE at short time scales to decreasing ITE at longer scales (figures 5(a1), (a4)). Thus, our multiscale analysis reveals that the postural stress decreases the redundancy between H and R while driving S at short scales, and increases such redundancy at longer scales; physiologically, this may be an indication of complex interactions between the slower rhythms of heart rate variability and the modulation of respiratory amplitude at these lower frequencies. The physiological interpretations made for postural stress, considering S as the target process, are not valid for mental stress, similar to what found in (Javorka et al 2017), where however the strength of the connection along the the direction from DBP to S was also decreased during mental stress.

When respiration was taken as the target of multiscale analysis of information transfer, we documented that postural stress induces at short temporal scales lower values of the information transfer from heart period together with higher values of the information transfer from systolic pressure and of the joint information (figure 6(a1)–(3)). These results are associated with increased synergy between H and S while they transfer information to R (figure 6(a1)–(3)). Previous studies have suggested that the heart rate dynamics are not only driven by respiration, but also influence it, resulting in bidirectional coupling between H and R that reflect a closed-loop cardiorespiratory interaction (Saul et al 1989, Yana et al 1993, Perrott and Cohen 1996, Porta et al 2013). Our results approach the study of these interactions from an entirely new point of view, and document the synergistic role played by heart rate and arterial pressure slow oscillations in influencing the respiratory dynamics during head-up tilt. The enhancement of synergy suggests that postural stress evokes separate mechanisms whereby H and S drive R during tilt, the first possibly related to afferent commands to the respiratory centers in the central autonomic network and the second to mechanical influences (Benarroch 1993).

On the other hand, during the mental stress, the low values of T_S→R and its increase for τ = 1, 2, suggests that S influences R only for faster time scales during MA. The lower values of T_S→R compared to T_R→S are in agreement with a previous study (Porta et al 2013) that underlined that systolic pressure variations at the respiratory rate are due to the respiratory-related fluctuations of intrathoracic pressure modulating right and left preloads and, in turn, stroke volume, while the negligible causality in the other direction is due to the fact that fast neural actions do not influence the relation between R and S.

Our simulation study revealed that long-range correlations have a substantial impact on the multiscale information transfer from one process to another, as well as on the synergistic/redundant interactions between two source processes sending information to the target. In particular, we noticed that increasing the strength of long-range correlations in a source process reduces the information transferred towards the target and enhances synergistic effects over redundant ones (figure 2(a)), while the opposite trends (i.e., higher information transfer and higher redundancy) are observed when long-range correlations become stronger in the target process (figure 2(b)). These findings can be translated to the results of our experimental analysis, also considering previous results showing that long-range correlations are almost absent in respiration variability while they are substantial in heart period and systolic arterial pressure variability and tend to increase with physiological stress (Martins et al 2020). In fact, the presence of long-range correlations in the target process H, and the strengthening during postural stress, can explain the increased information transfer and the increased redundancy observed during tilt in the analysis of cardiovascular and cardiorespiratory multivariate interactions (figure 5(a)). In the same way, long-range correlations are very important during tilt for the systolic pressure (Martins et al 2020), and this may explain the increased information transfer and the increased redundant transfer toward S at intermediate and long time scales observed during postural stress (figure 6(a)). Interestingly, opposite trends are observed when the respiration signal is taken at the target of the analysis of information transfer: the postural stress induces lower information transfer from H to R and higher synergistic transfer from H and S to R (figure 6(a)); in this case, long-range correlations are in the source process H and this explains the observed results in agreement with the simulation of figure 2(a).