Let us consider two possibly coupled dynamical systems X and Y such that their evolution over time is mapped by the continuous-time stochastic processes X = {X _t } and Y = {Y _t }, which are defined at each continuous-time instant t ∈ R . A well-known undirected measure of the dynamical interaction between X and Y is the mutual information rate (MIR), which quantifies the amount of information exchanged per unit of time by the two processes (Duncan, 1970). If the processes are stationary, the MIR is defined as I ̇ X ; Y = lim τ → ∞ 1 τ I X t − τ : t ; Y t − τ : t , (1) where I (⋅; ⋅) denotes mutual information (MI) and τ is the duration of the temporal window over which the MI is computed. The notation X _t−τ:t denotes the random function expressing the stochastic process evaluated along the time interval of duration τ ending at the time t, also referred to as path (Spinney et al., 2017), i.e. X _t−τ:t = {X _s : t − τ ≤ s < t} (the same holds for the process Y); note that the MI in Eq. 1, and consequently the MIR, are independent on t due to stationarity.

The MIR defined above, as any other information measure applied to continuous-time processes, cannot be readily formulated in terms of the probability mass functions or densities used for discrete and continuous random variables. In continuous time, a viable approach is to establish a generalized form for the information measures via measure-theoretic approaches that unify under one framework the methods specifically developed for discrete and continuous random variables (Spinney et al., 2017). In this framework, information measures can be expressed by using the Radon-Nykodim derivative between appropriate random density functions defined on paths in place of the ratio between probability distributions of random variables adopted in discrete-time (Gray, 2011). The MI measure in Eq. 1 can be then expressed in a generalized form as (Duncan, 1970) I X t − τ : t ; Y t − τ : t = E P ln d P x t − τ : t | y t − τ : t d P x t − τ : t , (2) where the expectation is taken over the path realizations x _t−τ:t and y _t−τ:t of the random functions X _t−τ:t and X _t−τ:t , and the argument of the logarithm is the Radon-Nykodim derivative of two probability measures defined on path functions. With a similar formalism, Spinney and colleagues have formalized different measures of information dynamics for continuous-time processes (Spinney et al., 2017; Spinney and Lizier, 2018). In particular, the transfer entropy rate (TER) from the ‘source’ process Y to the ‘target’ process X is defined as (Spinney et al., 2017) T ̇ Y → X t , τ = lim Δ t → 0 1 Δ t T Y → X t , Δ t , τ , (3) where T Y → X t , Δ t , τ = E P ln d P x t + Δ t | x t − τ : t , y t − τ : t d P x t + Δ t | x t − τ : t (4) is the transfer entropy (TE) formulated in terms of a Radon-Nykodim derivative of conditional probability measures similarly as in Eq. 2 for the MI, and the normalization by the time interval Δt ensures convergence of the TER in the limit of small Δt (Spinney et al., 2017). For stationary processes X and Y, the TER is independent on the time t; moreover, considering realizations of infinite duration yields the constant TER measure T ̇ Y → X = lim τ → ∞ T ̇ Y → X t , τ , (5) which quantifies the rate of information transferred along the causal direction from Y to X. By reversing the role of the two processes, the information transferred along the opposite causal direction can be quantified by the TER measure T ̇ X → Y .

The measures of the rates of information exchanged by X and Y defined in Eqs. 1, 5 are related to each other by a decomposition that expresses the MIR between X and Y as the sum of the TER along the two directions X → Y and Y → X, plus a term related to the instantaneous interaction between the two processes. Specifically, by using information-theoretic rules on Eq. 2 and recognizing Eq. 4 as a conditional MI, i.e., T _Y→X (t, Δt, τ) = I (X _t+Δt ; Y _t−τ:t |X _t−τ:t ), the MIR can be expanded as I ̇ X ; Y = T ̇ X → Y + T ̇ Y → X + I ̇ X ; Y 0 , (6) where the term I ̇ X ; Y 0 = lim Δ t → 0 lim τ → ∞ 1 Δ t I X t + Δ t ; X t + Δ t | X t − τ : t , Y t − τ : t (7) quantifies the rate of information instantaneously exchanged between the two processes conditioned to the knowledge of their past histories. The derivation of the important relation Eq. 6, where all three terms are quantified in [nats/s], is reported in the Appendix.

