Lag-specific transfer entropy as a tool to assess cardiovascular and cardiorespiratory information transfer

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication.

1

Lag-specific Transfer Entropy as a Tool to Assess Cardiovascular and Cardiorespiratory Information Transfer Luca Faes*, Member, IEEE, Daniele Marinazzo, Alessandro Montalto. and Giandomenico Nollo, Member, IEEE  Abstract— In the study of interacting physiological systems, model-free tools for time series analysis are fundamental to provide a proper description of how the coupling among systems arises from the multiple involved regulatory mechanisms. This study presents an approach which evaluates direction, magnitude and exact timing of the information transfer between two time series belonging to a multivariate data set. The approach performs a decomposition of the well known transfer entropy (TE) which achieves (i) identifying, according to a lag-specific information-theoretic formulation of the concept of Granger causality, the set of time lags associated with significant information transfer, and (ii) assigning to these delays an amount of information transfer such that the total contribution yields the aggregate TE. The approach is first validated on realizations of simulated linear and nonlinear multivariate processes interacting at different time lags and with different strength, reporting a high accuracy in the detection of imposed delays, and showing that the estimated lag-specific TE follows the imposed coupling strength. The subsequent application to heart period, systolic arterial pressure and respiration variability series measured from healthy subjects during a tilt test protocol illustrated how the proposed approach quantifies the modifications in the involvement and latency of important mechanisms of short-term physiological regulation, like the baroreflex and the respiratory sinus arrhythmia, induced by the orthostatic stress. Index Terms— autonomic nervous system, cardiovascular control, conditional entropy, dynamical systems, Granger causality, multivariate time series, mutual information.

I. INTRODUCTION

T

HE SHORT-TERM cardiovascular and cardiorespiratory regulation, which is manifested in the spontaneous beat-

Copyright (c) 2013 IEEE. Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending an email to [email protected] This work was supported by the Healthcare Research Implementation Program (IRCS), Provincia Autonoma di Trento and Bruno Kessler Foundation, Italy, by the University of Gent (Special Research Funds for visiting researchers) and by the Belgian Science Policy (IUAP VII project CEREBNET P7 11). L. Faes and G. Nollo are with the IRCS-FBK and BIOtech, Dept. of Industrial Engineering, University of Trento, Italy. * Corresponding author address: BIOtech, Università di Trento, via delle Regole 101, 38123 Mattarello, Trento, Italy (e-mail: [email protected]). D. Marinazzo and A. Montalto are with the Department of Data Analysis, University of Gent, Gent, Belgium.

to-beat variability of heart period (HP, measured as the inverse of the heart rate), arterial pressure (AP) and respiratory activity (RA), is known to be the result of the dynamical interaction among several physiological mechanisms, including both reflex pathways and central and peripheral oscillators [1]. In particular, the cardiovascular regulation is accomplished through a closed-loop interaction between HP and systolic AP (SAP), whereby the cardiac baroreflex acts as a delayed negative feedback of SAP on HP, and modifications of ventricular filling and diastolic runoff determine feedforward effects from HP to SAP [1-4]. Moreover, cardiac and vascular dynamics are continuously perturbed by RA, e.g. through the mechanical effects of breathing on intra-thoracic pressure and stroke volume [5] and through the coupling between respiratory centers and vagal outflow [6]. To address this complex picture of physiological interactions several methods for multivariate time series analysis, ranging from parametric models studied in time or frequency domains to model-free approaches working in the information domain [7], have been devised and applied in recent years to the short-term variability of HP, AP and RA [8]. While the usefulness of model-based approaches is indubitable, model-free techniques are gaining more and more attention due to the awareness that, especially in integrated system physiology, the likelihood that a predefined model class matches with good approximation the mechanism generating the observed dynamics can be very low. A natural framework for the model-free analysis of coupled systems is provided by information theory, where a well known measure of directional interaction between two time series of a multivariate data set is the so-called transfer entropy (TE) [9]. The popularity of TE stems from its inherent ability to incorporate directional and dynamical information, its sensitivity to both linear and nonlinear interactions, and its close connection with the ubiquitous concept of Granger causality [10]. Additionally, the recent introduction of dataefficient estimation procedures aimed at limiting the number of variables involved in entropy computation has made it possible to provide reliable TE estimates from short realizations of multivariate time series [11-15]. A limitation of the TE is that it is not lag-specific, i.e., it quantifies the information transfer between time series without detecting the timing through which such a transfer occurs.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication.

