On the Different Abilities of Cross-Sample Entropy and K-Nearest-Neighbor Cross-Unpredictability in Assessing Dynamic Cardiorespiratory and Cerebrovascular Interactions

Nonlinear markers of coupling strength are often utilized to typify cardiorespiratory and cerebrovascular regulations. The computation of these indices requires techniques describing nonlinear interactions between respiration (R) and heart period (HP) and between mean arterial pressure (MAP) and mean cerebral blood velocity (MCBv). We compared two model-free methods for the assessment of dynamic HP–R and MCBv–MAP interactions, namely the cross-sample entropy (CSampEn) and k-nearest-neighbor cross-unpredictability (KNNCUP). Comparison was carried out first over simulations generated by linear and nonlinear unidirectional causal, bidirectional linear causal, and lag-zero linear noncausal models, and then over experimental data acquired from 19 subjects at supine rest during spontaneous breathing and controlled respiration at 10, 15, and 20 breaths·minute−1 as well as from 13 subjects at supine rest and during 60° head-up tilt. Linear markers were computed for comparison. We found that: (i) over simulations, CSampEn and KNNCUP exhibit different abilities in evaluating coupling strength; (ii) KNNCUP is more reliable than CSampEn when interactions occur according to a causal structure, while performances are similar in noncausal models; (iii) in healthy subjects, KNNCUP is more powerful in characterizing cardiorespiratory and cerebrovascular variability interactions than CSampEn and linear markers. We recommend KNNCUP for quantifying cardiorespiratory and cerebrovascular coupling.

The dynamic interactions between R and HP and between MAP and MCBv feature common characteristics: (i) they are inherently nonlinear; (ii) their strength is between the maximum and the minimum values found in the presence of full coupling and perfect uncoupling, respectively. dimensional vector obtained by concatenating y i with y − i as y i = y i y i−1 . . . y i−m+1 . We remark that y − i and y i can be interpreted indifferently as patterns of length, respectively, and m − 1 and m in the time domain or points of, respectively, the (m − 1)-and m-dimensional state phase spaces built uniformly using the method of lagged coordinates in the state phase domain. Analogously, we define x j , x − j = x j−1 . . . x j−m+1 , and x j = x j x j−1 . . . x j−m+1 the equivalent quantities computed over x. In order to characterize the relationship, if present, between X and Y, we indicate the probability that y i and x j are close in the m-dimensional state phase space within a tolerance r with p y i − x j ≤ r and with p y − i − x − j ≤ r the probability that y − i and x − j are closer than r in the (m − 1)dimensional state phase space, where · is a metric to compute distance. In this study the adopted metric is the Euclidean norm.

CSampEn
CSampEn [40] is defined as the negative logarithm of the ratio of the averaged where · performs the time average over all the reference vectors built from x and log(·) is the natural logarithm. p y i − x j ≤ r and p y − i − x − j ≤ r are estimated by counting the number of y i closer than r to x j for i = m, . . . , N and the number of y − i closer than r to x − j for i = m − 1, . . . , N − 1 and by dividing them by Nm + 1. The unfortunate case was that p y i − x j ≤ r = 0 and/or p y − i − x − j ≤ r = 0 were never observed in our study given the adopted values of m and N. CSampEn is a measure of the negative logarithm of the conditional probability that two points remain close in the m-dimensional phase space given that they are close in the (m − 1)-dimensional phase space: the closer to 1 the conditional probability, the lower the uncertainty associated with the reciprocal position of patterns in the m-dimensional phase space, the smaller the CSampEn, the stronger the relationship between X and Y.