In this subsection we formulate the computation of MIR for point processes. A point process is a particular class of continuous-time process that is uniquely characterized by a series of indistinguishable events described by their time of occurrence. In a bivariate context, the statistical description of two point processes is provided in terms of the instants marking the event times, i.e., by writing X = {x _i }, i = 1, … , N _X , and Y = {y _j }, j = 1, … , N _Y , where x _i and y _j represent the times of the ith event in X and of the jth event in Y, respectively. For these point processes, the MIR can be computed by leveraging the decomposition provided in Eq. 6 (Mijatovic et al., 2021a) and making the assumption that simultaneous events are not possible, i. e, x _i ≠ y _j , ∀i = 1, … , N _X , j = 1, … , N _Y (Spinney et al., 2017; Mijatovic et al., 2021a; Shorten et al., 2021). This assumption implies that the measure I ̇ X ; Y 0 of instantaneous information exchange between X and Y defined in Eq. 7 is null, so that the MIR between two point processes simply becomes the sum of the two TER terms I ̇ X ; Y = T ̇ X → Y + T ̇ Y → X . (8)

Starting from Eq. 8, the MIR can be calculated by employing methods to define (Spinney et al., 2017) and compute (Shorten et al., 2021) the TER for point processes. Specifically, the TER from Y to X is formulated as T ̇ Y → X = λ ̄ X E p x ln λ X , x i | X x i − , Y x i − λ X , x i | X x i − , (9) where λ ̄ X = N X / T is the average event rate of X, N _X is the number of target events, and T is the duration of the target process; in Eq. 9, λ X , x i | X x i − and λ X , x i | X x i − , Y x i − are the instantaneous event rates of the target process X evaluated at the time of its ith event x _i , respectively conditioned on the history of X and on the histories of both X and Y. In general, the unconditioned instantaneous event rate of the process X, evaluated at the arbitrary time u, is given by λ X , u = lim Δ u → 0 p u N X , u + Δ u − N X , u = 1 / Δ u , where N _X (u) is the counting process that returns the number of events occurred up to time u. At this point it is worth noting that, while the probability p _u is defined at any time point u ∈ R , the expectation in Eq. 9 is taken over the probability p _x of observing a quantity precisely at the time of target events x _i , i = 1, … , N _X (Shorten et al., 2021). This important distinction, upon expressing the conditional event rates in terms of p _u , making a Bayes inversion and noting that lim Δ u → 0 p u ⋅ | N X , u + Δ u − N X , u = 1 = p x ( ⋅ ) , allows to reformulate the expression of the TER as (Shorten et al., 2021) T ̇ Y → X = λ ̄ X E p x ln p x X x i − , Y x i − p u X x i − , Y x i − ⋅ p u X x i − p x X x i − . (10)

Equation 10 shows that the TER depends on the probabilities of the process histories X x i − and Y x i − , evaluated at target events and at arbitrary time points (respectively, p _x and p _u ), whose statistical average is taken only at target events (i.e., over p _x ). The last expression constitutes the basis for the MIR estimation strategy presented in the next subsection.

The approach for MIR estimation, devised in (Mijatovic et al., 2021a; Shorten et al., 2021) and briefly presented in the following, relies on creating history embeddings that cover the past states of the two observed point processes, implementing an operational formulation of Eq. 10 to estimate the TER, and finally using Eq. 8 to obtain the MIR estimate.

In the estimation of the TER from the source process Y to the target process X, the procedure for building history embeddings approximates the past history of the two processes observed either at the times of target events x _i = 1, … , N _X , or at arbitrary time points u _i = 1, … , N _U , sampled in continuous time. In the first case, illustrated in Figure 1A), the history embedding of the target X referred to the event x _i is approximated by taking l inter-event intervals, i.e., X x i − ≈ X x i l = { x i − k + 1 − x i − k , k = 1 , … , l } ; the history embedding of the driver Y referred to x _i is approximated as Y x i − ≈ Y x i l = [ x i − y p , Y y p l − 1 ] , where y _p is the most recent driver event preceding x _i . In the second case (Figure 1B), the histories of both processes as observed from u _i are approximated by taking the interval from the most recent event to u _i followed by l − 1 inter-event intervals, i.e., X u i − ≈ X u i l = [ u i − x p , X x p l − 1 ] , Y u i − ≈ Y u i l = [ u i − y p , Y y p l − 1 ] .