However, evaluating the exact coupling delays among interacting time series is of great importance to gain a complete understanding of the function of complex networks [16]. Time delays arise naturally in coupled systems as a consequence of propagation effects, or may be introduced through closed-loop control schemes [17]. This is the case for the networks sub-serving cardiovascular and cardiorespiratory regulation, where assessing time-lagged information transfers would be desirable to complement the evaluation of direction and strength of a detected interaction with the evaluation of its timing, thus providing useful information about the latency of the physiological mechanisms under investigation. The present study introduces an approach for the detection of time-lagged information transfer between two time series belonging to a multivariate data set. The approach decomposes the TE between the series into a sum of terms, the so-called lag-specific TEs, each quantifying the information transfer at a given time delay. To implement the method, a procedure facing the problem of entropy estimation in high dimensions is devised, which also provides the basis for assessing the lags associated with significant transfer of information. The overall approach is first validated on simulations reproducing linear and nonlinear interactions among multiple time series. Then, it is applied to HP, SAP and RA time series measured from healthy subjects during an experimental protocol inducing expected modifications of the physiological mechanisms underlying cardiovascular and cardiorespiratory regulation, i.e., head-up tilt testing. Part of this paper has been presented in a preliminary form in a conference contribution [18].

information (CMI), can be interpreted as the reduction of uncertainty about Yn when learning the past of X, if the past of Y and the past of Z are already known. The TE measures aggregate Granger-causal influence of X at past lags, and thus is not lag-specific. In order to characterize Granger-causal influences between processes for specific time lags, the definition of G-causality can be intuitively itemized as follows [13]: G-causality exists from X to Y at lag u, Xn-u→Yn, if Xn-u contains information that helps predicting Yn above and beyond the information contained in ─ ─ ─ Yn , Zn and in Xn \ Xn-u (where \ denotes subtraction from a set). In the information domain, this definition can be formulated stating that Xn-u→Yn if and only if IX→Y|Z(u)>0 [13], where the CMI IX→Y|Z(u) is defined as ─

─

─

IX→Y|Z(u)=I(Yn ; Xn-u | Yn , Zn , Xn \ Xn-u).

(2)

The CMI in (2) can be interpreted as the reduction of uncertainty about Yn when learning the past of X at lag u, if the past of X at all other lags and the whole past of Y and Z is already known. While the condition IX→Y|Z(u)>0 allows assessing lag-specific G-causality, the modulus of IX→Y|Z(u) cannot be directly related to the information transfer measured by the TE. To relate lag-specific and overall information transfers, we propose to decompose the TE as

TEX Y | Z 

k 1TEX Y | Z (uk ) L

(3) ─

II. METHODS A. Lag-Specific Granger Causality and Transfer Entropy Let us consider an overall dynamical system composed of M interacting (sub)systems, and assume that the states visited by the systems are described as discrete-time stationary stochastic processes. Suppose that we are interested in evaluating the information flow from an assigned source system, described by the process X, to a destination system, described by the process Y, in the presence of the remaining systems, described by the set of processes Z={Z(k)}k=1,...,M-2. Let us further denote Xn, Yn and Zn as the random variables obtained sampling the ─ ─ processes at time n, and Xn ={Xn-1,Xn-2,...}, Yn ={Yn-1,Yn-2,...} ─ and Zn ={Zn-1,Zn-2,...} as the sets of variables describing the past of the processes. Then, according to the definition of ─ Granger causality [10], G-causality from X to Y, Xn →Yn, ─ exists if Xn contains information that helps predicting Yn ─ ─ above and beyond the information contained in Yn and Zn . A more general formulation of this notion is based on conditional probabilities, and can be formulated in the information domain in terms of the well known TE [9] stating ─ that Xn →Yn if and only if TEX→Y|Z>0 [19], where ─

