Disentangling cardiovascular control mechanisms during head-down tilt via joint transfer entropy and self-entropy decompositions

A full decomposition of the predictive entropy (PE) of the spontaneous variations of the heart period (HP) given systolic arterial pressure (SAP) and respiration (R) is proposed. The PE of HP is decomposed into the joint transfer entropy (JTE) from SAP and R to HP and self-entropy (SE) of HP. The SE is the sum of three terms quantifying the synergistic/redundant contributions of HP and SAP, when taken individually and jointly, to SE and one term conditioned on HP and SAP denoted as the conditional SE (CSE) of HP given SAP and R. The JTE from SAP and R to HP is the sum of two terms attributable to SAP or R plus an extra term describing the redundant/synergistic contribution to the JTE. All quantities were computed during cardiopulmonary loading induced by −25° head-down tilt (HDT) via a multivariate linear regression approach. We found that: (i) the PE of HP decreases during HDT; (ii) the decrease of PE is attributable to a lessening of SE of HP, while the JTE from SAP and R to HP remains constant; (iii) the SE of HP is dominant over the JTE from SAP and R to HP and the CSE of HP given SAP and R is prevailing over the SE of HP due to SAP and R both in supine position and during HDT; (iv) all terms of the decompositions of JTE from SAP and R to HP and SE of HP due to SAP and R were not affected by HDT; (v) the decrease of the SE of HP during HDT was attributed to the reduction of the CSE of HP given SAP and R; (vi) redundancy of SAP and R is prevailing over synergy in the information transferred into HP both in supine position and during HDT, while in the HP information storage synergy and redundancy are more balanced. The approach suggests that the larger complexity of the cardiac control during HDT is unrelated to the baroreflex control and cardiopulmonary reflexes and may be related to central commands and/or modifications of the dynamical properties of the sinus node.

A full decomposition of the predictive entropy (PE) of the spontaneous variations of the heart period (HP) given systolic arterial pressure (SAP) and respiration (R) is proposed. The PE of HP is decomposed into the joint transfer entropy (JTE) from SAP and R to HP and self-entropy (SE) of HP. The SE is the sum of three terms quantifying the synergistic/redundant contributions of HP and SAP, when taken individually and jointly, to SE and one term conditioned on HP and SAP denoted as the conditional SE (CSE) of HP given SAP and R. The JTE from SAP and R to HP is the sum of two terms attributable to SAP or R plus an extra term describing the redundant/synergistic contribution to the JTE. All quantities were computed during cardiopulmonary loading induced by −25 • head-down tilt (HDT) via a multivariate linear regression approach. We found that: (i) the PE of HP decreases during HDT; (ii) the decrease of PE is attributable to a lessening of SE of HP, while the JTE from SAP and R to HP remains constant; (iii) the SE of HP is dominant over the JTE from SAP and R to HP and the CSE of HP given SAP and R is prevailing over the SE of HP due to SAP and R both in supine position and during HDT; (iv) all terms of the decompositions of JTE from SAP and R to HP and SE of HP due to SAP and R were not affected by HDT; (v) the decrease of the SE of HP during HDT was attributed to the reduction of the CSE of HP given SAP and R; (vi) redundancy of SAP and R is prevailing over synergy in the information transferred into HP both in supine position and during HDT, while in the HP information storage synergy and redundancy are more balanced. The approach suggests that the larger complexity of the cardiac control during HDT is unrelated to the baroreflex control and cardiopulmonary reflexes and may be related to central commands and/or modifications of the dynamical properties of the sinus node.