KNNCUP
Cross-unpredictability (CUP) searches for the dependency f (·) of the future value y i+τ of the dynamics of the target Y on m − 1 past samples x − i of the dynamics of the driver X [41,47], where τ is the prediction horizon [14]. y i+τ is usually labelled as the image of x − i through f (·) and x − i is the reference vector for the search of its k nearest neighbors. Local CUP approach approximates f (·) in a region around x − i set by its k nearest neighbors [48]. Given the adopted metric for evaluating distances in the phase space, the region is a hypersphere and x − i is its center. The k nearest neighbors of the reference vector x − i were utilized to predict y i+τ . The prediction of y i+τ , namelyŷ i+τ , is defined as the weighted mean of the images of the k nearest neighbors of x − i , where the weights are the inverse of their distance from the reference vector. Vectors at zero distance from x − i were excluded from the set of k nearest neighbors.ŷ i+τ was computed as a function of time i, thus providing the predicted seriesŷ. The degree of unpredictability of Y given X is measured via the CUP function defined as 1 − ρ 2 , where ρ 2 is the squared correlation coefficient between y andŷ. CUP is bounded between 0 and 1, where 0 indicates that y can be perfectly predicted from x and 1 indicates complete unpredictability of y based on x. CUP depends on m. The course of CUP with m was the result of two opposite tendencies [49]: (i) prediction of y given x might improve and CUP might decrease because longer patterns built over x bring more information about the future behavior of y; (ii) at high m, the k nearest neighbors tend to be far away from the reference vector due to the spreading of the vectors in the phase space and this effect raises CUP. The minimum of the CUP function over m was taken as a measure of the uncoupling between x and y [41] and referred to as the CUP index (CUPI): the closer to 0 the CUPI, the greater the ability to predict y based on x, the stronger the relationship between X and Y.

Graded Unidirectional and Bidirectional Causal Couplings
We simulated bivariate autoregressive (BAR) processes whose components feature different degrees of coupling and causal relationships. The BAR process is defined as where W 1 and W 2 are Gaussian white noises with zero mean and variances assigned such that X and Y exhibit unit variance. If c 1 = 0 and c 2 = 0, X and Y are two uncoupled second-order autoregressive [AR (2)] processes. If c 1 = 0 and c 2 = 0, the directionality of the interactions is from X to Y (i.e., unidirectional causal model). If c 1 = 0 and c 2 = 0, the directionality of the interactions is from X to Y and vice versa (i.e., bidirectional causal model). In the uncoupled condition, the two processes X to Y were set to feature dominant rhythms according to two configurations: (i) ρ 1 = ρ 2 = 0.8 with phases ϕ 1 = ±3π/5 and ϕ 2 = ±3π/5 corresponding to an oscillation at normalized frequency f 1 = 0.3 cycles·sample −1 and f 2 = 0.3 cycles·sample −1 , thus simulating a dominant rhythm in the high frequency (HF) band typical of spontaneous respiration driving an HP rhythm at the same frequency with a mean HP equal to 1 s; (ii) ρ 1 = ρ 2 = 0.8 with phases ϕ 1 = ±π/5 and ϕ 2 = ±π/5 corresponding to an oscillation at normalized frequency f 1 = 0.1 cycles·sample −1 and f 2 = 0.1 cycles·sample −1 , thus leading to a dominant rhythm in the low frequency (LF) band typical of slow breathing driving an HP rhythm at the same frequency with a mean HP equal to 1 s. Coupling strengths between X and Y were varied according to the following setup: (i) c 1 = 0 and c 2 was varied incrementally from 0 to 1.0 in 0.1 steps, thus simulating unidirectional coupling from X to Y with incremental coupling strength (i.e., from the respiratory system to the heart or from the systemic to cerebral vasculature); c 1 = c 2 = c, and c was varied gradually from 0 to 1.0 in 0.1 steps, thus simulating bidirectional coupling from X to Y and vice versa with incremental coupling strength (i.e., from the respiratory system to the heart, or from the systemic to cerebral vasculature, and vice versa).