─

─

TEX→Y|Z = I(Yn ; Xn | Yn , Zn ).

(1)

The TE, which is formulated in (1) as a conditional mutual

2

─

where TEX→Y|Z (uk) = I(Xn-uk ; Yn | Xn-uk+1, ..., Xn-uL, Yn , Zn ), and the sum is extended to the set of lags ={u1,…,uL} with significant lag-specific G-causality (i.e., Xn-uk→Yn if and only if uk; here the lags are in ascending order). Thus, lagspecific TE analysis is performed first computing the CMI in (2) to determine the lags uk associated with significant Gcausality, and then computing the measures TEX→Y|Z (uk) to be used in (3) as the terms of the TE decomposition. As a result, this decomposition puts in evidence lag-specific information transfers, defined in a way such that their aggregate contribution yields the overall TE. The formal proof of the decomposition (3) is reported in the Appendix. B. Estimation Approach In this Section we describe a procedure for estimating the set  of the lags bearing G-causality according to (2), and then the overall TEX→Y|Z as well as the individual contributions TEX→Y|Z(u) according to (3). The procedure is based on recognizing that any CMI can be expressed as the difference between two conditional entropy (CE) terms; e.g., the TE can ─ ─ ─ ─ ─ be expressed as TEX→Y|Z=I(Yn ; Xn | Yn , Z n ) = H(Yn | Y n , Z n ) ─ ─ ─ − H(Yn | X n , Y n , Z n ), where H(Yn|V) is the CE measuring the entropy of the scalar variable Yn conditioned to the vector variable V. Then, we compute the CE according to a sequential procedure for non-uniform conditioning whereby the conditioning vector V is updated progressively taking all relevant processes into consideration at each step and selecting

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication.

3

After termination of the procedure for component selection, the conditioning vector is composed as V=[VX, VY, VZ], where VX, VY, and VZ denote the components of V belonging to X, Y, and Z. Then, the CMI in (2) can be estimated as:

I X Y | Z u   H (Yn V \ X n u )  H (Yn V ) ,

(4)

and, according to the definition of lag-specific G-causality, Xn-u→Yn if Xn-u is selected by the conditioning procedure, since IX→Y|Z(u)>0 if Xn-uV and IX→Y|Z(u)=0 if Xn-uV. Thus, the estimate of the set of lags bearing G-causality from X to Y is = { uk : Xn-uk  VX }. As a consequence, the lag-specific information transfer in (3) is estimated as TEX→Y|Z(u)=0 if u (i.e. if Xn-uV), and as Fig. 1. Example of lag-specific TE analysis performed for a realization of Simulation 1 with parameters 1=1, 2=2, 3=8, and C1=C2=C3=0.5. (a) Power spectral densities of the three processes generated according to (6). (b) left: CMI between the present component of the target process, Yn, and the component selected at the kth step of the conditioning procedure, Wk, conditioned to the vector of previously selected components Vk-1 (circles), plotted together with its corresponding threshold for significance (crosses); right: Entropy of Yn conditioned to the vector of components selected up to the kth step, Vk. The component selected at each step is indicated in the plot (the component in gray is discarded due to non-significant CMI). (c) Lagspecific TE (TEX→Y|Z(u), circles) and aggregate TE (TEX→Y|Z, dashed line); gray and white circles denote respectively uncoupled and coupled lags.