The history embeddings are then used to compute the entropy terms that compose the TER computed according to Eq. 10. Specifically, Eq. 10 can be expressed as T ̇ ̂ Y → X = λ ̄ X H ̂ p u X x i l , Y x i l − H ̂ p x X x i l , Y x i l + H ̂ p x X x i l − H ̂ p u X x i l , (11) where the estimates of the four entropies on the r.h.s. are obtained by approximating the past histories of infinite duration with the l − dimensional history embeddings, and computing the nearest neighbor entropy estimator (Vicente et al., 2011; Faes et al., 2015). Specifically, the terms H ̂ p x ( ⋅ ) and H ̂ p u ( ⋅ ) respectively refer to ‘standard’ differential entropy estimates where expectation is taken over the same probability distribution for which the log-likelihood is estimated, and to ‘cross-entropy’ estimates where the two distributions differ (a detailed procedure is given in Shorten et al. (2021); Mijatovic et al. (2021a)). The entropies are then estimated via the kNN estimator (Kozachenko and Leonenko, 1987), where the parameter k indicates the number of points used for searching the neighbors of each reference point; here, points are realizations of the history embeddings of dimension l or 2l specified in Eq. 11, and the search for neighbors is performed within the set of realizations taken at target events in the case of ‘standard’ entropy estimation, and within a set of realizations observed at arbitrary (randomly sampled) time points in the case of ‘cross-entropy’ estimation. The estimation algorithm, which is described in details in (Mijatovic et al., 2021a; Shorten et al., 2021), proceeds performing neighbor searches and range searches optimized to estimate together the four entropy terms in Eq. 11, in order to achieve compensation of the bias brought by the individual terms to the overall TER estimate. The TER from X to Y is estimated in the same way after reversing the role of the two point processes, and finally the MIR estimate is obtained by simply summing the two TER estimates in Eq. 8.

In this work, we face the issue of estimating the MIR from short realizations of coupled point processes. As any estimate of a measure computed on finite-length realizations of a process, the MIR exhibits bias and variance which typically depend on the system dynamics, the analysis parameters, and the time-series length. While the variance reflects random errors which cannot be corrected, the bias of an estimator is related to systematic errors that can be compensated by knowing the true value of the measure of interest and its average value computed over several repetitions of the analyzed process. However, unfortunately the true theoretical values are generally not known for the MIR of coupled point processes, as analytical results do not exist for the sampling distribution of kNN estimates of entropy quantities. Therefore, here we resort to an empirical procedure that follows previously proposed approaches using surrogate time series to reduce the bias of information-theoretic estimates (Marschinski and Kantz, 2002; Papana et al., 2011). Specifically, first we estimate the bias of the estimator computing its average over several realizations of uncoupled surrogate event series for which the expected MIR is zero, and then we use such average value to correct the MIR estimated on the original coupled processes. While this approach can be theoretically justified as a full correction of the bias only when the true coupling between the processes is zero, it has been shown to provide a reasonable compensation of the bias of coupling and causality measures even for coupled processes (Papana et al., 2011).

The correction procedure adopted in this work is based on the generation of surrogate time series that preserve the individual dynamics of a process while destroying any correlation between pairs of processes. While surrogates are typically used to set a significance threshold in the estimate of coupling measures (Faes et al., 2004; Lancaster et al., 2018), in our approach we do not apply a formal surrogate data test but rather correct the MIR for the bias estimated in the absence of coupling. To do this, after computing the MIR estimate I ̇ ̂ X ; Y for a given realization of two point processes, we generate M surrogate point processes, estimate the MIR over each surrogate pair, and finally compute the corrected MIR (cMIR) as I ̇ ̂ X ; Y c = I ̇ ̂ X ; Y − I ̇ ̂ X ; Y m , (12) where I ̇ ̂ X ; Y ( m ) is the median of the MIR estimated over the M surrogate pairs; we use the median instead of the mean to consider possible deviations of the MIR values from a symmetric distribution. The use of the corrected measure Eq. 12 aims at reducing the bias of MIR in the case of absence of coupling between the two analyzed processes. To generate surrogate data, we adopted the procedure proposed by Shorten et al. (2021) in the context of TER estimation. This procedure implements a local permutation of the patterns forming the history embeddings for the two processes under the null hypothesis of independence of the present of the target and the history of the source given the history of the target. This null hypothesis is related to a more conservative test than that typically performed in TER/MIR estimation; while standard shuffling procedures destroy any relation between the current and past states of the target and the past states of the source, the local permutation test maintains the relation between the target and source histories, by decoupling only the source histories from the target events (Shorten et al., 2021). Nevertheless, to test this approach in comparison with established methods for the generation of surrogate data, we also implemented the algorithm based on random shuffling of the inter-event intervals, which preserves the probability distribution of the series of inter-event intervals; the iterative amplitude-adjusted Fourier transform (IAAFT) procedure (Schreiber and Schmitz, 1996; Perinelli et al., 2020), which preserves both distribution and power spectrum of the intervals; and the JODI algorithm (Ricci et al., 2019; Perinelli et al., 2020), which is specifically designed to preserve amplitude distribution and inter-event autocorrelation in point process data.