INTRODUCTION
Head-down tilt (HDT) is an experimental maneuver inducing an increase of the venous return, central blood volume, and central venous pressure (London et al., 1983;Nagaya et al., 1995). The resulting loading of the cardiopulmonary receptors in the right atrium and pulmonary veins leads to a sympathoinhibitory response increasing the forearm blood flow (London et al., 1983;Nagaya et al., 1995;Tanaka et al., 1999) and decreasing the forearm vascular resistance (London et al., 1983;Nagaya et al., 1995;Tanaka et al., 1999), total peripheral resistance (Nagaya et al., 1995) and efferent muscle sympathetic nerve activity (Nagaya et al., 1995;Tanaka et al., 1999). Since the reduction of venous return and the consequent sinoaortic and carotid baroreceptor unloading during head-up tilt leads to a vagal withdrawal and sympathetic activation (Cooke et al., 1999;Marchi et al., 2013), to an engagement of the arterial baroreflex (Taylor and Eckberg, 1996;Porta et al., 2011) and to a reduction of the importance of the cardiopulmonary pathway (Porta et al., 2012b), it can be hypothesized that the acute central circulatory hypervoleamia induced by HDT produces the opposite cardiovascular response. Indeed, a reduced involvement of the cardiac baroreflex and an improved relevance of the cardiopulmonary circuits during HDT would explain the previously observed increase of respiratory sinus arrhythmia and the decrease of arterial pressure variability especially in the low frequency band (Porta et al., 2014b).
In the field of information dynamics multivariate tools have been recently devised that allow the quantitative description of the dynamical interactions among time series (McGill, 1954;Schreiber, 2000;Barnett et al., 2009;Faes et al., 2011Faes et al., , 2013Faes et al., , 2015Lizier et al., 2011;Wibral et al., 2011Wibral et al., , 2014Chicharro and Ledberg, 2012;Stramaglia et al., 2012;Kugiumtzis, 2013;Porta et al., 2014aPorta et al., , 2015Barrett, 2015). These tools have been successfully applied to disentangle physiological mechanisms from the spontaneous variability of the heart period (HP), systolic arterial pressure (SAP), and respiration (R). For example, the predictive entropy (PE) of HP assessed in the universe of knowledge = {HP,SAP,R}, measuring the decline of uncertainty about the present HP due to the knowledge of the past history of HP, SAP, and R series, was taken as a measure of the loss of complexity of the cardiac neural regulation. Indeed, the PE of HP increased during the vagal withdrawal induced by an orthostatic challenge Porta et al., 2015) and during an experimental maneuver imposing the regularization of the respiratory sinus arrhythmia (i.e., controlled respiration at slow breathing rate; Faes et al., 2015). The transfer entropies (TEs) from SAP to HP and from R to HP in , denoting the portion of PE of HP solely attributable to past values of SAP and R respectively, were taken respectively as measures of the degree of involvement of the cardiac baroreflex and cardiopulmonary reflexes in controlling HP. Indeed, the TE from SAP to HP in gradually augmented when the baroreflex was challenged in proportion to an orthostatic stimulus  and the TE from R to HP in progressively decreased with age (Nemati et al., 2013;Porta et al., 2014a) as a likely consequence of the gradual vagal withdrawal inducing a progressive decoupling between HP variations and respiratory centers (Seals and Esler, 2000;Eckberg, 2003). Unfortunately, the full exploitation of this approach in assessing cardiovascular control is limited by the incomplete decomposition of PE of HP (McGill, 1954;Chicharro and Ledberg, 2012;Stramaglia et al., 2012;Barrett, 2015;Faes et al., 2015). Indeed, usually the full decomposition of PE is given for bivariate interactions (e.g., HP and R)  or limited to the TE term (Chicharro and Ledberg, 2012;Stramaglia et al., 2012;Barrett, 2015), while the full decomposition of the self-entropy (SE) has never been investigated.
The aim of this study is to provide a decomposition of PE including that of both TE and SE terms when the assigned target dynamic is HP and the two exogenous signals are SAP and R. The proposed decomposition is applied to elucidate the response of the cardiovascular control to HDT. We utilized: (i) the PE of HP as a global descriptor of the HP dynamic; (ii) the TEs from SAP to HP and from R to HP in as markers of the degree of involvement of the cardiac baroreflex and cardiopulmonary reflexes, respectively; (iii) the SE of HP as a measure of the information stored into HP; (iv) the part of the SE of HP excluding the contributions of SAP and R as an index of the information stored into HP that cannot be explained by SAP and R; (v) the SAP-R interaction terms as measures of the redundant/synergistic contributions of SAP and R to the information transfer and storage. Results, over the same protocol exclusively relevant to the TEs from SAP to HP and from R to HP in were presented to the EMBC 2014 (Porta et al., 2014b).
In the following we consider HP = y = {y(n), n = 1, ..., N} as an effect series driven by a pair of exogenous series, SAP = x 1 = {x 1 (n), n = 1,...,N}, and R = x 2 = {x 2 (n), n = 1,...,N} where N is the series length and n is the progressive cardiac beat counter. The series y, x 1 , and x 2 form the universe of knowledge = {y,x 1 ,x 2 } about the short-term control of the HP variability exploited in this study (Porta et al., 2012a). The samples of all series in are normalized by subtracting the mean value and by dividing the result by the standard deviation, in such a way that y, x 1 , and x 2 have zero mean and unit variance. We also define the following restricted universes of knowledge obtained from by excluding both exogenous signals, i.e., \x 1 ,x 2 = {y}, or one exogenous signals i.e., \x 1 = {y, x 2 } and \x 2 = {y, The open loop autoregressive (AR) model with two exogenous (X) inputs (ARX 1 X 2 ) describes the dependence of the current value of y, y(n), on past values of the same signal and past values of the exogenous inputs x 1 and x 2 as where a j , b 1j , and b 2j , with 1 ≤ j ≤ p, are the constant coefficients of the regression of y on past values of y, x 1 , and x 2 respectively, p is the order of the regressions, τ 1 and τ 2 are the delays of the action from x 1 and x 2 to y respectively, and w ARX 1 X 2 (n) is zero mean white Gaussian noise with variance λ 2 ARX 1 X 2 . In addition to the ARX 1 X 2 model, four structures, all derived from the ARX 1 X 2 model, are of interest in this study: (i) the ARX 1 model in which the exogenous input x 2 is disregarded by ignoring the regression of y on x 2 ; (ii) the ARX 2 model in which the exogenous input x 1 is disregarded by ignoring the regression of y on x 1 ; (iii) the X 1 X 2 model in which the AR part is disregarded by ignoring the autoregression of y; (iv) the AR model in which only the AR part is considered by ignoring both regressions of y on x 1 and x 2 .
The one-step-ahead prediction of y,ŷ(n/n−1), based on the ARX 1 X 2 model is given bŷ (2) where the coefficientsâ j ,b 1j , andb 2j , with 1 ≤ j ≤ p, are estimated according to an optimization criterion (here the least squares approach minimizing the variance of w ARX 1 X 2 ) (Soderstrom and Stoica, 1988). According to the same optimization criterion the one-step-ahead prediction of the considered simplified versions of the ARX 1 X 2 model (i.e., the ARX 1 , ARX 2 , X 1 X 2 , and AR structures) can be analogously obtained after a new identification of the regression coefficients. The prediction error is defined as the difference between y(n) andŷ(n/n−1). The ability of the model to describe the dynamic of y is quantified by the variance of the prediction error. It is bounded between 0 and the variance of y, σ 2 . Given the normalization of the series, the variance of the prediction error actually ranges between 0 and 1, where 0 indicates perfect prediction (the entire σ 2 is explained by the model), and 1 indicates null prediction (no fraction of σ 2 is explained by the model). In the following we will indicate the variances of the prediction error of the ARX 1 X 2 , ARX 1 , ARX 2 , X 1 X 2 , and AR models as σ 2 ARX 1 X 2 , σ 2 ARX 1 , σ 2 ARX 2 , σ 2 X 1 X 2 , and σ 2 AR respectively.

Definition of the Information-theoretic Quantities
Under the hypothesis of linearity and Gaussianity of the dynamics, the Shannon entropy (ShE) of the normalized series y is ShE y = 0.5 · log(2πe) (McEliece, 2002), where log is the natural logarithm, and the conditional entropy (CE) of y in is CE y = 0.5 · log(2πeσ 2 ARX 1 X 2 ). ShE y and CE y measure, respectively, the total amount of information carried by y and its remaining portion that cannot be resolved using past samples of all signals present in (Barnett et al., 2009). We define the PE of y in (Chicharro and Ledberg, 2012;Faes et al., 2015) as quantifying the portion of uncertainty of y that has been resolved in (i.e., PE y = ShE y − CE y ); the SE of y in \x 1 ,x 2 (Lizier et al., 2012;Wibral et al., 2014) as measuring the part of uncertainty of y that has been resolved in \x 1 ,x 2 , (i.e., solely using past values of y); the conditional SE (CSE) of y given x 1 and x 2  as measuring the part of information carried by y that can be explained in , above and beyond the one that can be resolved using past values of x 1 and x 2 (Figure 1); the conditional selfentropy (CSE) of y given x 1 in \x 2  as measuring the part of uncertainty of y that has been resolved in \x 2 (i.e., without accounting for the possible influences of the exogenous signal x 2 ) above and beyond the one that can be resolved using past values of x 1 ; the cross-entropy (C) from x 1 to y in \x 2  as measuring the fraction of information carried by y that can be explained in \x 2 given past values of x 1 ; the joint TE (JTE) from x 1 and x 2 to y  as measuring the fraction of uncertainty of y that has been resolved in above and beyond the one that can be resolved in \x 1 ,x 2 ; the TE from x 1 to y in (i.e., by accounting for the possible influences of x 2 on y) (Schreiber, 2000;Lizier and Prokopenko, 2010;Kugiumtzis, 2013), denoted here as the conditional JTE (CJTE) from x 1 and x 2 to y given x 2 , as measuring the reduction of uncertainty of y that can be achieved when \x 1 is completed with the introduction of x 1 (Figure 1); the TE from x 1 to y in \x 2 (i.e., without accounting for the possible influences of x 2 on y)  as measuring the reduction of uncertainty of y that can be achieved when \x 1 ,x 2 is enlarged with the introduction of x 1 . By reversing the role between x 1 and x 2 , CSE y|x 2 , C x 2 →y , CJTE x 1 , x 2 →y|x 1 (Figure 1), and TE x 2 →y can be computed as well.
FIGURE 1 | A mnemonic Venn diagram, redrawn from Porta et al. (2015), of the information-theoretic quantities contributing to the decomposition of PE y in = {y,x 1 , x 2 }. The diagram is devised to represent in the information domain the terms whose sum gives the JTE x 1 ,x 2 →y according to Equation (21), i.e., CJTE x 1 ,x 2 →y x 1 (yellow area), CJTE x 1 ,x 2 →y x 2 (yellow area), and ITE x 1 ,x 2 →y (red area), and the terms whose sum gives the SE y according to Equation (22), i.e., CSE y|x 1 ,x 2 (green area), SE x 1 y x 2 (blue area), SE x 2 y x 1 (blue area), and ISE x 1 ,x 2 y (pink area).

Redundant/Synergistic Contribution of x 1 and x 2 to JTE
The interactive TE (ITE) of y is defined as the variation between the sum of the information individually transferred from x 1 to y and from x 2 to y and the quote jointly transferred (McGill, 1954;Stramaglia et al., 2012) (Figure 1). Therefore, the following equality JTE x 1 ,x 2 →y = TE x 1 →y + TE x 2 →y − ITE x 1 ,x 2 →y (11) provides a relation among the information jointly transferred from x 1 and x 2 to y, JTE x 1 ,x 2 →y , the quantities individually transferred from x 1 to y, TE x 1 →y , and from x 2 to y, TE x 2 →y and the ITE x 1 ,x 2 →y . ITE x 1 ,x 2 →y <0 implies synergy of x 1 and x 2 in contributing to JTE x 1 ,x 2 →y , indicating that the information jointly transferred from x 1 and x 2 to y is larger than the sum of the information transferred from x 1 to y and from x 2 to y when x 1 and x 2 are taken individually. Therefore, the ITE x 1 ,x 2 →y measures the redundant or synergistic contribution of x 1 and x 2 to JTE x 1 ,x 2 →y . The smaller and negative the ITE x 1 ,x 2 →y , the more relevant the synergy of x 1 and x 2 in reducing the uncertainty about the present of y above and beyond the contribution of past values of y. ITE x 1 ,x 2 →y >0 implies redundancy of x 1 and x 2 in contributing to JTE x 1 ,x 2 →y , indicating that the information jointly transferred from x 1 and x 2 to y is smaller than the sum of the information transferred from x 1 to y and from x 2 to y when x 1 and x 2 are taken individually. Therefore, the larger and positive the ITE x 1 ,x 2 →y , the more relevant the redundancy of x 1 and x 2 in reducing the uncertainty about the present of y above and beyond the contribution of past values of y. The ITE x 1 ,x 2 →y can be easily estimated as Equation (12) clearly indicates that the information transferred to y in \x 2 might be larger or smaller than that in depending on whether redundancy or synergy occurs. From Equation (12) it is clear that the sign of ITE x 1 ,x 2 →y depends on the balance between two quantities that are larger than 1: (i) σ 2 AR σ 2 ARX 1 ; (ii) σ 2 ARX 2 σ 2 ARX 1 X 2 . If σ 2 AR σ 2 ARX 1 >σ 2 ARX 2 σ 2 ARX 1 X 2 , redundancy is present. Conversely, if σ 2 AR σ 2 ARX 1 <σ 2 ARX 2 σ 2 ARX 1 X 2 , synergy is detected. Since ITE x 1 ,x 2 →y is symmetric in x 1 and x 2 , the abovementioned considerations hold by reversing the role of x 1 and x 2. The ITE x 1 ,x 2 →y can be also expressed as a percent value with respect to JTE x 1 ,x 2 →y . This quantity will be denoted as ITE% x 1 ,x 2 →y .

Redundant/Synergistic Contributions of x 1 and x 2 to SE
The SE y can be seen as represents the synergistic/redundant contribution of x 1 and x 2 to SE y (Figure 1). If SE x 1 ,x 2 y >0, then SE y > CSE y|x 1 ,x 2 . This indicates that past values of y and of the exogenous sources, when taken together, contribute redundantly to resolve the uncertainty of y because the joint consideration of past values of y and of the exogenous sources worsens predictability of y compared to their separated observations. Conversely, if SE x 1 ,x 2 y <0, the joint knowledge of past values of y and of the exogenous sources improves prediction compared to their separated consideration, thus indicating that past values of y and of the exogenous sources contribute synergistically to reduce the uncertainty of y. The SE x 1 ,x 2 y can be also expressed as a percent value with respect to SE y . This quantity will be denoted as SE% x 1 ,x 2 y . Defined the redundant/synergistic contributions of x 1 and x 2 to SE y when x 1 and x 2 are individually considered as and the equality formally corresponds to Equation (11) as far as the information storage in y is concerned, where ISE stands for the interactive SE and represents the redundant/synergetic contributions of x 1 and x 2 to the information stored in y. At difference from Equation (11), where the ITE x 1 ,x 2 →y is the unique synergistic/redundant term, in Equation (17) all parts might be positive or negative depending on whether x 1 and x 2 , taken individually or jointly, contribute redundantly or synergistically to the information storage of y.
By substituting Equations (13), (15), and (16) into Equation (17) it can be easily demonstrated that the thus leading to the following equality that provides a viable estimate of ISE x 1 ,x 2 y .

Full Decomposition of PE of y in
By following the definitions given in the previous sections it can be easily verified that the PE y in can be fully decomposed (Figure 1) as and SE x 2 where SE x 1 y x 2 and SE x 2 y x 1 represent the synergistic/redundant contribution of x 1 conditioned on x 2 and vice versa to SE y in (Figure 1). It is remarkable that the synergistic/redundant contribution of x 1 and x 2 to SE y , measured by SE x 1 ,x 2 y , depends on the balance among the synergistic/redundant contributions of one of the exogenous inputs conditioned on the other, i.e., SE x 1 y x 2 and SE x 2 y x 1 , and on the synergistic/redundant contribution of x 1 and x 2 when they are jointly considered (i.e., ISE x 1 ,x 2 y ). Even though Equation (23) formally corresponds to Equation (21) as far as the information storage in y is concerned, Equations (21) and (23) are structurally very different. Indeed, while all the terms in Equation (23) can be positive or negative, CJTE x 1 ,x 2 →y|x 2 and CJTE x 1 ,x 2 →y|x 1 in Equation (21) are positive (or null).

EXPERIMENTAL PROTOCOL AND DATA ANALYSIS Experimental Protocol
We studied 13 healthy men aged from 41 to 71 years (median: 59 years). A detailed medical history and examination excluded the evidence of any disease. The subjects did not take any medication nor did they consume any caffeine or alcohol containing beverages in the 24 h before the recording. The study adhered to the principles of the Declaration of Helsinki for medical research involving human subjects. The human research and ethical review board of the "L. Sacco" Hospital approved the protocol. All subjects gave their written informed consent. Electrocardiogram (ECG) and noninvasive finger blood pressure (Nexfin, BMEYE, Amsterdam, The Netherlands) were recorded during the experiments. Signals were sampled at 400 Hz. Each experiment consisted of 10 min of baseline recording at rest in supine position (REST) followed by 10 min of recording during HDT with a table inclination of −25 • . Before REST we allowed 15 min of stabilization. The recordings of the HDT session started 5 min after tilting the table. During the protocol, the subjects breathed according to a metronome at 16 breaths . min −1 to prevent modifications of the magnitude of the respiratory sinus arrhythmia owing to the changes of the breathing rate as much as possible (Hirsch and Bishop, 1981;Brown et al., 1993). All experiments were carried out in the afternoon in the same temperature-controlled room and the subjects were not allowed to talk during the protocol. Original data are available through the corresponding author.

Extraction of the Beat-to-Beat Variability Series
After detecting the R-wave on the ECG and locating the Rwave peak using parabolic interpolation, the temporal distance between two consecutive R-wave apexes was computed and utilized as an approximation of HP. The maximum of arterial pressure inside the n-th HP [i.e., HP(n)] was taken as the n-th SAP [i.e., SAP(n)]. R signal was obtained from the respiratoryrelated amplitude modulation of the ECG. The amplitude of the first QRS complex delimiting HP(n), taken as the difference between the R-wave apex and the isoelectric line, was taken as the n-th R [i.e., R(n)]. The occurrences of R-wave and SAP peaks were carefully checked to avoid erroneous detections or missed beats. If isolated ectopic beats affected HP and SAP values, these measures were linearly interpolated using the closest values unaffected by ectopic beats. HP, SAP, and R sequences of 256 consecutive synchronous measures were chosen inside the REST and HDT periods, thus focusing on short-term cardiovascular regulatory mechanisms (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996). The random selection of the onset of analysis within the overall REST and HDT periods made this preprocessing step operator-independent. The series were linearly detrended before multivariate linear regression analysis. If evident nonstationarities, such as very slow drifting of the mean or sudden changes of the variance, were visible despite the linear detrending, the random selection was carried out again. The HP and SAP means and the HP and SAP variances were indicated as µ HP , µ SAP , σ 2 HP and σ 2 SAP and expressed in ms, mmHg, ms 2 , and mmHg 2 respectively.

Power Spectral Analysis
The power spectrum was estimated according to an univariate parametric approach fitting the series according to the AR model (Kay and Marple, 1981). The Levinson-Durbin recursive algorithm was utilized to estimate the coefficients of the AR model and the variance of the white noise. The number of coefficients p was chosen according to the Akaike's figure of merit in the range from 8 to 14. The power spectral density was computed from the AR coefficients and from the variance of the white noise according to the maximum entropy spectral estimation approach (Kay and Marple, 1981). The AR spectral density was factorized into spectral components, the sum of which provides the entire power spectral density (Baselli et al., 1997). The AR spectral decomposition provided power and central frequency of the components of the AR spectral density. The central frequency of the components expressed in cycles . beats −1 was converted into Hz by dividing the values by µ HP . A spectral component was labeled as low frequency (LF) or high frequency (HF) if its central frequency ranged between 0.04 and 0.15 Hz or was in the range of ±0.05 Hz around the paced breathing rate respectively. The LF and HF powers were computed as the sum of the powers of all LF and HF spectral components respectively. The HF power of the HP series, expressed in absolute units (i.e., ms 2 ) and labeled as HFa HP , was utilized as a marker of vagal modulation directed to the heart (Akselrod et al., 1981), while the LF power of the SAP series, expressed in absolute units (i.e., mmHg 2 ) and labeled as LFa SAP , was utilized as a marker of sympathetic modulation directed to vessels (Pagani et al., 1986).

Construction of Surrogates and Surrogate Analysis
We tested the null hypotheses of HP-SAP and HP-R coupling without or with preservation of the HP information storage. This test was performed by creating two sets of surrogates.
The first set was composed by the original SAP and R series, while the HP sequence was substituted with a series obtained by randomly shuffling the HP samples (Palus, 1997). The shuffling procedure was performed according to one of the N! permutations of the HP samples. As a consequence the SAP and R series were fully uncoupled to the HP shuffled dynamics and the original HP information storage was destroyed, while preserving the distribution of all series and the repetitive dynamical structures of SAP and R. This surrogate will be referred to as HP-shuffled surrogate.
The second set was composed by time-shifted versions of the original series (Andrzejak et al., 2003). While the HP series was left unmodified, the SAP and R sequences were shifted according to a delay much larger than the maximal order of the multivariate model (i.e., 50 cardiac beats), thus destroying the short-term temporal correspondence of the SAP and R samples to HP values, while preserving the HP information storage. The delays from SAP and R to HP were independently chosen. The values at the end of the SAP and R sequences were wrapped to their onset. This surrogate will be referred to as time-shifted surrogate.
For each triplet of original HP, SAP, and R series we created a triplet of HP-shuffled and time-shifted surrogates. If the values derived from the original data were significantly different from those obtained from the surrogate sets the null hypothesis of HP-SAP and HP-R coupling without or with preservation of the HP information storage was rejected.
Surrogate data are available through the corresponding author.

Calculation of the Information-theoretic Quantities
The HP series was modeled via ARX 1 X 2 , ARX 1 , ARX 2 , X 1 X 2 , and AR models where X 1 = SAP and X 2 = R. The delays from SAP and R to HP, τ SAP and τ R , were set to 0 to allow the description of the fast vagal reflex (within the current HP) capable to modify HP in response to changes of SAP and R (Eckberg, 1976;Baselli et al., 1994;Porta et al., 2012a). The coefficients were identified via a traditional least squares approach and Cholesky decomposition method (Soderstrom and Stoica, 1988;Baselli et al., 1997). The AR and X parts had the same model order p. The model order was optimized in the range from 4 to 16 according to the Akaike figure of merit for multivariate processes (Akaike, 1974) over the most complex model structure (i.e., the ARX 1 X 2 , model). The whiteness of the HP prediction error and its mutual uncorrelation, even at zero lag, with the SAP and R series was checked over the same model (Baselli et al., 1997;Porta et al., 2012a). All remaining model structures were separately identified using the optimal order of the ARX 1 X 2 model. in both experimental conditions. Statistical analysis was carried out using a commercial statistical program (Sigmaplot, ver.11.0, Systat Software, San Jose, CA, USA). A p < 0.05 was always considered as significant. Table 1 summarizes the results relevant to the time and frequency domain analyses of the HP and SAP series. HDT did not affect the time domain parameters of both HP and SAP series (i.e., µ HP , σ 2 HP , µ SAP , and σ 2 SAP ). Conversely, the frequency domain indexes were significantly modified. Indeed, HFa HP increased during HDT, while LFa SAP significantly declined.

RESULTS
The bar graph shown in Figure 2 compares PE HP computed at REST and during HDT. PE HP was significantly decreased during HDT. The full decomposition of PE HP into more specific quantities is given in Figures 3-6.
The grouped bar graph of Figure 3 shows the two terms forming PE HP according to Equation (20) (i.e., SE HP , white bars, and JTE SAP,R→HP , black bars). SE HP was markedly larger than JTE SAP,R→HP in both conditions. SE HP significantly decreased during HDT, while JTE SAP,R→HP was not affected by the posture modification.
The grouped bar graph of Figure 4 depicts the three constituents of JTE SAP,R→HP according to Equation (21) (i.e., CJTE SAP,R →HP|R , white bars, CJTE SAP,R →HP|SAP , gray bars, and ITE SAP,R→HP , black bars). None of the terms forming JTE SAP,R→HP were modified by HDT. At REST CJTE SAP,R →HP|R , CJTE SAP,R →HP|SAP , and ITE SAP,R→HP were similar, while during HDT ITE SAP,R→HP was significantly smaller than CJTE SAP,R →HP|R . It is remarkable that ITE SAP,R→HP was larger than 0 in all subjects regardless of the experimental condition, suggesting that SAP and R contributed redundantly to JTE SAP,R→HP both at REST and during HDT.
The grouped bar graph of Figure 5 shows the two terms forming SE HP according to Equation (13) (i.e., CSE HP|SAP,R ,  white bars, and SE SAP,R HP , black bars). CSE HP|SAP,R significantly decreased during HDT, while SE SAP,R HP was not affected by the posture modification. SE SAP,R HP was significantly smaller than CSE HP|SAP,R in both experimental conditions. The average value of SE SAP,R HP was larger than 0 indicating that, on average, SAP and R contributed redundantly to the information storage in HP. However, SE SAP,R HP was negative in 31% and 46% of subjects at REST and during HDT respectively, thus indicating that at the level of the information storage in HP synergy between SAP and R occurred frequently.  The grouped bar graph of Figure 6 depicts the synergistic/redundant terms present in the decomposition of SE HP according to Equation (22)  subjects respectively during HDT. These results stressed that at the level of the information storage in HP synergy between SAP and R was commonly present both at REST and during HDT, even though it did not take priority over redundancy. The contributions of the redundant/synergistic terms to JTE SAP,R→HP and SE HP (i.e., ITE SAP,R→HP and SE SAP,R HP ) are compared in Figure 7 after expressing them as ITE% SAP,R→HP and SE% SAP,R HP . Since ITE% SAP,R→HP was larger than 20% at REST, ITE SAP,R→HP represented a sizable amount of JTE SAP,R→HP . Conversely, SE% SAP,R HP was significantly smaller (about 3%), being a negligible quantity compared to SE HP . Both ITE% SAP,R→HP and SE% SAP,R HP were not affected by HDT. Remarkably, while ITE% SAP,R→HP was consistently positive in all subjects, the variability of SE% SAP,R HP was much higher due to the presence of both negative and positive values.
For the surrogate analysis, the PE HP was significantly larger when computed over the original series than over both types of surrogates. This result held in the case of any term of the JTE SAP,R→HP decomposition according to Equation (21). The SE HP and, more specifically, CSE HP|SAP,R was significantly larger over the original data than HP-shuffled surrogates. Conversely, the terms of the SE SAP,R HP decomposition according to Equation (23) computed over the original data were similar to those calculated over both types of surrogates.

DISCUSSION
The methodological findings of this study can be summarized as follows: (i) the study proposes a full decomposition of PE of y in = {y,x 1 ,x 2 } that includes the decomposition of SE of y in addition to the known decomposition of JTE from x 1 and x 2 to y; (ii) both JTE and SE decompositions include a term Frontiers in Physiology | www.frontiersin.org describing redundancy/synergy of x 1 and x 2 in contributing to the information carried by y and the redundant/synergistic term of SE has a more complex structure; (iii) the utility of the JTE and SE decompositions is demonstrated in the field of cardiovascular control analysis to disentangle physiological mechanisms from spontaneous variations and clarify the origin of the increase of respiratory sinus arrhythmia during HDT.
The experimental findings of this study can be summarized as follows: (i) in agreement with the literature we confirm the increase of respiratory sinus arrhythmia after acute cardiopulmonary loading induced by HDT; (ii) the PE of HP decreases during HDT, thus suggesting a larger complexity of the cardiac control and a vagal enhancement; (iii) the SE of HP is larger than the JTE from SAP and R to HP and the SE of HP due to SAP and R is negligible compared to the CSE of HP given SAP and R both at REST and during HDT, thus emphasizing the relevance of physiological mechanisms other than those mediated by SAP and R changes in governing the HP dynamics; (iv) the invariance of the CJTE from SAP and R to HP given R and CJTE from SAP and R to HP given SAP during HDT suggests a limited involvement of the baroreflex and cardiopulmonary pathway in controlling the HP dynamic during HDT; (v) the decrease of the CSE of HP given SAP and R in response to HDT suggests that the increase of respiratory sinus arrhythmia might be the consequence of modifications of the sinus node dynamical properties and/or an enhanced activity of the respiratory centers occurring independently of the cardiac baroreflex and cardiopulmonary circuits; (vi) SAP and R contribute redundantly to the information jointly transferred to HP, while, at the level of the information storage into HP, redundancy and synergy are more balanced; (vii) the amount of synergy or redundancy of SAP and HP to the information transfer and storage in HP is not affected by HDT, thus again stressing that the response to HDT is not mediated by a strong involvement of the cardiac baroreflex and cardiopulmonary circuits; (viii) surrogate analysis indicates that the terms of the JTE decomposition detect significant HP-SAP and HP-R interactions and the CSE of HP given SAP and R measures a significant HP information storage both at REST and during HDT.

Methodological Findings
Decomposition of PE of y in = {y,x 1 ,x 2 } The study proposes a viable approach to tackle the issue of the PE decomposition of a target signal y affected by two exogenous signals, x 1 and x 2 . The PE, measuring the reduction of uncertainty about the present of y when past samples of y, x 1 , and x 2 are given (Chicharro and Ledberg, 2012;Faes et al., 2015), can be decomposed into the information jointly transferred from x 1 and x 2 to y and the information stored into y. The information stored into y estimates the amount of uncertainty about the present of y that can be resolved using only past values of y (Lizier et al., 2012;Wibral et al., 2014), while the JTE from x 1 and x 2 to y quantifies the reduction of uncertainty about the present value of y when past samples of x 1 and x 2 are given above and beyond the information stored into y . The decomposition of JTE from x 1 and x 2 to y into two terms considering the contribution of x 1 to y given x 2 and that of x 2 to y given x 1 plus an additional term describing the balance between redundancy and synergy to the joint information transfer (i.e., the ITE) was originally proposed in (Stramaglia et al., 2012) by generalizing to the conditional case the notion of interaction information (McGill, 1954). Conversely, the full decomposition of SE of y is original, thus completing the PE decomposition of y in = {y,x 1 ,x 2 }. The SE is first decomposed into two terms: (i) the first term is the SE of y conditioned on the exogenous sources (i.e., the CSE of y given x 1 and x 2 ) assessing the information storage not attributable to SAP and R influences; (ii) the second term is a synergistic/redundant term indicating whether past values of y and of both exogenous signals, when jointly considered, contribute synergistically or redundantly to resolve the uncertainty of y (i.e., the SE of y due to x 1 and x 2 ). The SE of y due to x 1 and x 2 can be decomposed further into two terms considering the unique contribution of x 1 and x 2 plus an extra term describing the joint contribution of x 1 and x 2 (i.e., the ISE). We stress that, while in the decomposition of JTE from x 1 and x 2 to y given in Equation (21) the synergistic/redundant term is only one (i.e., ITE), in the decomposition of SE of y given in Equation (22) only one term is definitely larger than or equal to zero (i.e., the CSE of y given x 1 and x 2 ). Indeed, the SE of y due to x 1 and x 2 can be positive or negative and, consequently, the SE of y larger or smaller than the CSE of y given x 1 and x 2 , depending on the balance among synergistic/redundant behaviors of past values of y and x 1 , past values of y and x 2 and past values of y and both x 1 and x 2 in reducing the uncertainty of y.
Both JTE and SE of y in = {y,x 1 ,x 2 } Feature Redundant/Synergistic Terms The JTE decomposition includes a term describing the redundant or synergistic contribution of x 1 and x 2 to JTE (McGill, 1954;Stramaglia et al., 2012;Barrett, 2015;Wibral et al., 2015). This term indicates that the single contributions of x 1 and x 2 to JTE do not invariably contain duplicate information (i.e., redundant contribution) about the present value of y, but, conversely, the contemporaneous knowledge of x 1 and x 2 might lead to extra information (i.e., synergistic contribution). The inclusion of any new additional exogenous signal in the universe of knowledge (e.g., x 2 in \x 2 ) leading to synergy means to improve the prediction of y well above the one obtained when the exogenous signals are taken individually. Through, the SE decomposition this study proves that the notion of synergy or redundancy applies to the information storage as well. Information storage depends on the action of exogenous inputs (Lizier et al., 2012;Wibral et al., 2014) and the significance of the contribution of the exogenous signals to SE of y has been proved experimentally in the context of cardiovascular control analysis . The SE decomposition indicates that x 1 and x 2 might contribute redundantly or synergistically to the information storage. The study of the synergistic/redundant contributions of x 1 and x 2 to the information stored in y is more complex because the assessment of synergy/redundancy should take into account not only the ability of past values of x 1 and x 2 to reduce the uncertainty of y but also that of past values of y, thus increasing the number of synergistic/redundant terms. Conversely, in the case of the information transfer synergy/redundancy is exclusively attributable to past values of x 1 and x 2 since the contribution of past values of y is conditioned out. The multiplicity of terms describing redundancy/synergy at the level of the information storage has been spelled out in Equation (23) and viable estimators for their computation were provided.

JTE and SE Decompositions in Network Physiology
Our approach is framed into the emerging field of network physiology describing the complexity of aggregates of parts and their interactions as a network of nodes with interconnections (Bashan et al., 2012). This feature is in common with other approaches adopting the same logic for representing complex interactions among subsystems regardless of the scale (David et al., 2006;Bressler and Seth, 2011;Bashan et al., 2012;Iatsenko et al., 2013;Kralemann et al., 2014;Stankovski et al., 2015;Porta and Faes, 2016). This description might involve the utilization of raw data (David et al., 2006;Bressler and Seth, 2011), realizations of point processes or series of events (Porta and Faes, 2016) or phase evolutions estimated from raw data or series of events (Iatsenko et al., 2013;Kralemann et al., 2014;Stankovski et al., 2015). The functionals exploited to assess the strength of the interconnections among nodes might be fully adherent to the Wiener-Granger principle (Granger, 1963) if their calculation is based on a direct comparison of indexes computed in the unrestricted and restricted universes of knowledge via metrics assessing the predictability improvement (Bressler and Seth, 2011;Porta and Faes, 2016) and/or uncertainty decrement (Schreiber, 2000;Hlavackova-Schindler et al., 2007;Porta and Faes, 2016), or based on the explicit computation of coupling functions (Iatsenko et al., 2013;Kralemann et al., 2014;Stankovski et al., 2015), or the estimation of coupling coefficients of an assigned model (David et al., 2006). The approach devised in this study is fully consistent with the Wiener-Granger principle in the information domain, where functionals assess the uncertainty decrement and account for conditioning variables according to a multivariate approach. As such, some analogs can be found with fully multivariate methods based on the explicit calculation of coupling functions (Kralemann et al., 2014;Stankovski et al., 2015) and on its decomposition into self-, direct, and indirect components (Stankovski et al., 2015). Nevertheless, in the present study the proposed decomposition is achieved in a completely different framework (i.e., the Wiener-Granger one) and it is expressively devised for the identification of synergistic/redundant components, rather than for the exclusive separation of the direct influences from the indirect ones.

Time and Frequency Domain Analyses of HP Dynamics during HDT
Time and frequency domain parameters confirmed that HDT does not significantly affect the HP and SAP means (Harrison et al., 1986;Nagaya et al., 1995;Kardos et al., 1997;Tanaka et al., 1999) but it increases the HF power of HP (Kardos et al., 1997) and decreases the LF power of SAP (Weise et al., 1995). These findings were interpreted as a sign of the involvement of the autonomic nervous system in adjusting HP and SAP in response to the posture challenge and, more specifically, as an indication of the increased vagal modulation directed to the sinus node and the decreased sympathetic modulation directed to the vessels during HDT (Nagaya et al., 1995;Weise et al., 1995;Tanaka et al., 1999). Unfortunately, time and frequency domain analyses exploited in this study, and traditionally utilized to understand the physiological adaptation to acute central circulatory hypervolaemia (Weise et al., 1995;Kardos et al., 1997), are not helpful to clarify the origin of the increase of respiratory sinus arrhythmia during HDT. This limitation can be primarily attributable to the univariate nature of these classical time and frequency domain analyses and to their inability to interpret causality, thus preventing the possibility of disentangling the HP response to HDT driven by changes of SAP and R from the one independent of them. This limitation is tackled by the proposed multivariate approach grounded in the framework of information dynamics.

Information Dynamics Approach to the Assessment of the Cardiovascular Control during HDT
The involvement of the cardiovascular control in regulating the HP dynamic during HDT is clearly suggested by the significant decrease of PE of HP. Since we reported earlier that the level of predictability of HP based on past samples of HP, SAP and R is under vagal control being increased during head-up tilt and high dose administration of atropine (Porta et al., 2012c), the decrease of the amount of uncertainty about the present of HP that can be resolved by past values of HP, SAP, and R, as measured by the PE of HP, suggests a larger complexity of the cardiac control and an increased vagal regulation during HDT. Even though based on multivariate analysis, this finding is useless in explaining the mechanisms underpinning vagal activation and the increase of respiratory sinus arrhythmia because PE is a global parameter vaguely linked to physiological mechanisms. We need to directly exploit the decomposition of JTE of SAP and R to HP and SE of HP to try to elucidate the origin of the increase of respiratory sinus arrhythmia during HDT.
The information stored into HP, as measured by the SE of HP, is significantly larger than the JTE of SAP and R to HP. This finding suggests that, even though the knowledge of SAP and R is really helpful in resolving the uncertainty of HP, the contribution of these two signals is significantly smaller compared to the ability of past HP values in predicting the current HP. The relevance of the information storage in HP as the likely consequence of the importance of signals driving HP dynamics independently of SAP and R, such as modulation of efferent cardiac neural activity driven by central commands, originating from respiratory and vasomotor centers in the brainstem, and independent of afferent inputs (Preiss and Polosa, 1974;Valentinuzzi and Geddes, 1974;Preiss et al., 1975;Koepchen, 1984;Dick et al., 2009). Also the dynamical properties of the sinus node (i.e., how it responds to changes of sympathetic and vagal inputs) might play an important role in setting the magnitude of the information stored in HP because they directly affect the memory of HP over its past values (Chess and Calaresu, 1971;Berger et al., 1989;Porta et al., 2003). In addition to the dynamical response of receptors, the self-storage of information into HP depends on the type of neurotransmitters, their concentration, rate of release, degradation and removal and sympatho-vagal interactions (Kawada et al., 1996;Nakahara et al., 1998Nakahara et al., , 1999Porta et al., 2003).
Remarkably, we found that the SE of HP and, more specifically, the CSE of HP given SAP and R decreased during HDT, while the JTE of SAP and R to HP remained unmodified. The reduction of the CSE of HP given SAP and R supports the central drive hypothesis as a possible explanation for the increase of respiratory sinus arrhythmia during HDT. Indeed, the activation of a central mechanism, independent of the cardiac baroreflex and cardiopulmonary stimulation, could limit the ability of past values of HP in reducing the uncertainty of the current HP value. For example, if the respiratory centers improved their activity (Valentinuzzi and Geddes, 1974;Dick et al., 2009), the resulting augmented modulation of the cardiovagal motorneuron responsiveness would produce an increase of respiratory sinus arrhythmia (Eckberg, 2003) and make cardiac regulation more complex (Porta et al., 2012c).
The invariance of JTE from SAP and R to HP during HDT is not a sufficient condition to exclude the role of the cardiac baroreflex and cardiopulmonary reflexes in the rise of respiratory sinus arrhythmia during HDT. Indeed, during graded head-up tilt we found that the invariance of JTE from SAP and R to HP hides the progressive increase of the information transferred from SAP to HP in , as quantified by the CJTE from SAP and R to HP given R, and the decrease of the information transferred from R to HP in , as quantified by the CJTE from SAP and R to HP given SAP, with the magnitude of the orthostatic challenge and baroreflex unloading . Therefore, it is necessary to check the trend of CJTE from SAP and R to HP given R and CJTE from SAP and R to HP given SAP during HDT to better characterize the involvement of the baroreflex and cardiopulmonary pathway in controlling HP dynamics. Given the invariance of the CJTE from SAP and R to HP given R and CJTE from SAP and R to HP given SAP during HDT we conclude that the amount of information transferred along the cardiac baroreflex and cardiopulmonary reflexes is not significantly different from that observed at REST and, thus, we exclude again the cardiac baroreflex control and cardiopulmonary reflexes as possible physiological mechanisms underpinning the observed increase of respiratory sinus arrhythmia.
The opposite influences on the venous return, central blood volume and central venous pressure during head-up tilt and HDT, leading to baroreflex unloading, sympathetic activation, and vagal withdrawal in the case of the head-up tilt (Montano et al., 1994;Cooke et al., 1999;Furlan et al., 2000;Marchi et al., 2013) and cardiopulmonary loading and sympathetic inhibition in the case of HDT (Nagaya et al., 1995;Tanaka et al., 1999), might suggest opposite effects on the degree of involvement of the baroreflex control of HP and cardiopulmonary neural circuits. Contrary to this expectation, while head-up tilt maneuver led to an augmented involvement of the baroreflex control of HP and a reduced participation of the cardiopulmonary reflexes (Porta et al., 2012b, the invariance of the information transferred along the cardiac baroreflex and cardiopulmonary pathway observed in the present study suggests that the physiological response to acute central circulatory hypervolaemia during HDT cannot be simply deduced from the knowledge of the response to acute central circulatory hypovolemia during head-up tilt.

Redundancy and Synergy of SAP and R in
Contributing to the Information Transferred and Stored in HP Dynamics SAP and R contribute redundantly to the information jointly transferred to HP. This means that SAP and R hold common information about the present value of HP above and beyond that derived from past values of HP. Remarkably, this quantity is important since it explains more than 20% of JTE from SAP and R to HP. The redundant nature of the SAP and R contributions to the information transferred to HP is not surprising. Indeed, R can directly modulate SAP by modifying venous return, pressure gradients over large arteries in the thorax and stroke volume via respiratory-related changes of the intrathoracic pressure (Innes et al., 1993;Toska and Eriksen, 1993;Caiani et al., 2000). However, the redundant contribution of SAP and R to JTE from SAP and R to HP might also come from more complex interactions and integrations between vasomotor and respiratory centers occurring at the brain stem level. This amount of redundancy might accomplish a principle of fault tolerance and harmonization of neural responses. In this specific experimental protocol the amount of redundancy of SAP and R to the information transferred to HP was not significantly varied during HDT, again confirming that HDT did not affect quantities closely linked to the functioning of the cardiac baroreflex and cardiopulmonary reflexes.
Even though on average SAP and R contribute redundantly to the information storage into HP, we cannot conclude that redundancy is prevailing over synergy as far as the information storage of HP is concerned. Indeed, the SE of HP due to SAP and R, measuring the balance between redundancy and synergy at the level of the information storage at REST, is quite small (i.e., 3% of the SE of HP) and in 31% of subjects synergistic behaviors between SAP and R in contributing to the SE of HP were observed. As a result of the presence of both redundancy and synergy inside the group of subjects in both experimental conditions, indexes describing the synergistic/redundant behavior of SAP and R to the information storage of HP are characterized by greater variability compared to the index describing the synergistic/redundant behavior of SAP and R to the information transferred into HP. The amount of redundancy of SAP and R to the information storage of HP was not significantly varied during HDT, thus again stressing that HDT did not affect quantities linked to the functioning of the cardiac baroreflex and cardiopulmonary circuits even when the action of these reflexes is mediated by memory effects of HP on its own past.

Surrogate Analysis
Surrogate data were constructed with the main aim to test the significance of the proposed indexes as markers of the strength of the HP-SAP and HP-R coupling in absence or presence of a significant amount of information stored in HP. According to this idea two types of surrogates, both destroying the HP-SAP and HP-R coupling are generated. The first type, the HPshuffled surrogates, wiped out the HP autocorrelation function, while the second type, the time-shifted surrogates, preserved it. Remarkably, all indexes derived from the decomposition of JTE from SAP and R to HP both at REST and during HDT were significantly larger from those derived from surrogates, regardless of the type. This result indicates that both at REST and during HDT the HP-SAP and HP-R interactions are significant as well as the detected redundancy of SAP and R in contributing to the JTE from SAP and R to HP, suggesting that indexes derived from the JTE decomposition are helpful to detect physiological interactions from spontaneous variations. The SE of HP and, more specifically, the CSE of HP given SAP and R, was significantly larger in the original data than in the HP-shuffled surrogates both at REST and during HDT. This result suggests that the information stored into the HP dynamics is significant in both experimental conditions. Conversely, the terms of the decomposition of SE of HP due to SAP and R computed over the original data were indistinguishable from those calculated over surrogate data regardless the type of surrogate both at REST and during HDT. We suggest two possible explanations for this finding: (i) the contributions of SAP and R to the SE of HP did not reach the level of significance both at REST and during HDT; (ii) the surrogate analysis proposed in the present study is not suitable to test the significance of the causal interactions from SAP and R to HP at the level of the SE decomposition.

Significance of the Study and Future Perspectives
A full decomposition of the amount of uncertainty about a target signal that can be resolved based on two presumed driving signals is provided. The decomposition is relevant to the information jointly transferred from the two driving signals to the target one and to the information stored into the destination signal.
Terms describing the balance between redundancy and synergy of the two driving series in resolving the uncertainty of the target signal have been highlighted and viable estimators have been proposed. The application to the experimental data suggests the relevance of the approach in dissecting out cardiovascular control mechanisms with the aim of accepting or rejecting physiological hypotheses. Since the proposed quantities are highly specific and take the form of indexes that can be computed very efficiently and robustly via a traditional multivariate regression analysis of spontaneously varying variables, they appear to be suitable candidates for large scale applications to clinical databases recorded even under uncontrolled conditions. Due to the generality of the approach it might be applied not only to cardiovascular physiology and neuroscience, but also in any field of science in which interactions among systems, or constituents of the same system, are under evaluation. Future studies should extend the decomposition to model-free frameworks to account for the possible presence of nonlinear dynamics disregarded by the present approach. In addition, given that in the present contribution the interaction terms actually represent the balance between redundancy and synergy, future studies might test different information decomposition strategies (Barrett, 2015;Wibral et al., 2015) and extend the proposed decomposition of SE to allow the coexistence of both redundancy and synergy as independent positive quantities. We also advocate studies devoted to the improvement of the physical/physiological interpretation of the parts of the JTE and SE decompositions that can be achieved by extending the application of these decompositions to new experimental conditions, proposing new experiments aimed at modulating the terms of the decompositions, comparing this approach to different techniques for the quantification of the coupling strength, developing new strategies for the construction of ad-hoc surrogate sets and designing specific simulation studies.