Graded Lag-Zero Noncausal Coupling
We simulated a bivariate process whose components featured different degrees of coupling in absence of any causal relationship. This bivariate process is obtained by corrupting an AR(2) process X with a Gaussian white noise W 2 with zero mean and standard deviation taken as a fraction of the standard deviation of X. The bivariate process is defined as where X is an AR(2) process with zero mean and unit standard deviation and W 2 is a Gaussian white noise with zero mean and standard deviation 1 − c 2 . As in Section 3.1, X exhibits dominant rhythms according to two pole configurations: (i) ρ 1 = ρ 2 = 0.8 with phases ϕ 1 = ±3π/5 and ϕ 2 = ±3π/5; (ii) ρ 1 = ρ 2 = 0.8 with phases ϕ 1 = ±π/5 and ϕ 2 = ±π/5. Since interactions between X and Y occur at lag-zero, directionality is not set (i.e., the coupling is instantaneous), thus simulating noncausal interactions between X and Y. If c 2 = 1, X and Y are coincident and, assigned the outcome of the random experiment, the points [x(n), y(n)] lie exactly on the diagonal line for n = 1, . . . , N, and the two processes X and Y are fully coupled. c 2 was varied incrementally from 0 to 1.0 in 0.1 steps, thus progressively decreasing the standard deviation of W 2 and increasing the coupling strength between X and Y.

Unidirectionally-Coupled Identical Logistic Maps
Nonlinear dynamics were simulated via logistic maps. We considered two unidirectionallycoupled identical logistic maps [50], described as with Thus, X and Y are deterministic signals obtained by iterating (4) starting from the initial conditions randomly chosen within the interval from 0.0 to 1.0. When c 2 = 0, X and Y evolve independently of each other because their initial conditions are different according to the equations of the two logistic maps. Given that r = 3.7, the two logistic maps are in chaotic regime. The coupling is unidirectional because Y does not affect X. The strength of the interactions from X to Y increases with c 2 . The parameter c 2 was varied from 0 to 1.0 in 0.1 steps.

Ethical Statement
The R and HP series were extracted from a historical database built to evaluate changes in the cardiorespiratory coupling with the breathing rate in healthy subjects [42,43]. The MAP and MCBv series belonged to another historical database built to study the dCA during postural stimulus in healthy subjects [44,45]. All the original protocols were in keeping with the Declaration of Helsinki. The protocols were approved by the local ethical review board of the L. Sacco Hospital, Milan, Italy, and Sacro Cuore Don Calabria Hospital, Negrar, Italy, respectively, and authorized by the same structures. Written signed informed consent was obtained from all subjects. Physical examination and full neurological evaluation certified the healthy status of all the subjects. The subjects were not under pharmacological treatments interfering with cardiovascular and cerebrovascular controls. Subjects avoided caffeinated and alcoholic beverages and heavy physical exercises for 24 h before the study.

CB Protocol
Data were acquired from 19 healthy subjects (age: 27-35 years, median = 31 years; 8 males, 11 females) at rest in a supine position during spontaneous breathing (SB) and during CB at 10, 15, and 20 breaths·minute −1 , labelled as CB10, CB15, and CB20, respectively. The period of SB always preceded the session of CB. The respiratory frequency of the CB was selected randomly. The subjects performed all the CB sessions. The timing of inspiratory and expiratory onsets was provided via a metronome and reinforced verbally by the experimenter. All the experimental sessions lasted 10 min. The subjects were not allowed to talk during the entire protocol. The electrocardiogram (ECG) was acquired from lead II via a bioamplifier (Marazza, Monza, Italy) and respiratory flow via a nasal thermistor (Marazza, Monza, Italy). Both signals were digitalized synchronously at 300 Hz by an analog-to-digital 12-bit board (National Instruments, Austin, TX, USA) plugged into a personal computer. From the ECG, we derived the beat-to-beat variability series of the HP and a downsampled version of the R signal. After identifying the R-wave from the ECG, the inter-heartbeat interval between the nth and (n + 1)th R-wave peaks was taken as the nth HP. The R signal was sampled at the nth R-wave peak.