In all simulations and real data analyses, we implemented the nearest neighbor entropy estimator by using k = 30 neighbors and the maximum norm to compute distances (Faes et al., 2015), and generating a number of random time points equal to the number of target events (N _U = N _X ) (Mijatovic et al., 2021a). Analyses were repeated varying the length of the history embedding in the range l ∈ {1, 2, 3, 4, 5}. In the simulation study, the dependence of MIR and cMIR on the coupling parameter, type of distribution of the inter-event intervals, and time series length was also analyzed.

In the first simulation, the MIR computed according to Eq. 8, where the two TER terms are estimated as in Eq. 11, was evaluated in pairs of uncoupled point processes by varying the type of inter-event interval distribution of the processes and the distribution parameters. Since for these processes the true value of the index is I ̇ X ; Y = 0 , the values of the MIR estimate I ̇ ̂ X ; Y highlight the bias of the estimator. The results reported in Figure 3 indicate the presence of a negative bias in all simulations, as documented by the negative values of I ̇ ̂ X ; Y measured by varying the type and parameters of the distribution of the uncoupled processes.

For Poisson processes, the bias tends to increase with the event rate and with the mismatch between the rates of the two processes (Figure 3A). For Gaussian processes, the bias increases when the standard deviation of the inter-event intervals is decreased, and is not substantially affected by the mean (Figure 3B). In the case of uncorrelated inverse Gaussian inter-event intervals, the bias is inversely related both to the mean and to the shape parameter of the interval distribution (Figure 3C); the dependence on the shape parameter becomes direct when the inverse Gaussian intervals are correlated in HDIG processes (Figure 3D). Overall, these results indicate that, in the presence of short realizations of point processes as in the present case where N = 300 spikes are simulated, the MIR estimates are strongly biased, and therefore strategies are needed for the compensation of such bias in the practical analysis of the information shared between point processes.

The procedure for compensating the bias of MIR estimates, as well as the performance of the corrected cMIR estimator, are illustrated in Figure 4 for several runs of simulation 2 generated by varying the intensity of the interaction between the HDIG processes modulated by the parameter σ _PAT.

As shown in Figure 4A, the progressive de-coupling of the interactions between X and Y obtained by increasing σ _PAT is reflected by a progressive decrease of the MIR estimates; this behavior is observed for all the analyzed values of the history embedding length l. However, the analysis also confirms the presence of a substantial bias in the estimates of MIR, which take on negative values when the coupling between the two processes decreases. Figure 4B reports the bias-corrected MIR estimates, showing how the correction leads to non-negative values of cMIR even when the processes approach the uncoupled states for high values of σ _PAT. The correction brings the cMIR values in the range 0 − 0.6 nats/s for l = 1, which extends to ∼ 0.7 nats/s for l = 5, and evidences the appropriateness of using higher embedding lengths in the simulated process when the inter-event intervals are modeled by an AR model of order p = 5.

The benefit of longer history embeddings is documented also in Figure 4C, where we employ the standard procedure for testing coupling significance based on surrogate data. This procedure tests the null hypothesis of uncoupling between the two analyzed point processes and is based on generating, from each pair of original realizations of the processes, a suitable number of pairs of surrogate event series using the local permutation method, and then on deeming the original pair as significantly coupled if the MIR value was above the 95th percentile of the MIR surrogate distribution. The percentage of realizations for which the MIR/cMIR values were detected as statistically significant is reported in Figure 4C, showing that the rate of detection of weakly coupled point processes (higher values of σ _PAT) increases for higher embedding lengths.

