On the validity of the state space correspondence strategy based on k-nearest neighbor cross-predictability in assessing directionality in stochastic systems: Application to cardiorespiratory coupling estimation

Understanding the causal direction is of paramount importance for the comprehension of the mechanisms governing events and for setting adequate interventions. Thus, tools inferring causality are flourishing, even when considering techniques defined solely in the time domain^1–4 and in very simple situations featuring two interacting systems X and Y. The relevance of assessing causal interactions in bivariate systems comes from the observations that (1) pairwise interactions provide the basic model of causality; (2) the behavior of complex network is frequently analyzed in terms of pairwise interactions; (3) dynamic behaviors of fundamental importance are often reduced in terms of interactions between two systems; (4) due to the limited accessibility and/or restricted possibility of sensing devices, only two observable variables are monitored; and (5) the issue of assessing causality between two systems still presents open questions.

The state space correspondence (SSC) class is a family of methods that has been exploited for the assessment of the causal interactions between two, potentially linked, systems X and Y. The method is grounded on the application of Takens' theorem allowing the separate reconstruction of the dynamic behavior of X and Y via lagged embedding procedure from, respectively, a realization x = {x_n, n = 1, …, N} of an observable variable $X n$ taken from X and followed over time n up to N, and a realization y = {y_n, n = 1, …, N} of an observable variable $Y n$ taken from Y and evolving over n as well.⁵ After reconstructing the dynamic behavior of X and Y, a metric was applied to check the existence of a smooth function mapping the reconstructed behavior of X, taken as a driver, to that of Y, taken as a responder. If a functional relationship exists from X to Y, close points in the reconstructed phase space of X should correspond to close points in the reconstructed phase space of Y and recurrent patterns of X should correspond to recurrent features of Y. The SSC class was originally proposed to assess generalized synchronization and degree of similarity via functionals symmetric under the reversal of the role of X and Y^6–8 and later exploited to infer causality via asymmetry of the drive–response relation when the role of X and Y is reversed.^9–22 The SSC techniques designed to assess the directionality of the interactions can be roughly classified into categories based on cross-predictability assigned sequence length,^9–13 evolution of the cross-prediction while increasing the sequence length,^2,14,15 distance among states of one system conditioned to neighboring states of the other,^15–21 and on the probability of observing recurrent behaviors of one system given recurrences in the activity of the other.²² 

Under the umbrella of the SSC class, the techniques based on cross-predictability^2,9–15 attract a predominant interest in the context of interactions among physiological systems.^{9–13,23–25} However, there are two open questions using SSC methods based on cross-predictability. The first open question is linked to contradictory suggestions on the use of SSC metrics based on cross-predictability because some studies indicated that directionality from X to Y can be decided via the better prediction of the current state of X based on the past states of Y,^2,9,15 while others suggest the use of the better prediction of the current state of Y based on past states of X.^10–14 It has been suggested that this misunderstanding is associated to applications over deterministic non-separable systems where the Granger causality paradigm¹ fails because the effect of the X on Y cannot be canceled after eliminating X from the universe of knowledge.² However, it remains unelucidated which cross-predictability strategy should be utilized in stochastic systems, especially in the presence of weak-to-moderate coupling commonly observed in physiological systems and brain activity interactions.^{9–13,23–25} The second open question is whether the different dynamic properties of driver and responder, e.g., assessed through a marker of irregularity, could bias directionality assessment. It was suggested¹⁶ that, when the activity of the system is assessed via the estimation of its degrees of freedom, the behavior of the more active system tends to be more informative about the dynamics of the less active system than vice versa, as a result of the higher likelihood of observing ambiguous conditions when predicting future behaviors of a complex system using the more limited amount of informative states provided by the less active one, thus biasing the direction of coupling from the more active system to the less active one. Addressing the second open question over stochastic systems should facilitate designing normalization procedures useful to limit the effect of possible biases in the analysis of inherently stochastic interactions such as those among physiological control systems.

The aim of this study is to establish, in the context of nonlinear stochastic processes, the strategy that should be followed to detect causality using the SSC class and to understand whether different dynamical properties of the interacting systems might bias estimation of the causal direction. Among the SSC class, we utilized the established technique²⁶ based on k-nearest neighbor cross-predictability (KNNCP) because it has been exploited in contributions setting the two open questions.^2,9–15 Nonlinear bivariate autoregressive (NBAR) processes will be utilized to simulate gradually varying strength of the action of one stochastic process on the other, while modulating the degree of coupling on the reverse causal direction. Both unidirectional and bidirectional stochastic interactions will be modelled. To identify possible biases in assessing the directionality of coupling, the KNNCP markers will be discussed in connection with indexes assessing the complexity of the two processes.^27–29 The approach is then applied to characterize the interactions between the heart and the respiratory system that occur in a closed loop: indeed, the pathway from respiration (R) to heart period (HP), inducing rhythmical changes of HP at the respiratory rate known as the respiratory sinus arrhythmia (RSA),³⁰ results from the activity of the respiratory centers modulating the vagal activity directed to the sinus node,³¹ while the reverse pathway, namely, from HP to R, usually referred to as cardioventilatory coupling,^32,33 has been suggested by the constancy of the time interval between the heartbeat just preceding the respiratory onset and the respiratory onset.^32–34 Directionality in cardiorespiratory interactions is quantified during experimental maneuvers known to affect cardiorespiratory control and complexity of both HP variability and R such as controlled respiration (CR) and active standing (STAND).