the components that mostly reduce the uncertainty about the target variable Yn [11,12,14]. Specifically, a set of initial candidate components is first defined including the past of all relevant processes up to a maximum lag Lmax, i.e., Ω={Xn-1,...,Xn-Lmax, Yn-1,...,Yn-Lmax, Zn-1,...,Zn-Lmax}. Then, the entropy of Yn conditioned to Ω is computed starting from an empty conditioning vector, V0=[·], and proceeding as follows: (i) at each step k1, form the candidate vector [W,Vk-1], where WΩ\Vk-1, and compute the CE H(Yn|W,Vk-1); repeat the computation for all possible candidates; (ii) retain the candidate component for which the estimated CE is minimum, i.e., set Vk=[Wk ,Vk-1] such that Wk=arg minW H(Yn|W,Vk-1); (iii) terminate the procedure when an irrelevant component has been selected, i.e. when I(Yn;Wk|Vk-1)=H(Yn|Vk-1)−H(Yn|Vk) is not significant. In this study, statistical significance of the reduction in the CE brought by the component Wk selected at the kth step is assessed empirically using surrogate data. Specifically, the CMI I(Yn;Wk|Vk-1) is compared with the 100(1−)th percentile of its distribution computed over Ns surrogates of the last selected component Wk, each generated shifting the original Wk of a randomly selected lag (larger than a minimum lag min set to exclude autocorrelation effects) with respect to Yn and Vk-1. This independent shift generated surrogates with uncoupling between Wk and Yn. Then, if the original CMI is above the percentile threshold, Wk is included in the conditioning vector, otherwise it is discarded and the procedure terminates including k−1 components in the final vector V=Vk-1. Note that, since surrogates are derived only on the basis of the new component Wk, the significance test is actually performed on the entropy term H(Yn|Vk).

TE X Y |Z u k   H (Yn | V Y ,V Z , X nuk 1 , , X nuL )  H (Yn | V Y ,V Z , X nuk , , X nuL )

,

(5)

if u=uk ( i.e. if Xn-uV ). Finally, the aggregate TE results either summing up all lag-specific terms or computing TEX→Y|Z=H(Yn|VY,VZ)−H(Yn|V). In this study, practical estimation of CE from time series data was performed by the classical histogram-based method, that consists in coarsegraining the observed dynamics using Q quantization levels, and computing entropies by approximating probability distributions with the frequencies of occurrence of the quantized values [20]. III. VALIDATION In this Section we test the ability of the proposed estimator of lag-specific TE on simulation examples reproducing different types of interaction between multivariate time series. We considered both linear stochastic and nonlinear deterministic systems interacting at different time lags and with different strength. For each simulation scheme, 100 realizations of the simulated processes were generated, lasting either N=300 or N=600 points. The conditioning procedure was executed using Q=6 quantization levels in the histogrambased entropy estimation and including Lmax=10 (Simulation 1) or Lmax=5 (Simulation 2) past components from each series in the set of candidates. The termination criterion was run over Ns=100 surrogates of the selected component, each generated by circular shift of the component of at least min=20 points [14], and setting =0.05 as statistical significance level. A. Linear Stochastic Processes In the first simulation, we considered M=3 stochastic systems defined by the linear equations: X n  a1 X n 1  a2 X n  2  0.07Yn 1  0.5Z n 1  U n Yn  0.2Z n 1  C1 X n 1  C 2 X n  2  C3 X n  3  Vn , Z n  b1Z n 1  b2 Z n  2  Wn

(6)

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication.

Fig. 2. Overall results on 100 realizations of Simulation 1 performed with parameters 1=1, 2=2, 3=8, C1=C2=C3=0.5, with series length N=300 (a) and N=600 (b). Plots depict the distribution of lag-specific TE (TEX→Y|Z(u); gray boxes: median and 25-75th quantiles, whiskers: 10-90th percentiles; crosses: 5-95th percentiles) and the corresponding number of realizations for which lag-specific G-causality was detected (n(IX→Y|Z(u)>0), gray bars).