Figure 5 has the same structure of Figure 4A, and shows alternative approaches to generate the surrogate data consistent with the null hypothesis of uncoupling between the two analyzed point processes. The figure shows that the analysis of cMIR is rather stable at varying the type of surrogate data. The most remarkable difference is that using the shuffling surrogates, similarly to the local permutation surrogates employed in Figure 4A even though with a lower extent, the MIR estimates partially overlap with the distribution of the MIR for the original process realizations when the de-coupling parameter is high (σ _PAT = 210 ms and particularly σ _PAT = 235 ms); such an effect is not observed using IAAFT and JODI surrogates. This suggests that surrogates which preserve autocorrelation properties of the inter-event intervals are more prone to detect weak but significant amounts of information shared by two point processes and to return higher values of the cMIR measure. On the other hand, the use of the local permutation method for the generation of surrogate time series resulted in higher values of the MIR assessed on the surrogates (see the gray areas in Figure 4A vs those in Figure 5). This result is expected, as the local permutation method maintains the relationship of the source history embeddings with the history embeddings of the target, thus allowing to keep a low rate of false positive detection of information transfer (Shorten et al., 2021). Thus the comparison between Figures 4A, 5 evidences that the local permutation surrogates adopted as a main solution in our work tends to favor specificity in the detection of coupled point process dynamics, while surrogates preserving autocorrelation structure of the inter-event intervals tend to favor sensitivity. These considerations are of practical relevance for the analysis of real-world data.

To show that the bias of the MIR estimates is due to the small sample size of the point process realizations analyzed, in Figure 6 we show the MIR computed as a function of the decoupling parameter σ _PAT for different lengths of the simulated processes, N ∈ {150, 300, 1,000, 5,000, 10,000}, together with the cMIR obtained using either the local permutation method or the JODI algorithm to generate surrogate point processes. We observe that increasing the number of simulated events progressively reduces the bias, as documented by the progressively higher values observed for the MIR and by the absence of negative values for N ≥ 5,000. As expected, also the variance of the MIR estimates decreases while increasing N, confirming that larger sample sizes reduce not only the bias, but also the variability of the estimates. We also note that the median of MIR over the surrogate distribution (gray dotted line in Figures 6A,C) is not a constant function of σ _PAT and differs for the two methods for surrogate generation. As a consequence, the cMIR does not represent a simple translation of MIR toward positive values and depends on the adopted surrogates. In particular, the use of local permutation surrogates results in lower values of cMIR compared to that based on JODI surrogates (see Figures 6B,D), suggesting a better bias compensation for the latter approach. Moreover, the non-monotonic behavior of MIR estimated for small sample size (N = 150 and N = 300) is accentuated in cMIR when local permutation surrogates are used (Figures 6A,B), while it is smoothed when JODI surrogates are used (Figures 6C,D).

Figure 7 reports the results of Simulation 3, where coupled HDIG processes are generated so that increasing the variability of the propagation delay from X to Y may also increase the coupling between the two processes. This effect is verified in our simulations by observing that the cMIR measure increases with the parameter σ _PAT, which in this case modulates the variability of both the LF component of the inter-event intervals of X and the propagation delays; the increase of the information shared between the two processes at increasing σ _PAT is observed consistently for all the analyzed history embedding lengths, l ∈ [1, 5].

The analyzed data belong to an historical database previously used to study the effects of physiological stress and cognitive workload on cardiovascular variability (Javorka et al., 2017; Pernice et al., 2019). The data were acquired on 76 young healthy subjects (age: 18.4 ± 2.7 years, 32 males), normotensive and with a normal body mass index (21.3 ± 2.3 kg/m²), and consisted of electrocardiographic (ECG) and blood pressure (BP) recordings acquired synchronously with a sampling frequency of 1 kHz. ECG and BP signals were recorded by using CardioFax ECG-9620 (Nihon Kohden, Japan; horizontal bipolar thoracic leads) and the Finometer Pro devices (FMS, Netherlands; volume-clamp continuous BP measurement), respectively. The experimental protocol foresaw the acquisition of the signals in different physiological states, going from resting conditions to different types of stress (orthostatic or mental). For the analyses carried out in this work, we have taken into account the following states: 1) baseline state (B), with subjects resting in the supine position for 15 min; 2) head-up tilt state (T), obtained by passively tilting the subjects by 45° to the upright position and maintaining them in that state for 8 min in order to produce orthostatic stress; 3) mental arithmetic state (M), obtained with subjects in the supine position and by asking them to sum up as fast as possible 3-digit numbers projected on the ceiling until reaching a 1-digit number and to decide whether the resulting number was even or odd (PMT test, Psycho Soft Software, s. r.o., Brno, Czech Republic), where this task was repeated over a period of 6 min to elicit cognitive load. Further details on the experimental protocol can be found in (Javorka et al., 2017; Pernice et al., 2019).