Let us indicate with X and Y two dynamic systems and with $X n$ and $Y n$ two random variables observed from X and Y at time n up to N. The time evolution of $X n$ and $Y n$ is described by the two stochastic processes X = {X_n, 1 ≤ n ≤ N} and Y = {Y_n, 1 ≤ n ≤ N} collecting the evolution of $X n$ and $Y n$ overtime. Let us indicate with x = {x_n, 1 ≤ n ≤ N} and y = {y_n, 1 ≤ n ≤ N} two realizations of X and Y, respectively, acquired during an experimental session designed to observe the functioning of X and Y and their possible interactions. Let us denote with $x n$ the current value of x and with $x n −= [ x n − 1 … x n − m + 1 ]$ the pattern formed by m − 1 past values of x. Analogously, we define $y n$ and $y n −= [ y n − 1 … y n − m + 1 ]$⁠, the equivalent quantities computed over y. According to Takens's theorem,⁵  $x −={ x n −,m≤n≤N}$ and $y −={ y n −,m≤n≤N}$ provide the reconstruction of the dynamic behavior of X and Y in the (m − 1)-dimensional space using the procedure of delay coordinate embedding. In addition, $x −$ and $y −$ provide the totality of patterns $x n −$ and $y n −$ of length m − 1 that can be extracted from x and y, respectively. The SSC class assume that there is a smooth function $f(⋅)$ mapping a point in the reconstructed phase space of X onto a point in the reconstructed phase space of Y and that this function reflects the link from X to Y.

A traditional method²⁶ assessing KNNCP of $y n$ based on $x n −$ is utilized to assure the existence of $f(⋅)$⁠. At odds with the original method,²⁶ immediate actions of x on y were disregarded, thus focusing only on lagged causality. To briefly recall the method, let us indicate with $x n −$ the reference vector and with $y n$ the image of reference vector $x n −$⁠. KNNCP computes the prediction of $y n$⁠, namely, $y ^ n$⁠, as the weighted mean of the images of the k-nearest neighbors of $x n −$ detected in (m − 1)-dimensional embedding space built over x. The weights are the natural exponential of the inverse of their distance from the reference vector.³⁵ The region of the embedding space around $x n −$ including the k-nearest neighbors is set according to the Euclidean norm. It is worth noting that the adopted predictor is equivalent to using the set of time indexes of the k-nearest neighbors of $x n −$ to find the associated vectors in (m − 1)-dimensional embedding space built over y, namely, those featuring the same time indexes as the k-nearest neighbors of $x n −$ and utilizing their images in y to set prediction.^10–14 Since in causality analysis elimination of vectors close in time³⁶ to the reference vector from the set of k-nearest neighbors might affect results, we excluded only the reference vector. The degree of KNNCP of Y given past states of X is measured as the square correlation coefficient ρ² between y and $y ^$⁠, and ρ² is referred to as the cross-predictability function (CPF) of Y based on past states of $X(CP F X → Y)$⁠. $CP F X → Y$ is bounded between 0 and 1, where 0 indicates complete unpredictability of y based on x and 1 suggests that y can be perfectly predicted from x. $CP F X → Y$ depends on the embedding dimension m. The course of $CP F X → Y$ with m was the result of two opposite drifts^12,29: (i) $CP F X → Y$ might increase with m because longer patterns built over x might be more informative about future behavior of y; (ii) at high m, the k-nearest neighbors spread out due to noise and, consequently, $CP F X → Y$ declines toward 0. The maximum of the $CP F X → Y$ over m was taken as a measure of the coupling strength from X to Y^12,29 and referred to as cross-predictability index (CPI) from X to $Y(CP I X → Y)$⁠: the closer to 1, the $CP I X → Y$⁠, the greater the ability to predict the behavior of Y based on X, the stronger the relationship between X and Y.

We reversed the role of X and Y to test the existence of a relationship $g(⋅)$ linking Y to X. We utilized $CP F Y → X$ to assess the degree of predictability of X based on Y. $CP F Y → X$ and $CP I Y → X$ were defined in complete analogy with $CP F X → Y$ when using $y n −$ to predict $x n$ and $CP I X → Y$⁠, respectively. We argue that X → Y is the dominant direction of causality if $CP I X → Y>CP I Y → X$⁠. Conversely, if $CP I Y → X>CP I X → Y$⁠, we conclude that Y → X is prevalent.

The complexity of Y was assessed as the degree of predictability of $y n$ based on $y n −$⁠.^27,28 An approach completely analogous to that exploited for the computation of KNNCP was utilized with the sole difference that the reference vector to predict $y n$ was $y n −$ and the images of the k-nearest neighbors of $y n −$ in y were utilized to compute $y ^ n$⁠. The degree of predictability of Y given past states of Y, measured as the square correlation coefficient ρ² between y and $y ^$⁠, was referred to as the predictability function (PF) of Y based on past states of $Y(P F Y → Y)$⁠. Similarly to CPF, $P F Y → Y$ exhibited a maximum over m and this maximum was taken as a measure of the degree of predictability of Y^28,29 and referred to as predictability index (PI) of $Y(P I Y → Y)$⁠: the closer to 1 the $P I Y → Y$⁠, the greater the ability to predict Y based on its own past states, the lower the complexity of the behavior of Y. Similarly, we defined the degree of predictability of X given past states of $X( PF X → X)$⁠, and its maximum $P I X → X$ measured the level of predictability of X and the degree of regularity of X activity.