where Un, Vn and Wn are independent white noises with zero mean and unit variance. The parameter design in (6) is explicitly chosen to simulate the rhythms and interactions commonly observed in cardiovascular and cardiorespiratory variability [21]. The autoregressive parameters a1, a2 and b1, b2, are set to mimic the self-sustained RA (process Z) and the activity of the so-called Mayer waves in AP (process X). These autonomous oscillations are generated at the typical frequencies for AP and RA oscillations, respectively fx=0.1 Hz and fz=0.3 Hz, and were obtained designing two complexconjugate poles with modulus x,z and phases x,z=±2fx,z; in this study we set x=0.8 and z=0.9 to get the values: 2 a1=2xcosx=1.2944, a2=−x =−0.64, b1=2zcosz=−0.5562, 2 b2=−z=−0.81. The other parameters in (6) identify causal effects between pairs of processes, which are imposed with fixed strength and delay along the directions Z→X, Z→Y, and Y→X, and with strength and delay dependent on the parameters C1,2,3 and 1,2,3 along the direction X→Y. The effects Z→X and Z→Y simulate the known respiratory-related modulations of AP and HP, while the effects X→Y and Y→X simulate the closed loop cardiovascular regulation [22-26]. Fine tuning of all parameters in (6) was performed to reproduce as much as possible the spectral properties of cardiovascular and cardiorespiratory variability series; this is seen in Fig. 1a, showing that the power spectral densities of the simulated SAP, HP and RA series closely resemble those typically observed in real short-term variability [21]. The analysis was performed focusing on the evaluation of lag-specific G-causality from X to Y at varying the coupling and delay parameters. First, we consider the case of Gcausality imposed from X to Y at lags 1=1, 2=2, 3=8 with equal coupling C1=C2=C3=0.5. An example of the analysis performed for a single realization is reported in Fig. 1. The conditioning procedure starts at k=0 with an empty vector of selected components, V0=[·], and a CE value initialized at the entropy of the target series, H(Yn|V0)=1.49 (Fig. 1b). At the first step, the candidate component minimizing the CE is Xn-2, which is selected for entering the conditioning vector because the CMI is above its corresponding threshold for significance; then, the procedure selects progressively the components W2=Xn-8 and W3=Xn-1, and terminates at the fourth step because

4

the component with minimum CE, i.e. W4=Zn-1 such that H(Yn|Zn-1,V3)0, is high for high C and decreases for lower C at lag 1, is high for low C and decreases for higher C at lag 2, and is always very low for the other time lags (Fig. 3, right). These results indicate that the modulus of the lag-specific TE reflects the coupling imposed for specific time delays, and suggest a good sensitivity and specificity in the detection of lag-specific G-causality. B. Nonlinear Multivariate Systems In the second simulation, we considered the system of M=5 coupled Henon maps defined as [12]:

Ym,n  1.4 Ym2,n1  0.3Ym,n2

for m  1, M

Ym,n  1.4  [0.5C(Ym1,n1  Ym1,n 2 )  (1  C )Ym,n1 ]2  0.3Ym,n2

(8)

for m  1, M

where the parameter C modulates the strength of coupling from the (m-1)th and (m+1)th systems to the mth system. In this study we varied C from 0 to 0.8 in steps of 0.1; this range of values assures that realizations of the driven processes do not

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication.

5

Fig. 4. Example of lag-specific TE analysis performed for a realization of Simulation 2 with parameters 1=4, 2=5, and C=0.4. (a) CMI (left) and CE (right) computed progressively by the conditioning procedure; symbols are as in Fig. 1b. (b) Lag-specific TE computed with Y3 as destination process and taking either Y1, Y2, Y4, or Y5 as source process, collecting the remaining four processes in the vector Z; gray and white circles denote respectively uncoupled and coupled lags.