HUT Protocol
Data were acquired from 13 healthy subjects with no history of postural syncope (age: 27 ± 8 years; 5 males, 8 females). Subjects were instrumented to continuously monitor the ECG from lead II and noninvasive continuous arterial pressure (AP) from the middle finger of the nondominant arm (Finapres Medical Systems, Enschede, The Netherlands). The cerebral blood velocity (CBv), measured from the right, or left, middle cerebral artery through a transcranial Doppler device (Multi-Dop T, DWL, 2 MHz, Compumedics, San Juan Capistrano, CA, USA), was taken as a surrogate of cerebral blood flow [51]. The CBv signal was low-pass filtered with a sixth-order Butterworth filter with a cut-off frequency of 10 Hz. The signals were acquired synchronously at a sampling rate of 1000 Hz. The subjects underwent 10 min of recording at rest in a supine position (REST) followed by HUT with a tilt table inclination of 60 • . Both sessions were under SB. Analyses were carried out 3 min after the HUT onset and within the first 10 min of HUT. None of the subjects exhibited presyncope signs. We utilized the detection of the R-wave on the ECG to trigger the process of identification of systolic AP (SAP) and diastolic AP (DAP). The nth SAP was taken as the maximum of the AP signal within the nth HP. The (n − 1)th DAP was defined as the minimum of the AP signal preceding the nth SAP. The nth MAP was computed as the ratio of the definite integral of AP between the occurrences of the (n − 1)th and nth DAP to the inter-diastolic interval. The same procedure was applied to CBv to calculate the MCBv. The fiducial points for the computation of MAP were utilized for the computation of the MCBv.

Time Domain Analysis
The HP, SAP, DAP, MAP, and MCBv series were manually checked and corrected in case of missing beats or misdetections. The effects of ectopic beats or isolated arrhythmic events were mitigated via linear interpolation. Synchronous sequences were randomly selected within the whole recordings. The length of the series was kept constant regardless of protocol and experimental condition. As to the CB protocol, we monitored the mean and variance of the HP series. These markers were labeled as µ HP and σ 2 HP , respectively, and expressed in ms and ms 2 . In the CB protocol, µ HP did not vary compared to SB, being 1010 ± 168, 989 ± 157, 1023 ± 162, and 1028 ± 162 during SB, CB10, CB15, and CB20, respectively, while σ 2 HP increased solely during CB10 compared to SB, being 3368 ± 2622, 4784 ± 3356, 3705 ± 3091, and 2813 ± 2158 during SB, CB10, CB15, and CB20, respectively. As to the HUT series, we monitored the mean and variance of the MAP and MCBv. These indices were denoted as µ MAP , σ 2 MAP , µ MCBv , and σ 2 MCBv and expressed in mmHg, mmHg 2 , cm·s −1 and cm 2 ·s −2 . In the HUT protocol, µ MCBv decreased during HUT compared to REST, being 72 ± 23 and 62 ± 21, respectively, while σ 2 MCBv increased from 19 ± 12 to 26 ± 16. In the HUT protocol, µ MAP and σ 2 MAP did not vary with HUT: µ MAP was 99 ± 17 and 95 ± 12 at REST and during HUT, respectively, while σ 2 MAP was 18 ± 21 and 19 ± 12.

Computation of a Linear Marker of Association between Time Series
Squared coherence function K 2 x,y ( f ) provides an estimation of the degree of linear association between x and y as a function of the frequency f [4]. It is computed as the ratio of the square cross-spectrum modulus between x and y divided by the product of their power spectra. By definition, K 2 x,y ( f ) is bounded between 0 and 1, where 0 and 1 indicate the minimum and maximum correlation between x and y. A linear marker of the strength of the cardiorespiratory coupling is commonly computed by sampling K 2 R,HP ( f ) at its peak in the HF band (i.e., from 0.15 to 0.4 Hz) [4,7]. This marker is indicated as K 2 R,HP (HF) in the following. A linear marker of the strength of the cerebrovascular coupling is routinely computed by sampling K 2 MAP,MCBv ( f ) at its peak in the very LF (VLF, from 0.02 to 0.07 Hz), LF (from 0.07 to 0.15 Hz), and HF (i.e., from 0.15 to 0.4 Hz) bands [8,9]. These markers are indicated as K 2 MAP,MCBv (VLF), K 2 MAP,MCBv (LF), and K 2 MAP,MCBv (HF) in the following. The superior limit of the LF band and the inferior limit of the HF were modified compared to the original definition [8] to account for possible slow respiratory rhythms [12]. K 2 R,HP ( f ) and K 2 MAP,MCBv ( f ) were estimated according to a bivariate AR model [7]. The model order was fixed to 10, and the coefficients of the bivariate AR model were identified via the least squares approach [7].