The data belong to a historical database designed to evaluate the impact of R and its dominant frequency on HP variability in healthy humans.^38,39 Full details about population and experimental setup were reported in Refs. 38 and 39. The protocol adhered to the principles of the Declaration of Helsinki. The ethical review boards of the “L. Sacco” Hospital, Milan, Italy, approved the protocol (protocol code: 1999-3; date of approval: 1/2/1999). Written informed consent was obtained from all subjects. After a resting period of 10 min necessary to ensure stabilization of the physiological variables, we acquired surface electrocardiogram (ECG) from lead II via a bioamplifier (Marazza, Monza, Italy) and respiratory flow via a nasal thermistor (Marazza, Monza, Italy) in 19 healthy humans (aged from 27 to 35 years, median = 31 years; 11 females). Signals were sampled at 300 Hz through a 12-bit analog-to-digital board (National Instruments, Austin, TX, United States) plugged into a personal computer. Signals were acquired at rest in supine position during spontaneous breathing (SB). We controlled as much as possible physical activity of the subjects and feeding habits the day before the experiment and conditions of the room where the experimental protocol was carried out. The period of SB was followed by CR sessions. Breathing rhythm was imposed by a metronome whose sound was reinforced by the physician's voice. Respiratory rate was randomly selected with the set 10, 15, and 20 breaths/min (CR10, CR15, and CR20). Each subject underwent recordings at different breathing frequencies. All sessions lasted 10 min. All the subjects were familiar with the paced breathing procedure and were able to follow without discomfort the rhythm set by the metronome.

The data belong to ahistorical database designed to evaluate whether the influence of a postural challenge on the link between R and HP could depend on the level of physical training of the subjects.⁴⁰ Full details about characteristics, criteria of inclusion and exclusion and experimental setup were reported in Ref. 40. The protocol adhered to the principles of the Declaration of Helsinki. The study was approved by the Human Research Ethics Committee of the Federal University of São Carlos, São Carlos, Brazil (Protocol: 1.558.731 and 73/2011). Written informed consent was obtained from all subjects. Briefly, we considered 50 apparently healthy men comprising 25 athletes (ATHLETES, aged from 20 to 40 yrs, 30 ± 6 yrs) who practiced cycling for at least 6 uninterrupted months and for at least 150 min per week and non-athletes (NON-ATHLETES, aged from 20 to 40 yrs, 30 ± 6 yrs) who performed physical activity for less than 150 min per week. Peak oxygen uptake assessed during cardiopulmonary exercise test confirm the trained status of ATHLETES compared to NON-ATHLETES. After a resting period of 10 min necessary for ensuring the stabilization of the physiological variables, we recorded ECG from lead MC5 via a bioamplifier (BioAmp FE132, ADInstruments, Australia) and respiratory thoracic movements via a piezoresistive thoracic belt (Marazza, Monza, Italy). Both signals were sampled at 1000 Hz using a commercial analog-to-digital converter (Power Lab 8/35, ADInstruments, Australia). Level of control on the physical activity of the subjects and feeding habits the day before the experiment and conditions of the room where the experimental protocol was carried out were like those of the CR protocol. Signals were recorded for 15 min at rest in supine position (REST) and for additional 15 min during STAND. STAND session always followed the REST one. During the procedure, subjects were advised to breathe normally but were not allowed to talk.

After detecting the R-wave peaks on the ECG using a traditional method based on a threshold on its first derivative, the time distance between two consecutive R-wave peaks was denoted as the nth HP. The jitters in locating the R-wave peak were minimized using parabolic interpolation. The nth value of the R series was derived by sampling the R signal at the first R-wave peak delimiting the nth HP. Detections of the R-wave peaks were visually checked: missed R-waves were manually inserted and erroneous identification was deleted. ECG recordings were free of complex arrhythmic episodes. The effect of isolated arrhythmic beats on the HP series was mitigated through linear interpolation using the two most adjacent HPs computed between sinus beats. The limit of 5% of corrections was never reached.

PI and CPI markers were computed over HP and R series. The series was normalized by subtracting to each value the mean and by dividing the result by the standard deviation, thus obtaining sequences with zero mean and unit variance. All the sequences underwent linear detrending before applying a normalization procedure. We set N = 256, thus focusing on short-term mechanisms governing cardiorespiratory neural interactions.⁴¹ The embedding dimension m was optimized between 1 and 15. The parameter k was fixed to 20 according to our previous experience that suggests limited variability of the results when the value of k is not small compared to the series length.^12,28,29 The parameter τ was fixed to 1 in such a way that the vector of past values corresponds to a pattern of consecutive samples extracted from the series. Since the SSC methods requires stationarity, the sequences, randomly chosen within the experimental sessions, were tested for the stability of the mean and variance.⁴² If the selected series did not pass the test of stationarity of the mean and variance, a new selection was taken until the test was fulfilled.

The significance of the CPI markers was tested against a situation of full uncoupling between the pair x and y generated via a surrogate data approach.²⁰ Surrogate couples were built by a procedure preserving the distribution, power spectral density of the original series, and eventual nonlinear dynamics, but destroying their cross-correlation.²⁰ The surrogate pairs were built as time-shifted versions of the original series. While the cause series was left unmodified, the effect series was shifted according to a delay selected at random between 40 and 160 cardiac beats, thus destroying the short-term temporal correspondence of the two sequences. The values at the end of the effect sequences were wrapped to their onset. We generated 100 surrogate pairs from any original couple. CPI markers were computed over each set of surrogates, and the 95th percentile was extracted. If the marker computed over the original series was above the 95th percentile of the CPI distribution built over surrogates, the null hypothesis of uncoupling was rejected and the alternative hypothesis, namely, the series was significantly associated in the considered time direction, was accepted. The percentage of subjects with a significant CPI was monitored as well and indicated as %CPI. Both %CPI from X to Y $(%CP F X → Y)$ and vice versa $(%CP I Y → X)$ were computed.