Fig. 3. Overall results on 100 realizations of Simulation 1 performed with time series length N=300, with coupling parameters C1=C, C2=1-C, C3=0, and with delay parameters 1 and 2 such that their absolute difference increased progressively: (a) |1-2|=1; (b) |1-2|=3; (c) |1-2|=5; (d) |1-2|=7. Left panels depict the mean of TE from X to Y computed for the coupled lags (u=1, white circles; u=2, black triangles) and for the uncoupled lags (u≠1, u≠2, gray triangles), as well as the aggregate TE from X to Y (black circles). The right panels depict the number of realizations for which lag-specific Gcausality was detected at lag u=1 (white bars), at lag u=2 (black bars), and at lags u≠1,2 (mean over lags, gray bars).

explode and that complete synchronization is avoided for any pair of processes [12]. In this simulation, the imposed Gcausality relations are Ym-1,n-1→Ym,n and Ym+1,n-2→Ym,n. Different realizations of the simulation were obtained varying the initial conditions and taking series of length N after a transient of 103 points. At each realization, the lags 1 and 2 were chosen randomly between 1 and 5. The analysis was performed assuming Y3 as target process and taking one of the other maps as source process. Fig. 4 reports an example of the analysis performed for a single realization of the simulation of length N=300, with coupling C=0.4 and imposed delays 1=4 and 2=5. As seen in Fig. 4a, in this example the sequential conditioning procedure selects the lagged components of the target process with lags 1 and 2 at the first and fourth steps, and the components of other processes with delays matching the imposed delays (i.e. Y2,n-4 and Y4,n-5) at the second and third step, terminating with an embedding vector that fully corresponds with the terms set in (8) for Y=Y3. As a consequence, the estimated lag-specific TE closely reflects the expected profiles, i.e. TEX→Y|Z(u) is always zero for X=Y1 and X=Y5, while it is nonzero only at lag u=1=4 for X=Y2 and only at lag u=2=5 for X=Y4 (Fig. 4b).

Fig. 5 reports the overall results of lag-specific TE analysis with target process Y3, performed separately for the uncoupled source processes (X=Y1, X=Y5; Fig. 5a,b) and the coupled source processes (X=Y2, X=Y3), in the latter case separating the values relevant to the coupled lags (u=1,2; Fig. 5c,d) from those relevant to the uncoupled lags (u≠1,2; Fig. 5e,f). Along the uncoupled directions, the lag-specific TE was flattened at the zero level, and the corresponding relative frequency of detected G-causality was very low (Fig. 5a,b). The same results were essentially found when considering the TE evaluated along the coupled directions but for the uncoupled lags (Fig. 5e,f), thus demonstrating the high specificity of the approach for this simulation. When computed for the lags at which causality was imposed in the simulation, the TE was distributed at progressively higher values at increasing the parameter C; the corresponding number of significant detected links also increased with C (Fig. 5c,d). To investigate the dependence of lag-specific coupling detection on the significance level adopted for candidate selection, the analysis was repeated for different values of the upper percentile of the surrogate CMI distribution, obtained setting =0.1,  =0.05, and  =0.01. The results obtained over sets of 100 process realizations are reported in Table I in terms of binary classification. As expected, decreasing the significance  led to lower sensitivity and higher specificity, reflecting the tendency of the conditioning procedure to select less components due to the higher threshold. Moreover in this simulation a lower  determined a higher accuracy. However, given that the specificity is always high, a higher  might be preferred to improve the sensitivity to weakly coupled dynamics. The high specificity and the relatively lower sensitivity suggest a conservative behavior of the procedure for component selection, which may be due to the dimensiondependent bias of the adopted histogram-based entropy

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication.

Accuracy Sensitivity Specificity

TABLE I CLASSIFICATION RESULTS FOR SIMULATION 2 N=300 N=600 =0.1 =0.05 =0.01 =0.1 =0.05

=0.01

90.3% 61.0% 93.2%

96.2% 77.5% 98.1%

92.2% 57.8% 95.5%

94.1% 57.3% 97.7%

92.8% 80.9% 93.9%

94.3% 78.1% 95.9%

6