Computation of CSampEn and KNNCUP
Since the aim of the study was to explore the physiological mechanisms responsible for short-term control of HP and MCBv, the length N of the series was set to 256 [8,52]. The series were first linearly detrended and then normalized to have zero mean and unit variance. Cardiorespiratory coupling was assessed with x = R and y = HP, while the cerebrovascular link was evaluated using x = MAP and y = MCBv. CSampEn was computed over the normalized series with m = 3 and r = 0.2 [14]. KNNCUP was performed with k = 30 [41]. Over simulated data, KNNCUP was computed with time horizon τ = −1 in the case of lag-zero noncausal model and with time horizon τ = 0 in the case of unidirectional and bidirectional causal models. Over experimental data, KNNCUP was carried out with τ = −1 in agreement with the fast vagal actions responsible for the respiratory sinus arrhythmia [23,53] and fast resistive component of the MCBv-MAP relationship [54,55].

Statistical Analysis
One-way repeated measures analysis of variance (Dunnett's test for multiple comparisons), or Friedman repeated measures analysis of variance on ranks if appropriate (Dunnett's test for multiple comparisons), was utilized to check the effect of CB versus SB. A paired t-test, or a Wilcoxon signed rank test when appropriate, was applied to check the effect of HUT. Statistical analysis was carried out using a commercial statistical program (Sigmaplot, v.14.0, Systat Software, Inc., Chicago, IL, USA). A p < 0.05 was always considered statistically significant.