As to the CR protocol, two-way repeated measures analysis of variance vs a reference (Holm–Sidak test for multiple comparisons) was utilized to assess the significance of the differences between PI, or CPI, markers within the same experimental condition (i.e., SB, CR10, CR15, and CR20) and between experimental conditions within the same index. The references were markers computed when the target series was HP and during SB, respectively. If the hypothesis for the utilization of the Holm–Sidak test was not fulfilled, the Wilcoxon signed rank test was applied when appropriate. The level of significance of each test was lowered according to the number of comparisons (i.e., 10) to account for the multiple comparison issue. The χ² test (McNemar's test) was applied to the proportion of subjects featuring the rejection of the null hypothesis of uncoupling to assess the effect of CR within the assigned PI, or CPI, marker, and the difference between indexes within the assigned experimental condition. Even in this case, comparisons were set as previously. The level of significance of each test was lowered according to the number of comparisons (i.e., 10) to account for the multiple comparison issue.

As to the STAND protocol after pooling together REST and STAND data, two-way repeated measures analysis of variance (one factor repetition, Holm–Sidak test for multiple comparisons) was utilized to assess the significance of the differences between PI, or CPI, markers within the same group (i.e., ATHLETES and NON-ATHLETES) and between groups within the same index. If the hypothesis for the utilization of the Holm–Sidak test was not fulfilled, nonparametric tests (i.e., the Mann–Whitney rank sum and the Wilcoxon signed rank tests) were applied when appropriate. The level of significance of each test was lowered according to the number of comparisons (i.e., 4) to account for the multiple comparison issue. The χ² approach was applied to the proportion of subjects featuring the rejection of the null hypothesis of uncoupling to assess the between-group difference within the assigned direction of causality (traditional χ² test) and the difference between directionality within the assigned group (McNemar test). The level of significance of each test was lowered according to the number of comparisons (i.e., 4) to account for the multiple comparison issue. The same statistical approach applied to the CR protocol was applied to the STAND protocol when the two experimental conditions were considered separately.

Statistical analysis was performed with a commercial statistical software (Sigmaplot v.14.0, Systat Software, San Jose, CA, USA). The level of statistical significance of all the tests was set to 0.05. A type-I error probability p smaller than the level of significance, eventually corrected according to the number of comparisons, was always taken as significant.

Figure 1 shows the trends of PI [Figs. 1(a)–1(c)] and CPI [Figs. 1(d)–1(f)] in type-I NBAR simulations as a function of d. Results are obtained with b = 0.0 [Figs. 1(a) and 1(d)], b = − 0.3 [Figs. 1(b) and 1(e)], and b = − 0.6 [Figs. 1(c) and 1(f)]. $P I X → X$ and $CP I X → Y$ are represented as black lines, while $P I Y → Y$ and $CP I Y → X$ as red lines. The dashed line represents the mean, and the two solid lines correspond to two standard deviations about the mean computed over the 50 simulated pairs. $P I X → X$ and $P I Y → Y$ were similar regardless of d and b [Figs. 1(a)–1(c)]. The course of PI with d was flat when b = 0.0 and b = −0.3 [Figs. 1(a) and 1(b)]. Only with b = − 0.6, both $P I X → X$ and $P I Y → Y$ exhibited a robust increase with d [Fig. 1(c)]. $CP I Y → X$ remained stable with d when b = 0.0 [Fig. 1(d)] and b = −0.3 [Fig. 1(e)], while it increased with d when b = −0.6 [Fig. 1(f)]. Regardless of b, $CP I X → Y$ rose with d [Figs. 1(d)–1(f)]. $CP I X → Y$ was above $CP I Y → X$ at large d both when b = 0.0 and b = −0.3 [Figs. 1(d) and 1(e)], while the reverse situation was observed at small d when b = −0.6 [Fig. 1(f)].

Figure 2 has a similar structure as Fig. 1, but it shows the trends in type-II NBAR simulations as a function of $c 2$ when $ρ 1= ρ 2=0.60$⁠. As expected, the degree of complexity of X and Y was similar at $c 1= c 2=0.0$ [Fig. 2(a)]. The course of $P I X → X$ with $c 2$ was flat when c₁ = 0.0 [Fig. 2(a)] and c₁ = 0.3 [Fig. 2(b)], while it tended to decrease when c₁ = 0.6 [Fig. 2(c)]. Regardless of $c 1$⁠, $P I Y → Y$ exhibited a tendency to decrease with $c 2$ [Figs. 2(a)–2(c)]. $P I X → X$ was higher than $P I Y → Y$ at high values of $c 2$ with c₁ = 0.0 [Fig. 2(a)] and c₁ = 0.3 [Fig. 2(b)]. The reverse situation was observed at low values of $c 2$ with c₁ = 0.6 [Fig. 2(c)]. The course of $CP I Y → X$ remained stable with $c 2$ at c₁ = 0.0 [Fig. 2(d)], while it exhibited $a$ gradual decrease when both c₁ = 0.3 [Fig. 2(e)] and c₁ = 0.6 [Fig. 2(f)]. Conversely, $CP I X → Y$ increased with $c 2$ regardless of the values of $c 1$ [Figs. 2(d)–2(f)]. $CP I X → Y$ was above $CP I Y → X$ at high values of $c 2$ when c₁ = 0.0 [Fig. 2(e)], while the reverse situation was detected at low values of $c 2$ at c₁ = 0.6 [Fig. 2(f)]. Remarkably, both situations were present with c₁ = 0.3 [Fig. 2(e)]. Trends of PI and CPI shown in Fig. 2 held when $ρ 1= ρ 2=0.85$ as well (not shown), namely, in the presence of type II NBAR simulations with a greater regularity of both X and Y at $c 1= c 2=0.0$⁠.

Figure 3 has the same structure as Fig. 2, but it shows the trends in type-II NBAR simulations as a function of $c 2$ when $ρ 1=0.60$ and $ρ 2=0.85$⁠. As expected, the degree of complexity of X was significantly larger than that of Y, namely, $P I Y → Y$ was above $P I X → X$ at $c 1= c 2=0.0$ [Fig. 3(a)]. The course of $P I X → X$ with $c 2$ was flat when c₁ = 0.0 [Fig. 3(a)], while it exhibited a tendency toward a slow decrease with c₁ = 0.3 [Fig. 3(b)] and c₁ = 0.6 [Fig. 3(c)]. Regardless of $c 1$⁠, $P I Y → Y$ gradually decreased with $c 2$ [Figs. 3(a)–3(c)]. $P I Y → Y$ was above $P I X → X$ at low values of $c 2$ [Figs. 3(a)–3(c)]. The course of $CP I Y → X$ remained stable with $c 2$ at c₁ = 0.0 [Fig. 3(d)], while it decreased both c₁ = 0.3 [Fig. 3(e)] and c₁ = 0.6 [Fig. 3(f)]. Conversely, $CP I X → Y$ increased with $c 2$ regardless of the values of $c 1$ [Figs. 3(d)–3(f)]. $CP I X → Y$ was above $CP I Y → X$ at high values of $c 2$ when c₁ = 0.0 [Fig. 3(d)], while the reverse situation was detected at low values of $c 2$ when c₁ = 0.6 [Fig. 3(f)]. Remarkably, both situations were present with c₁ = 0.3 [Fig. 3(e)].

Figure 4 has the same structure as Fig. 2, but it shows the trends in type-II NBAR simulations as a function of $c 2$ when $ρ 1=0.85$ and $ρ 2=0.60$⁠. Given the opposite setting of the pole moduli when X and Y compared to Fig. 3, $P I X → X$ was above $P I Y → Y$ at $c 1= c 2=0.0$ [Fig. 4(a)]. The trends of $P I X → X$ and $P I Y → Y$ [Figs. 4(a)–4(c)] were similar to those reported in Figs. 3(a)–3(c) with starting points at $c 2=0.0$⁠, completely different reflecting the modifications of the regularity of the processes. Regardless of $c 2$⁠, $P I X → X$ remained above $P I Y → Y$ when c₁ = 0.0 and c₁ = 0.3 [Figs. 4(a) and 4(b)], while $P I X → X$ and $P I Y → Y$ were similar when c₁ = 0.6. Regardless of $c 1$⁠, the course of $CP I X → Y$ increased with $c 2$ [Figs. 4(d)–4(f)]. When c₁ = 0.0, $CP I Y → X$ increased as a function of $c 2$ [Fig. 4(d)], while a slight decline was observed when c₁ = 0.3 and c₁ = 0.6 [Figs. 4(e) and 4(f)]. $CP I X → Y$ was above $CP I Y → X$ at high values of $c 2$ when c₁ = 0.0 [Fig. 4(d)]. At c₁ = 0.3 and c₁ = 0.6, we observed again that $CP I X → Y$ was above $CP I Y → X$ at high values of $c 2$ but also the opposite situation at low values of $c 2$ [Figs. 4(e) and 4(f)].

The vertical error bar graphs of Fig. 5 show PI [Fig. 5(a)] and CPI [Fig. 5(b)] computed in the CR protocol as a function of the experimental condition (i.e., SB, CR10, CR15, and CR20). The black and white bars represent, respectively, $P I HP → HP$ and $P I R → R$ in Fig. 5(a) and $CP I HP → R$ and $CP I R → HP$ in Fig. 5(b). $P I HP → HP$⁠, $P I R → R$⁠, $CP I HP → R$⁠, and $CP I R → HP$ increased during CR10 and CR15 compared to SB, while no significant differences were observed during CR20 [Figs. 5(a) and 5(b)]. $P I R → R$ and $CP I HP → R$ tended to be higher than, respectively, $P I HP → HP$ and $CP I R → HP$ [Figs. 5(a) and 5(b)]. This result was significant during SB, CR15, and CR20 in the case of PI [Fig. 5(a)] and during CR15 and CR20 in the case of CPI [Fig. 5(b)].

The vertical bar graphs of Fig. 6 show %CPI computed in the CR protocol as a function of the experimental condition (i.e., SB, CR10, CR15, and CR20). The black and white bars represent, respectively, $%CP I HP → R$ and $%CP I R → HP$⁠. The percentage of rejection of the null hypothesis of uncoupling did not vary either across experimental conditions or coupling directions.

Figure 7 has the same structure as Fig. 5 but shows PI [Fig. 7(a)] and CPI [Fig. 7(b)] computed in the STAND protocol as a function of the group (i.e., ATHLETES and NON-ATHLETES). Data were pooled together regardless of the experimental condition. $P I HP → HP$⁠, $P I R → R$⁠, $CP I HP → R$⁠, and $CP I R → HP$ did not vary across the groups [Figs. 7(a) and 7(b)]. $P I HP → HP$ and $CP I HP → R$ were significantly higher than, respectively, $P I R → R$ and $CP I R → HP$ [Figs. 7(a) and 7(b)].

Figure 8 has the same structure as Fig. 6 but show %CPI computed in the STAND protocol as a function of the group (i.e., ATHLETES and NON-ATHLETES). Data were pooled together regardless of the experimental condition. The percentage of rejection of the null hypothesis of uncoupling did not vary either across groups or coupling directions.

Figure 9 has the same structure as Fig. 7 but shows PI [Figs. 9(a) and 9(c)] and CPI [Figs. 9(b) and 9(d)] in ATHLETES [Figs. 9(a) and 9(b)] and NON-ATHLETES [Figs. 9(c) and 9(d)] as a function of the experimental condition (i.e., REST and STAND). In ATHLETES, $P I HP → HP$ increased during STAND compared to REST, and during STAND, it was higher than $P I R → R$ [Fig. 9(a)]. The same results held in NON-ATHLETES [Fig. 9(c)]. In ATHLETES, $CP I HP → R$ was larger than $CP I R → HP$ at both REST and during STAND, and solely $CP I HP → R$ decreased during STAND compared to REST [Fig. 9(b)]. Similar conclusions held in NON-ATHLETES, even though the difference between $CP I HP → R$ and $CP I R → HP$ during STAND was not significant [Fig. 9(d)].

Figure 10 has the same structure as Fig. 8 but show %CPI in ATHLETES [Fig. 10(a)] and NON-ATHLETES [Fig. 10(b)] as a function of the experimental condition (i.e., REST and STAND). In ATHLETES, the percentage of rejection of the null hypothesis of uncoupling did not vary across coupling directions, but STAND reduced solely $%CP I HP → R$ compared to REST [Fig. 10(a)]. In NON-ATHLETES, the percentage of rejection of the null hypothesis of uncoupling did not vary across either coupling directions or experimental conditions [Fig. 10(b)].

The main findings can be summarized as follows: (i) results of the simulations suggest that the SSC method based on KNNCP can detect causality in the context of stochastic NBAR processes as the time direction setting the better predictability from the past of the driver to the future state of the target; (ii) simulations suggest that this ability holds in any condition of complexity of the interacting series; (iii) slowing the breathing rate increases the strength of the causal relationship in both temporal directions in a balanced modality; (iv) STAND is more powerful in modulating the coupling strength on the pathway from HP to R; (v) regardless of the protocol and experimental condition, the strength of the link from HP to R is stronger than the one of the reverse pathway; (vi) significant causal relationships in both temporal directions are found regardless of the level of complexity of HP variability and R.

Controversial statements are present in literature about the detection of causality from X to Y via the SSC method based on KNNCP. Table I summarizes contributions using the SSC technique based on KNNCP for causality detection and lists the criterion to decide the dominant causal direction from X to Y. It was suggested that the criterion for detecting causality from X to Y is the better predictability of Y using past values of X,^10–14 while the same conclusion was reached in the presence of a better predictability of X using past values of Y.^2,9,15 These controversial suggestions were given according to simulations of nonlinear deterministic systems^2,9–12,15 and motivated by the need of facing situations of non-separability typical of nonlinear deterministic systems² and different complexity of the driver and responder behaviors.^9,16 However, given that the criterion about the use of the dominant predictability of X based on Y to detect causality from X to Y is validated only over simulations generated by nonlinear deterministic systems,^2,9 the issue of elucidating whether this criterion is valuable in stochastic systems is relevant. Addressing this issue is important because the abovementioned criterion was applied in typical applications in which dynamical interactions are unlikely to be fully deterministic and are more likely to be inherently stochastic such as in the case of the assessment of neural activity interactions.^9,23–25 Remarkably, in the context of stochastic interactions, a paradigm exploiting the opposite information flow, namely, Granger causality,¹ is commonly applied. In addition, it remains to be established whether conclusions might be biased by the different dynamic properties of X and Y^16–21 because a less complex system is more prone to produce extra ambiguities when using it to predict the behavior of the more complex one, and this effect might bias the detection of directionality in different way according to the strategy followed to decide causality,^{2,10,11,14,15} namely, whether X is mapped onto Y or vice versa. This issue was again suggested in the context of nonlinear deterministic systems,¹⁶ but its relevance when stochastic interactions are considered was not elucidated. As an additional difficulty that can make interpretation of directionality more cumbersome is that it depends on whether cross-unpredictability⁹ or cross-predictability^{2,10,11,14,15} was assessed. Indeed, metrics detecting the dominance of the causal direction from X to Y according to the better ability of X in predicting Y (i.e., cross-predictability of Y on X) are equivalent to those based on the larger cross-unpredictability of X given Y than vice versa.

It can be argued that in stochastic systems, Granger causality paradigm is sufficient to detect causality.¹ However, Granger causality test is based on the predictability improvement when the restricted universe of knowledge comprising the effect and all cofounding factors is completed with the cause. Such a definition requires that the cause carries a unique contribution to the effect that cannot be derived from the dynamic behavior of the effect per se. The SSC method based on KNNCP does not rely on the principle of unique contribution of the cause to the effect because eventual self-contributions of past values of the effect on future dynamics of the effect are disregarded. As such it might be useful in testing causality in contexts where past samples of the cause might inflate memory of the target leading to null values of Granger causality in the presence of significant coupling between the driver and the responder. This situation is more likely in control physiological systems exhibiting autoregulatory mechanisms and reflex pathways depending on the activity of sensors measuring controlled variables such as in cardiorespiratory interactions.⁴³ 

Thus, among the SSC class, we assessed the KNNCP technique over stochastic NBAR processes designed to gradually vary the strength of the forward action (i.e., from X to Y) while modulating that of the backward one (i.e., from Y to X). The two components of the NBAR processes were devised to feature equal or different complexity in the uncoupling condition to favor testing directionality when the driver and target behaviors have balanced, or unbalanced, complexity. Simulations were devised to indicate whether the SSC method based on KNNCP could be exploited in the context of stochastic systems and in association with the criterion detecting causality from X to Y when the predictability of Y based on past values of X is better than that of X based on past values of Y.

While increasing d in type-I NBAR simulations (Fig. 1) and $c 2$ in type-II NBAR simulations (Figs. 2–4), the strength of the coupling from X to Y rises, and it is expected that the causal direction from X to Y could become dominant, at least when d and $c 2$ are large. This expectation is confirmed by the CPI markers given that $CP I X → Y$ increases with d [Figs. 1(d)–1(f)] and $c 2$ [Figs. 2, 3, 4(d), 4(e), and 4(f)], and $CP I X → Y$ is significantly above $CP I Y → X$ at high values of d [Figs. 1(d)–1(f)] and $c 2$ [Figs. 2, 3, 4(d), 4(e), and 4(f)]. Depending on the setting of the parameters b and $c 1$⁠, the strength of the coupling on the reverse causal direction (i.e., from Y to X) is expected to be dominant over the causal direction from X to Y at small d and $c 2$⁠. This expectation is confirmed by the CPI markers in all the simulations because $CP I X → Y$ is significantly below $CP I Y → X$ at values of d close to 0.0 when b = −0.6 [Fig. 1(f)] and at $c 2$ close to 0.0 when $c 1=0.6$ [Figs. 2, 3 and 4(f)]. Remarkably, in situations where the strength of the coupling from Y to X is intermediate, namely, when b = −0.3 [Fig. 1(e)], and at $c 1=0.3$ [Figs. 2, 3 and 4(e)], we observed that $CP I X → Y$ is significantly below $CP I Y → X$ at $c 2$ close to 0.0, indicating a prevalent causal direction from Y to X [Figs. 2, 3 and 4(e)], while the opposite situation was visible at d [Fig. 1(e)] and $c 2$ [Figs. 2, 3 and 4(e)] close to 1.0, indicating the dominance of the reverse causal direction (i.e., from X to Y). These results stress that when considering stochastic NBAR processes, the better prediction of the target behavior based on that of the driver could be utilized for the detection of the causal direction from the driver to the target,^10–14 while the opposite setting to decide the same causal direction, namely, the use of the past of the target to condition the behavior of the driver,^2,9,15 should be limited to applications where determinism is evident.

Directionality was tested in the presence of no significant differences in the dynamical properties of X and Y, as it is suggested by the similar predictability levels [Figs. 1 and 2(a)–2(c)] and when remarkable differences between the predictability of X and Y [Figs. 3 and 4(a)–4(c)]. The expected causal direction was detected when complexity was similar, or different, and in both situations when the complexity of the target was larger than that of the driver [e.g., in Figs. 3(b) and 3(e) at small $c 2$ and Figs. 4(a) and 4(d) at large $c 2$] as well as in the opposite situation [Figs. 4(b) and 4(e) at small $c 2$]. Therefore, results of our simulations support the notion that the SSC method based on KNNCP could detect the expected direction of causality regardless of the dynamical properties of X and Y, as measured via PI. We stress that this conclusion holds for stochastic processes like those utilized in the present study, while it might not work in deterministic systems,^16–19 thus requiring additional normalization procedures.^20,21

Cardiorespiratory coupling is defined as the series of physiological mechanisms that allows a certain degree of coordinated activity between the heart and the respiratory system.⁴³ The main consequence of cardiorespiratory interactions is the rhythmical decrease in HP during inspiration and increase during expiration, usually termed RSA.³⁰ Several physiological mechanisms contribute to the link from R to HP: (i) respiratory centers modulating responsiveness of vagal motoneurons impinging the sinus node;³¹ (ii) the venous return-stroke volume synchronization leading to an inverse relationship between respiratory changes of the venous return to the right atrium and left ventricular filling, or stroke volume;^44,45 (iii) the cardiac arm of the baroreflex converting modifications of systolic arterial pressure at the breathing rate into HP fluctuations;^46,47 (iv) respiratory modifications of the intrathoracic pressure modulating venous return and activating low pressure stretch receptors present in the atria;⁴⁸ (v) the Hering–Breuer reflex activated by pulmonary stretch receptors during lung inflation.⁴⁹ The relevance of the link from R to HP is mainly under vagal control, being limited in situations of vagal withdrawal,^50–52 but also the sympathetic control can contribute given that sympathetic bursts are inhibited in the late half of inspiration and early phase of expiration⁵³ and facilitated when diastolic arterial pressure hits the nadir during the respiratory cycle via the activation of the sympathetic arm of the baroreflex.^54,55 In addition to the link from R to HP, even the reverse pathway has been repetitively proven.^32–34 The link from HP to R accounts, usually referred to as the cardioventilatory coupling,^32,33 for the possibility that cardiac activity could contribute to the respiratory rhythm. The cardioventilatory coupling accounts for the observation that the time interval between the cardiac beat just preceding the inspiratory onset and the inspiratory onset is independent of HP mean and with a tight distribution centered about 0.5 s.³⁴ The relevance of this pathway has been suggested via model-based causality analysis as well.^56–58 The more likely explanation of the link from HP to R is that the activation of the fast vagal arm of the cardiac baroreflex could facilitate the initiation of the inspiratory phase by acting on the respiratory centers.^59,60 The complexity of the HP-R dynamic relationship explains why HP beat-to-beat fluctuations and R are correlated at the respiratory rate with a degree of association usually significant, but below 1,^40,52,61,62 and stress the necessity to disentangle closed-loop interactions via a directional approach.⁶³ 

It is well known that RSA and breathing regularity increase while decreasing the respiratory rate.^38,62,64,65 Even the degree of association between HP and R at the respiratory rate, as measured via squared coherence or a synchronization marker, increases while slowing the breathing rate,^38,62 and this rise might be responsible for the increased regularity of HP.³⁸ The SSC method utilized in the present study confirmed the dependence of complexity of HP and R on the breathing rate that has been suggested in the information domain via model-free approaches.³⁸ In addition, we observed an increase in the coupling strength from R to HP and vice versa while slowing the breathing rate, and this finding makes more specific the notion that slowing the respiratory frequency increased the strength of the association between R and HP at the respiratory rate.⁶² Results about the coupling strength from R to HP and vice versa stress that any analysis involving the computation of cardiorespiratory coupling markers carried out over different groups and/or experimental conditions must be performed at similar respiratory rates. In addition, the present study suggests that the coupling strength along both causal directions can be modulated by varying the respiratory rate. As an additional finding of the present analysis, the coupling strength from HP to R was larger than that along the more usual time direction exploited in modeling applications,^50–52 thus suggesting the stronger effect of cardioventilatory coupling compared to the relevance of the autonomic pathway comprising several different reflexes all targeting HP. This result might suggest the greater complexity of the pathway from R to HP compared to the reverse causal link resulting from several reflexes achieving different goals. However, the stronger impact of cardioventilatory coupling compared to that of the autonomic influences directed to the sinus node did not produce any significant difference on the percentage of subjects exhibiting the rejection of the null hypothesis of HP-R uncoupling that remains extremely high. Therefore, this result supports the remarkable presence of closed-loop interactions between R and HP regardless of the type of breathing (i.e., spontaneous or paced) and respiratory rate. Since the coupling strength from HP to R was generally higher than that in the reverse direction and this result was detected in the presence of the greater complexity of the driver compared to the target, we cannot exclude that this configuration might contribute to the dominance of the link from HP to R in the CR protocol. However, the surrogate data approach indicated that, although the relationship from HP to R is stronger than the opposite one, the causal direction from R to HP is significant as well, thus suggesting that the lower complexity of R compared to HP does not prevent the detection of significant causal interactions from R to HP.

It is well known that RSA in ATHLETES is larger than that in NON-ATHLETES,^66,67 as a likely result of increased vagally-mediated influences,^68,69 even though the dependence of the magnitude of HP variability on the mean HP might distort the perception of the relevance of the phenomenon.^70,71 Remarkably, this result is not associated with a higher HP-R squared coherence function.⁴⁰ In agreement with the CR protocol, in the STAND protocol, we observed that the strength of the causal direction from HP to R is larger than that in the reverse causal direction. Again, like in the CR protocol in the STAND one, the null hypothesis of HP-R uncoupling was rejected in most subjects on both the causal directions, and this conclusion holds in ATHLETES and NON-ATHLETES. Remarkably, in the STAND protocol, the higher strength of the relationship from HP to R compared to the reverse one was observed in the opposite situation of complexity of HP series and R compared to the CR protocol, namely, with the complexity of R larger than that of HP variability, thus stressing again the independence of the conclusions about directionality on the complexity of the driver and target dynamics. Conclusion about the dominance of the strength of the link from HP to R compared to the opposite causal direction is confirmed after separating data recorded at REST from those during STAND. It is worth noting that this finding held in both ATHLETES and NON-ATHLETES, and even when the complexity of the HP series was lower than that of R during STAND. Sympathetic activation and vagal withdrawal induced by STAND, known to limit the complexity of HP variability,^67,72 is associated with a decrease in the percentage of the rejection of the null hypothesis of HP-R uncoupling, and this tendency is significant, especially in ATHLETES, in the time direction from HP to R, thus stressing the possibility of modulating HP-R dynamic interactions even in a separate way according to their causal directions.

The SSC method based on KNNCP can be exploited to assess the directionality of coupling in the context of bivariate applications of time series with stochastic properties. The better KNNCP of the behavior of Y given the activity of X compared to the reverse situation can be utilized to detect the dominant causal action from X to Y. This conclusion is at odds with that proposed in deterministic systems,^2,10 and this criterion should be preferred in applications featuring stochastic interactions such as in assessing physiological control mechanisms and brain activity interactions. Directionality assessment appears to be reliable regardless of the complexity of the dynamics of the driver and target systems. Applications to the characterization of the dynamic interactions between the heart and the respiratory system based on spontaneous variations of HP and R activity suggest that they interact in closed loop and the strength of the dynamic link on both causal directions, namely, from HP to R and vice versa, can be modified by slowing respiratory rate and modifying posture. We conclude that the SSC method based on KNNCP can be fruitfully utilized to describe cardiorespiratory control and its modifications imposed by physiological maneuvers with promising applications in physiology, sports medicine, and clinics. Future study should test whether conclusions drawn through the SSC method based on KNNCP held even when using SSC metrics grounded on cloud size,^16–21 conditional probability of recurrence,²² transfer entropy markers²⁹ or different indexes measuring the dynamical properties of X and Y.^73,74

This work was partially supported by the Italian Ministry of University and Research (MUR) via the Progetto di Rilevante Interesse Nazionale (PRIN) 2022SLB5MX to Alberto Porta (CUP code: G53D23000410006; UGOV code: PRIN202223APORT_01) and by the Italian Ministry of Health via the Ricerca Corrente program to the IRCCS Policlinico San Donato.

The authors have no conflicts to disclose.

Alberto Porta: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Resources (lead); Software (lead); Supervision (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Raphael Martins de Abreu: Data curation (lead); Writing – review & editing (supporting). Vlasta Bari: Formal analysis (supporting); Writing – review & editing (supporting). Francesca Gelpi: Formal analysis (supporting); Writing – review & editing (supporting). Beatrice De Maria: Writing – review & editing (supporting). Aparecida Maria Catai: Data curation (supporting); Writing – review & editing (supporting). Beatrice Cairo: Data curation (supporting); Formal analysis (supporting); Supervision (supporting); Validation (supporting); Writing – review & editing (supporting).

The data relevant to the CR protocol have been utilized for the European Study Group on Cardiovascular Oscillations (ESGCO) 2022 challenge and is available in “ESGCO 2022 challenge: Characterizing cardiorespiratory interactions from spontaneous fluctuations of heart period and respiratory flow during controlled breathing” at https://www.esgco.org/challenges. The data that support the findings of this study are available from the corresponding author upon reasonable request.