The data analyzed consisted of sequences containing the timings of the consecutive R peaks in the ECG (event series of the R times) and of the following maxima in the BP signals (event series of the systolic times), previously extracted by means of LabChart 8 (ECG analysis, blood pressure modules) toolbox from ADInstruments (Javorka et al., 2017; Pernice et al., 2019). Moreover, the time series of the RR and PAT intervals were measured respectively as the sequences of the difference between two consecutive R times, and of the difference between each systolic time and the preceding R time. The event series and time series analyzed for each subject and experimental condition consisted of N = 300 events, which were extracted starting respectively ∼8 min after the beginning of the phase B, ∼3 min after the beginning of the phase T, and ∼2 min after the beginning of the phase M; the corresponding RR and PAT time series were checked for stationarity by using a test targeting a restricted form of weak stationarity (Magagnin et al., 2011).

Starting from the interval series of RR and PAT, the mean and standard deviation of the two series, respectively computed as the average interval duration and the interval variability, were computed for each subject and experimental condition. Starting from the corresponding event series of R times and systolic times, the intervals forming the history embeddings were extracted as displayed in Figure 1, and employed as described in Section 2 to estimate first the TER along the two directions of interaction, then the MIR, and finally the cMIR. To test the statistical significance of the differences in the median of the distributions of each measure (mean, standard deviation and cMIR) evaluated across conditions (B, T, M), we used the non-parametric Kruskal-Wallis test, followed by post-hoc paired Wilcoxon signed rank test to assess pairwise differences (B vs. T, B vs. M, T vs. M) with 5% significance and employing the Bonferroni-Holm correction for multiple comparisons.

The results of the real data analysis are summarized in Figure 8, reporting the distributions of the basic statistics (mean and standard deviation of RR and PAT intervals) in the upper panels and of cMIR in the lower panels.

The mean RR interval decreased significantly moving from B to T and from B to M; the effect was more pronounced during tilt than during mental arithmetic. Similarly, both postural stress and mental stress induced a decrease of the variability of the RR intervals, with a larger effect during head-up tilt, as documented by the statistically significant decrease of the standard deviation of the RR intervals moving from B to T and from B to M and by its higher values during M compared to T.

The physiological stressors induced also statistically significant variations in the mean and variability of the propagation delays of the sphygmic wave from the heart to the periphery. Specifically, the mean PAT decreased progressively and significantly while moving from B to T and from T to M, and the standard deviation of PAT increased during T compared to B, and decreased during M compared to T.

The analysis of the cMIR measure indicated that the postural stress tends to increase the information shared between the R times and the systolic times, while mental stress does not have significant effects. In fact, cMIR was significantly higher during T compared to B, and significantly lower during M compared to T when a history embedding l = 1 was used. These variations were less evident when l = 2, as the Kruskal-Wallis test reported statistically significant differences among the three distributions despite the post-hoc tests did not reach statistical significance (B vs T, p = 0.070; B vs M, p = 0.195, T vs M, p = 0.333), and were reduced to non-significant trends when l = 3 and l = 4.