Results are obtained over 100 realizations of Simulation 2, considering all four analyzed directions, nine values of C and five lags. Positives collect the detections at coupled lags (u=1,2, C>0); negatives collect the detections at all uncoupled lags (u≠1,2, all C), plus those at u=1,2 with C=0.

estimator. Moreover, the reported robustness in specificity of the lag-specific TE with respect to  may be explained by the narrow CMI distribution always observed for the resampled CMI at the terminating step compared with the range of variation of CMI (e.g., see Figs. 1,3,6 where the rejection at the last step is always marginal). IV. APPLICATION TO CARDIOVASCULAR AND CARDIORESPIRATORY VARIABILITY

The proposed lag-specific TE analysis was tested on shortterm cardiovascular and cardiorespiratory variability studied during a protocol which is known to have an impact on the cardiovascular control systems, i.e., head-up tilt. The considered dynamical systems were the heart, vascular and respiratory systems, with descriptive processes respectively represented by the variability of HP (process H), SAP (process P) and RA (process R) measured on a beat-by-beat basis. A. Experimental Protocol and Data Analysis Fifteen young healthy subjects (8 males, 25 ± 3 years old), in sinus rhythm and breathing spontaneously, were considered for the analysis. The experimental protocol consisted of 15 minutes of data collection in the resting supine position, followed by head-up tilting performed using a motorized table, and by another 15 minutes of acquisition in the 60° upright position. The acquired signals (1 KHz sampling rate, 12 bit precision) were the surface ECG (lead II), the finger photopletismographic arterial blood pressure (Finapres), and the respiratory nasal flow (by differential pressure transducer). The beat-to-beat time series of HP, SAP and RA were measured offline respectively as the sequences of the temporal distances between consecutive R waves of the ECG, the local maxima of the arterial pressure signal inside each detected HP, and the values of the respiratory signal sampled at the onset of each HP. After removing artifacts and slow trends, two stationary time series of N=300 points were selected for each subject, one in the supine position (~2 min before tilting) and one in the upright position (~2 min after tilting). Weak stationarity was checked for each series verifying the stability of mean and variance in sub-portions the selected window [27]. Each series was normalized to zero mean and unit variance, obtaining the dimensionless series related to HP, Hn, SAP, Pn, and RA, Rn (n=1,...,N). Non-uniform conditioning was performed quantizing the series with Q=6 levels, setting the parameters of the criterion

Fig. 5. Overall results on 100 simulations of Simulation 2 performed with series length N=300 (a,c,e) and N=600 (b,d,f). Plots depict the distribution of lag-specific TE (TEX→Y|Z(u)) and the corresponding number of realizations for which lag-specific G-causality was detected (n(IX→Y|Z(u)>0)), computed as a function of the coupling parameter C. Results are presented grouping X=Y1 and X=Y5, averaging over all lags u=1,2,3,4,5 (a,b); grouping X=Y2 and X=Y4 and averaging over the coupled lags u=1, u=2 (c,d); grouping X=Y2 and X=Y4 and averaging over the uncoupled lags u≠1, u≠2 (e,f).

for termination of candidate selection at Ns=100, min=20 and =0.05, and including Lmax=10 past samples for each series into the set of candidates. As the adopted measurement convention allows instantaneous (i.e., not delayed) effects from RA to systolic SAP and to HP, as well as from SAP to HP, the zero-lag term, i.e. Rn or Pn, was included into the set of initial candidates when appropriate [15]. The statistical significance between the distributions of a given TE measure (aggregate or lag-specific) computed in the supine and upright positions was assessed by the Wilcoxon signed rank test for paired data. For an assigned condition (supine or upright), the statistical significance among the TE distributions computed at different lags was assessed by the non-parametric one-way ANOVA (Kruskall-Wallis test); when the null hypothesis of equal median for all lag-specific TEs was rejected, the statistical significance between pairs of distributions was assessed by a post-hoc analysis (Wilcoxon test) using the Tukey’s honestly significant difference criterion to correct for multiple comparisons. A p