Results on Simulations
The line plot of Figure   Like Figure 1, Figure 2 shows the results of simulations generated via linear unidirectional causal (Figure 2a   Like Figure 1, Figure 2 shows the results of simulations generated via linear unidirectional causal (Figure 2a (Figure 3b) computed over dynamics generated via unidirectionally-coupled identical logistic maps while varying c 2 . The curves were built over 20 pairs of signals generated according to different initial conditions. CUPI decreased to 0 with c 2 and reached 0 when X and Y synchronized (Figure 3b). Conversely, CSampEn remained stable with c 2 , and this result was the consequence of the inability of CSampEn to detect the situation of uncoupling when c 2 was close to 0 (Figure 3a). while varying c 2 . The curves were built over 20 pairs of signals generated according to different initial conditions. CUPI decreased to 0 with c 2 and reached 0 when and synchronized (Figure 3b). Conversely, CSampEn remained stable with c 2 , and this resul was the consequence of the inability of CSampEn to detect the situation of uncoupling when c 2 was close to 0 (Figure 3a).  Figure 4 shows CSampEn ( Figure 4a) and CUPI (Figure 4b) computed in the CB pro tocol. CSampEn decreased significantly during CB10 compared to SB, thus suggesting tha cardiorespiratory coupling increased at the slowest breathing rate. CUPI exhibited a sim ilar trend, but its variation compared to SB was significant during both CB10 and CB15 thus making more evident the dependence of the cardiorespiratory coupling strength on the respiratory rate during CB.     (Figure 4b) computed in the CB protocol. CSampEn decreased significantly during CB10 compared to SB, thus suggesting that cardiorespiratory coupling increased at the slowest breathing rate. CUPI exhibited a similar trend, but its variation compared to SB was significant during both CB10 and CB15, thus making more evident the dependence of the cardiorespiratory coupling strength on the respiratory rate during CB. while varying c 2 . The curves were built over 20 pairs of signals generated according to different initial conditions. CUPI decreased to 0 with c 2 and reached 0 when and synchronized (Figure 3b). Conversely, CSampEn remained stable with c 2 , and this result was the consequence of the inability of CSampEn to detect the situation of uncoupling when c 2 was close to 0 (Figure 3a).  Figure 4 shows CSampEn ( Figure 4a) and CUPI (Figure 4b) computed in the CB protocol. CSampEn decreased significantly during CB10 compared to SB, thus suggesting that cardiorespiratory coupling increased at the slowest breathing rate. CUPI exhibited a similar trend, but its variation compared to SB was significant during both CB10 and CB15, thus making more evident the dependence of the cardiorespiratory coupling strength on the respiratory rate during CB.       . The vertical box-and-whisker plots show the K 2 marker computed between R and HP in the HF band in the CB protocol (a) and the K 2 markers computed between MAP and MCBv in the VLF (b), LF (c), and HF (d) bands in the HUT protocol. The height of the box represents the distance between the first and third quartiles, with the median marked as a horizontal segment, and the whiskers denote the 5th and 95th percentiles. The symbol § indicates p < 0.05 versus CB or REST.

Discussion
The main findings of this study can be summarized as follows: (i) CSampEn and KNNCUP can quantify the modifications of coupling strength but with different abilities; (ii) KNNCUP is more reliable than CSampEn when interactions occur according to a causal structure, while performances are similar over a lag-zero linear noncausal model; (iii) in healthy subjects KNNCUP is more powerful in identifying the changes in cardiorespiratory and cerebrovascular coupling strength in response to experimental challenges than CSampEn and linear markers.     . The vertical box-and-whisker plots show the K 2 marker computed between R and HP in the HF band in the CB protocol (a) and the K 2 markers computed between MAP and MCBv in the VLF (b), LF (c), and HF (d) bands in the HUT protocol. The height of the box represents the distance between the first and third quartiles, with the median marked as a horizontal segment, and the whiskers denote the 5th and 95th percentiles. The symbol § indicates p < 0.05 versus CB or REST.

Discussion
The main findings of this study can be summarized as follows: (i) CSampEn and KNNCUP can quantify the modifications of coupling strength but with different abilities; (ii) KNNCUP is more reliable than CSampEn when interactions occur according to a causal structure, while performances are similar over a lag-zero linear noncausal model; (iii) in healthy subjects KNNCUP is more powerful in identifying the changes in cardiorespiratory and cerebrovascular coupling strength in response to experimental challenges than CSampEn and linear markers. Figure 6. The vertical box-and-whisker plots show the K 2 marker computed between R and HP in the HF band in the CB protocol (a) and the K 2 markers computed between MAP and MCBv in the VLF (b), LF (c), and HF (d) bands in the HUT protocol. The height of the box represents the distance between the first and third quartiles, with the median marked as a horizontal segment, and the whiskers denote the 5th and 95th percentiles. The symbol § indicates p < 0.05 versus CB or REST.

Discussion
The main findings of this study can be summarized as follows: (i) CSampEn and KNNCUP can quantify the modifications of coupling strength but with different abilities; (ii) KNNCUP is more reliable than CSampEn when interactions occur according to a causal structure, while performances are similar over a lag-zero linear noncausal model; (iii) in healthy subjects KNNCUP is more powerful in identifying the changes in cardiorespiratory and cerebrovascular coupling strength in response to experimental challenges than CSampEn and linear markers.

Assessing the Coupling Strength between Dynamic Systems via CSampEn and KNNCUP
The degree of interaction between dynamic systems X and Y is assessed in the present study via two model-free methods, namely CSampEn and KNNCUP.
CSampEn reconstructs the dynamics of X and Y in two separate embedding spaces built using the technique of delayed components and assesses the strength of the relationship between the embedding spaces by computing the negative logarithm of the conditional probability that two patterns generated, respectively, by X and Y that are close in the (m − 1)-dimensional embedding space remain close in the m-dimensional one. In other words, according to the philosophy of state space correspondence methods, if two patterns occupy the same region in the (m − 1)-dimensional embedding spaces, in presence of a significant coupling between X and Y, they do not move apart when a new component is added.
Therefore, although the model-free nature of CSampEn allows for the theoretical description of nonlinear interactions, strong nonlinearities that imply links among different regions of the phase space, such as those responsible for interactions among rhythms at different frequencies, cannot be reliably described and the method seems to be more suitable for describing relations occurring according to a 1:1 coupling ratio.
However, simulations provided in this study prove that, even in the presence of a 1:1 coupling ratio typically occurring when two processes with the same dominant oscillation are considered, CSampEn cannot detect the progressive modification of the coupling strength when interactions between X and Y are generated by a linear causal model. This conclusion held even in the case of a model simulating the interactions between two nonlinear deterministic signals generated by logistic maps in a chaotic regime. This limitation is related to the inability of CSampEn to interpret causality given that it is a symmetric metric under the reversal of the role between X and Y. The most regrettable feature of CSampEn is that it can suggest even the opposite trend with the coupling strength as shown in Figures 1c and 2a,c. Solely in trivial simulations where two stochastic processes interact in absence of a causal structure (i.e., immediate interactions), CSampEn can detect the expected modifications with the coupling strength. The limited ability of CSampEn confirms that the exploitation of a metric based on conditional probability does not necessarily assure a causal approach [56].
KNNCUP assesses the relationship between X and Y using a completely different logic, namely the cross-predictability of Y from patterns taken from the past history of X. KNNCUP recalls model-based cross-conditional entropy even though the metric utilized to assess irregularity of future behaviors of Y given the activity of X is different from the evaluation of the logarithm of the variance of the innovation of Y given X [57]. Therefore, the pattern is created over the activity of X and the ability of this pattern to set future behaviors of Y is tested by measuring the ability of predicting y based on the knowledge of x. In addition to having the possibility to describe nonlinear interactions among rhythms at different frequencies according to the number of past samples of x utilized to predict y, KNNCUP could account for the causal structure of the mechanisms generating the interactions between X and Y. In addition to accounting for the directionality of the interactions, KNNCUP optimizes the embedding dimension and coarse graining via the search for the minimum over m [41,49] and the coverage of the embedding space with cells of different size according to the k-nearest-neighbor strategy [58], respectively. The most remarkable feature of KNNCUP is its robustness in indicating the expected modifications of the coupling strength regardless of the model structure (i.e., causal or noncausal), type of causal interactions (i.e., unidirectional or bidirectional), and signal feature (i.e., linear stochastic or nonlinear deterministic).
In conclusion, KNNCUP should be preferred to CSampEn as a model-free technique for the assessment of coupling strength even in the presence of interactions occurring at the same frequency.

Superior Ability of CUPI Compared to CSampEn in Evaluating Cardiorespiratory Coupling Strength
CSampEn and CUPI exhibited similar results in the CB protocol. Indeed, both methods detected an increase in cardiorespiratory coupling strength during CB10. This result is not surprising given that the same conclusion was reached via cross-conditional entropy based on a uniform quantization strategy for the computation of probabilities [42]. A linear approach based on K 2 R,HP ( f ) was able to reach the same conclusion, but the trend with the breathing rate was not evident. Therefore, given the dependence of the cardiorespiratory coupling strength on the breathing rate, we recommend the use of a respiratory ratematched control group in applications to pathological subjects. At first sight, this conclusion might be the consequence of the increase of the transfer function magnitude from R to HP while slowing the respiratory rate [4,21,22] resulting from the shape of the sinus node transfer function [59,60]. Conversely, since a rise in the transfer function magnitude is not necessarily accompanied by an increase of the coupling strength, the observed finding is more likely to be the consequence of a firmer pacing of the neural efferent activity operated by the respiratory centers [27][28][29][30]. Given that the influence of R on HP is very rapid [4,23,53] and produces dominant HP oscillations at the same frequency [22,23], the difficulties of CSampEn in assessing coupling strength are less evident in the CB protocol. However, CSampEn and CUPI are not fully equivalent. As a matter of fact, CSampEn exhibited a slightly lower statistical power in detecting the dependence of the cardiorespiratory coupling strength on the respiratory rate than CUPI. This finding suggests that some difficulties of SampEn might be linked to the presence of slower components in the responses of HP to R [5,61] and to the significant action in the reverse causal direction (i.e., from HP to R) [3,26,53]. In addition, since the likelihood of the cross-frequency coupling, for example, between the main respiratory rhythm and oscillations generated by baroreflex control loop [62,63] increases while slowing the breathing rate, the resulting impairment of the state space correspondence might have, more importantly, reduced the effectiveness of CSampEn compared to that of CUPI.

Superior Ability of CUPI Compared to CSampEn in Evaluating Cerebrovascular Coupling Strength
CUPI declined during HUT, thus indicating that postural challenge augmented the cerebrovascular coupling strength. This result might be the consequence of the raise of AP variability [64,65] that is not buffered with suitable changes to the vessel diameter and drives changes of MCBv via the pressure-to-flow relationship [66,67]. This consideration suggests a tendency toward a worsening of the dCA given that an increased association between MAP and MCBv variations has been observed during fainting [68] and in subjects with impaired dCA [9,[69][70][71]. A linear approach based on K 2 MAP,MCBv ( f ) suggested an increased cerebrovascular coupling strength as well. However, it is worth noting that the observed increase is likely to occur in a physiological range that does not compromise the dCA given that the quality of the MCBv response to a sustained step increase of MAP was found to be preserved during HUT [72,73]. Regardless of the magnitude of the changes, the robust assessment of the MCBv-MAP coupling strength is a critical issue. Indeed, traditional markers of MCBv-MAP association based on K 2 MAP,MCBv ( f ) failed to detect an increased cerebrovascular coupling strength in pathological subjects with impaired dCA [74].
Remarkably, CSampEn did not decrease during HUT, thus stressing its limited ability to quantify the modifications of the cerebrovascular coupling strength during orthostatic challenge. This result might be the consequence of the causal structure of the dynamic MCBv-MAP interactions occurring in a closed loop [45,54,55,75] according to the pressureto-flow link [76] and the Cushing-like pathway [77] mediated by the activity of the mechanoreceptors of the brainstem [78] and fully operative under physiological modifications of intracranial pressure [79].

Conclusions
This study proves the different abilities of CSampEn and KNNCUP metrics in assessing the coupling strength between two time series and clarifies the conditions under which their behaviors diverge. In addition to the trivial disruption of the state space correspondence resulting from nonlinear interactions occurring according to the n:m coupling ratio, with m and n different from 1, even interactions occurring according to a causal structure are sufficient to limit the performance of CSampEn. We recommend the avoidance of the use of CSampEn to assess cardiorespiratory and cerebrovascular coupling from spontaneous variability, because CSampEn might be inadequate to describe the complex causal nature of these relationships [4][5][6]42,44,45,55,62,75,80]. Conversely, KNNCUP assures a more reliable nonlinear model-free framework to quantify cardiorespiratory and cerebrovascular variability interactions because its performances are less dependent on the structure of the underlying mechanism generating the dynamics, thus allowing the detection of a trend toward an increase of the cardiorespiratory coupling strength while slowing the breathing rate and the raise of the cerebrovascular coupling strength during the sympathetic activation and vagal withdrawal induced by an orthostatic challenge.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author (i.e., A.P.) upon permission of the IRCCS Policlinico San Donato. The data are not publicly available because they contain information that could compromise the privacy of research participants.