The alterations observed in the basic cardiovascular parameters during the two physiological stressors are in agreement with a large body of literature in cardiovascular variability analysis, and document the involvement of several physiological mechanisms in the elicitation of these stressors. In particular, the lower mean and variability of the RR intervals during tilt and mental arithmetic reflect well-known effects such as the tachycardia and the shift of the cardiac autonomic balance towards sympathetic activation and parasympathetic inhibition induced by postural and mental stress (Montano et al., 1994; Carnethon et al., 2002; Garde et al., 2002; Wood et al., 2002; Martinelli et al., 2005; Javorka et al., 2017, 2018; Kim et al., 2018; Pernice et al., 2019). The interpretation of the shortening of PAT and of the increase of its variability observed during tilt is less straightforward. The PAT is composed by the pre-ejection period (PEP), i.e., the interval from the electrical depolarization of the ventricles to the ejection of the blood from the heart, and by the pulse transit time (PTT), i.e., the time that it takes for the blood pressure wave to reach the body periphery; the PEP depends mainly on the strength of left ventricular contraction, influenced by the Frank-Starling law and by sympathetic control (Krohová et al., 2017), while the PTT is mostly affected by arterial compliance, reflecting (on a short time scale) modulation of blood pressure and vasomotion (Mukkamala et al., 2015; Czippelova et al., 2019). In accordance with our previous research in a related database (Krohová et al., 2017), we expect an increase in PEP during orthostasis as an effect of decreased diastolic filling of the heart via the Frank-Starling mechanism leading to a lower strength of the cardiac contraction. Therefore, the decrease of the mean PAT observed during tilt should reflect mostly a decrease in PTT related to an augmented arterial stiffness caused by peripheral vasoconstriction, which is in turn evoked by the vascular baroreflex response associated with a decrease of blood pressure due to pooling of blood in the lower extremities (Czippelova et al., 2019); the concomitance of these opposite trends (i.e., increase of PEP and decrease of PTT) and the complexity of the related physiological mechanisms including autonomic reflexes and mechanical effects (Rapalis et al., 2017; Czippelova et al., 2019; Pernice et al., 2021) may - together with an increased systolic blood pressure variability associated with tilt - explain the higher variability of PAT observed during tilt. During cognitive load, induced in our protocol by the mental arithmetic task, the more prominent decrease of PAT likely reflects–in addition to vasoconstriction driven by commands stemming from the central nervous system which reduces the PTT–also a reduction of PEP associated with an increased cardiac contractility mediated by the sympathetic nervous system (Martin et al., 2016); in this case, the presence of common trends (i.e., decrease of both PEP and PTT) may explain both the lower PAT and its lower variability measured during mental arithmetic.

According to our results, the physiological mechanisms described above are associated with an increase of the rate of information exchange between the point processes marking the R times of the ECG and the times of arrival of the sphygmic wave in the body periphery. Higher values of MIR are expected when the variations of the propagation delay from one process to another are small, or when such variations occur in phase due to the effect of some common driver mechanism. Since we observe an increase in cMIR simultaneously with a shortening of the mean PAT and an increase of the PAT variability, we conclude that the presence of a common driver oscillation is the mechanism underlying the higher exchange of information. This mechanism was synthetically reproduced in our third simulation (see Figure 7), and can be physiologically explained by the sympathetic activation induced by head-up tilt (Montano et al., 1994; Carnethon et al., 2002; Martinelli et al., 2005). The “common driver” nature of this mechanism can be explained by observing that during postural stress the sympathetic activation is related to the baroreflex mechanism and, as such, it simultaneously involves the variability of the heart period (and thus that of the R times) and the variability of the arterial pressure (and thus that of the PAT) (Porta et al., 2011; Faes et al., 2013), thereby determining a more intense exchange of information between the two processes. In fact, vasoconstriction in the arterioles in systemic circulation is modulated almost exclusively by the sympathetic part of the autonomic nervous system (Krohova et al., 2020) whose oscillations mostly occur in the LF band; a similar effect is mimicked in our simulations in Section 3. On the other hand, the less evident variations of cMIR observed during the mental arithmetic test may be associated with the fact that the sympathetic activation evoked by mental stress is of a different type, likely involving central commands from the upper brain centers (cortex) which control more independently the heartbeat and the arterial compliance without prominent synchronization effects related to the baroreflex (Fauvel et al., 2000).

The observation of statistically significant differences across conditions of the cMIR index only for small values of the history embedding length (variations from B to T and from T to M are detected for l = 1 and, to a lower extent, for l = 2) suggests that the cardiovascular interactions altered by physiological stress occur mostly as a consequence of the variability of the propagation time of the sphygmic wave from the heart to the body periphery, and that the use of longer memory effects may confound the detection of such altered interactions. This result can be expected by considering that the largest part of the analyzed type of interactions is due to the PAT, whose effects are fully captured with l = 1 (note that, within the point process framework, effects explained with l = 1 are not immediate but rather indicative of time-lagged effects with short memory). The result is in agreement with previous observations reporting that the latency of cardiovascular information transfer is typically limited to zero-lag or one-beat interactions, especially during postural stress (Faes et al., 2014). Nevertheless, we remark that the type of cardiovascular interactions studied using time-series based methods (Faes et al., 2014; Porta and Faes, 2015) reflect different mechanisms than those reflected by the event-based method employed here, the former being related mainly to the baroreflex control of heart rate, while the latter being related to blood pulse propagation and arterial contractility.