Reliable detection of directional couplings using cross-vector measures

Consider two stationary dynamical systems $X$ and $Y$ that exhibit a self-sustained motion in the absence of a coupling between them. Assume that if there is a coupling, it is unidirectional from $X$ to $Y$⁠. This means that the state of the driver $X$ influences the motion of the response $Y$⁠, but $Y$ has no influence on $X$⁠. The aim of this work is to further improve techniques that allow for the characterization of such driver-response dynamics from the analysis of pairs of time series. Here, one time series is measured from the driver, and the other is measured from the response. Such techniques must be sensitive to and specific for couplings, i.e., they should detect a coupling if and only if it is actually present. In particular, it is essential to detect the correct coupling direction for a broad range of non-zero coupling strengths. We, therefore, use the following notions. Given a coupling $X→Y$⁠, for a technique to be direction-sensitive, it should detect a coupling in this direction. Moreover, to be direction-specific, the technique should not falsely detect a coupling $Y→X$⁠.

Čenys et al.¹⁵ gave the key to achieving these goals based on state-space reconstructions from the time series of $X$ and $Y$⁠. They pointed out that if the states that the response dynamics $Y$ attains at two times are similar, the states of the driver dynamics $X$ at these times are likely also to be similar. The strength of this conditioning increases with increased coupling strength. The conditioning in the other direction, i.e., similarity of states in $X$ implying similarity of states in $Y$⁠, holds to a significantly lesser degree. This principle is utilized by numerous state-space-based approaches that aim at characterizing the driver-response dynamics.^15–24 We here follow Ref. 25 and refer to it as the asymmetric state similarity conditioning principle.

Importantly, the way in which the different state-space-based approaches quantify the asymmetric state similarity conditioning can have a decisive impact on their direction-sensitivity and direction-specificity (see, for example, Refs. 20, 22, and 24). Aiming to optimize this performance, we recently introduced the approach of cross-distance vectors and defined a novel measure, denoted here by $c d$⁠, for the detection of directional couplings from the prominence of the initial tails of these vectors.²⁴ We compared this cross-distance vector measure $c d$ against a state-space-based approach that was introduced and referred to as the nonlinear interdependence measure $M$ in Ref. 21. This comparison provided evidence for the advantages of $c d$ over $M$⁠, including a better direction-specificity.²⁴ However, for both $M$ and $c d$⁠, the similarity between states is quantified by distances between the corresponding points in reconstructed state-spaces of $X$ and $Y$⁠. This use of distances can lead to bias caused by asymmetries between $X$ and $Y$ that are unrelated to couplings, such as different effective dimensions of these individual dynamics.²² To diminish this bias, Ref. 22 introduced the nonlinear interdependence measure $L$⁠, which uses ranks of distances instead of the values of the distances. A further very prominent state-space-based approach to characterize coupled dynamics is the convergent cross mapping measure $ρ$⁠,²³ in which distances are used to construct weights in cross prediction estimates between $X$ and $Y$⁠. Therefore, the following question arises. Can the performance of the cross-vector approach be improved by using vectors of ranks or vectors of cross prediction estimates instead of vectors of distances? To address this question, we here introduce the cross-rank vector measure $c r$ and the cross-estimate vector measure $c e$⁠. We compare these novel measures against the cross-distance vector measure $c d$⁠, the nonlinear interdependence measures $M$ and $L$⁠, and the convergent cross mapping measure $ρ$⁠.

In Sec. II, we start by reviewing the definitions of the previously established measures, $M$⁠,^21, $L$⁠,^22, $ρ$⁠,²³ and $c d$⁠.²⁴ We then introduce the definitions of $c r$ and $c e$⁠. In Sec. III A, we apply all six measures to noise-free unidirectionally coupled nonlinear deterministic dynamics. Since this application shows the advantages of the rank-based measures $L$ and $c r$ over the other measures, we restrict the remaining analysis to these two measures. In Sec. III B, we compare the noise robustness of $L$ and $c r$⁠. Section III C explores the performance of these measures on short time series. Section III D provides an application to stochastic dynamics and illustrates the effect of linear cross correlation. To compare the measures’ performance on real-world data, an application to electroencephalographic (EEG) recordings from epilepsy patients²⁶ is given in Sec. IV. For this part of our analysis, we combine $L$ and $c r$ with the method of surrogates.^21,27,28

Now, we introduce the central concept of the $k$th conditioned nearest neighbor of $x i$⁠, defined as $x P i Y ( k )$⁠. If there is coupling from $X$ to $Y$⁠, due to the asymmetric state similarity conditioning principle, the conditioned nearest neighbors of $x i$ are commonly also its true neighbors in the sense that they are close in space to $x i$⁠. In contrast, in the absence of this coupling, the conditioned neighbors of $x i$ are typically not close in space to $x i$⁠. Keep in mind that we have the pruned distance matrices of size $N 1× N 2$ that take into account the Theiler correction and our additional constraint. Therefore, for all $N 1$ points, we have the same number of $N 2$ candidates for both the true nearest neighbors and conditioned nearest neighbors.

All six state-space-based measures used in this work share the initial steps described up to this point. Based on these definitions, they quantify the closeness of the conditioned nearest neighbors in different ways. However, the four measures $M$⁠, $L$⁠, $c d$⁠, and $c r$ also share the subsequent steps. The permutations $P i X$ and $P i Y$ are used to construct the following ordered matrices. At first, the permutation $P i X$ is applied to the $i$th row of $D X$ across $i=1,…, N 1$ resulting in the matrix $D X → X$⁠, where $D i j X → X= D i , P i X ( j ) X$ is the distance from $x i$ to its $j$th nearest neighbor. Then, in analogy, the permutation $P i Y$ is applied to the $i$th row of $D X$ across $i=1,…, N 1$ resulting in the cross-distance matrix $D X → Y$⁠, where $D i j X → Y= D i , P i Y ( j ) X$ is the distance from $x i$ to its $j$th conditioned nearest neighbor. Analogous procedures are used to determine $D Y → Y$ and $D Y → X$⁠.

All three measures defined thus far reach the maximal value of 1 when the two time series are identical. In the other limiting case of independent dynamics, their values are distributed around 0.

Reference 24 reported a substantially better direction-specificity of $c d$ compared to the distance-based nonlinear interdependence measure $M$⁠. In this paper, we expand the idea of quantifying the vector tails to ranks [analogous to $L$⁠, Eq. (3)] and to time series estimates [analogous to $ρ$⁠, Eq. (6)]. In what follows, we accordingly introduce two new cross-vectors, namely, the cross-rank vector and the cross-estimate vector.

Let us briefly address the differences between the three variants of the cross-vectors. Despite that the ordinates in the different panels of Fig. 1 have different orders of magnitude, all three variants of vectors result in very similar shapes. However, the relative length of the initial tail and, therefore, the values of the obtained measures $c d , X → Y, c r , X → Y$⁠, and $c e , X → Y$ differ substantially. This difference will be thoroughly discussed in Sec. III A.

In this section, we compare the direction-sensitivity and direction-specificity of the nonlinear interdependence measures $M$ and $L$⁠, the convergent cross mapping measure $ρ$⁠, and the three cross-vector approaches. We first analyze noise-free deterministic dynamics in Sec. III A, then study the impact of noise in Sec. III B, explore the performance on short time series in Sec. III C, and, finally, consider stochastic dynamics and the effects of linear cross correlation in Sec III D.

We begin by analyzing the dependence of the six measures on the coupling strength between noise-free deterministic dynamics. Recall that each measure has two complementary variants: one to detect couplings from $X$ to $Y$ and the other to detect couplings from $Y$ to $X$⁠. We analyze three pairs of deterministic chaotic dynamical systems with a mismatch in their parameters, unidirectionally coupled in the direction $X→Y$ with diffusive coupling, with strength $ε$⁠:

1. Rössler systems:

   $x ˙ 1 = − 1.01 x 2 − x 3 , x ˙ 2 = 1.01 x 1 + 0.15 x 2 , x ˙ 3 = 0.2 + x 3 ( x 1 − 10 ) , y ˙ 1 = − 0.99 y 2 − y 3 + ε ( x 1 − y 1 ) , y ˙ 2 = 0.99 y 1 + 0.15 y 2 , y ˙ 3 = 0.2 + y 3 ( y 1 − 10 ) .$

   (13)
2. Lorenz systems:

   $x ˙ 1 = 10 ( x 2 − x 1 ) , x ˙ 2 = 39 x 1 − x 2 − x 1 x 3 , x ˙ 3 = x 1 x 2 − 8 3 x 3 , y ˙ 1 = 10 ( y 2 − y 1 ) + ε ( x 1 − y 1 ) , y ˙ 2 = 35 y 1 − y 2 − y 1 y 3 , y ˙ 3 = y 1 y 2 − 8 3 y 3 .$

   (14)
3. Duffing systems:

   $x ¨ + 0.3 x ˙ − x + x 3 = 0.5 cos ⁡ ( 1.2 t ) , y ¨ + 0.2 y ˙ − y + y 3 = 0.3 cos ⁡ ( t ) + ε ( x − y ) .$

   (15)

The exclusion window is $W=50$ samples, the $M,L$⁠, and $ρ$ measures’ nearest neighbor parameter from Eqs. (2)–(5) is $k=10$⁠, and the cross-vector measures’ parameters from Eqs. (8), (10), and (12) are $k 1=10 and k 2=100$⁠. We choose $k= k 1$⁠, which allows for an equitable comparison across all six measures, as $k$ and $k 1$ both determine the number of nearest neighbors. Following the previous work,²⁴ we set $k 2=10⋅ k 1$⁠. This previous work showed that setting a too high value of $k 2$ generally decreases the direction-specificity of $c d$⁠, and setting $k 2$ to the maximal possible value of $k 2= N 2$ makes $c d$ perform almost identically to $M$⁠. Accordingly, an adequate choice of the parameter $k 2$ is key for good direction-specificity. For simplicity, we use the embedding parameters $τ=1 andm=20$⁠. As long as $m$ is sufficiently large, the product $(m−1)τ$ matters more than both individual parameter values. This is because this product determines the temporal distance between the first and the last entry in the embedding vectors [Eq. (1)]. All parameters are kept the same for the analysis of all three deterministic dynamics.

For each system, we vary the coupling strength $ε$ equidistantly on a logarithmic scale, from a negligibly small value up to and above the threshold that leads to generalized synchronization³³ in different dynamics. We, furthermore, include the uncoupled case, i.e., $ε=0$⁠. For each value of $ε$⁠, we compute the distribution of the measures’ values across $100$ realizations of the dynamics. The six measures are always calculated for the same set of realizations so that the randomness of the initial conditions has no impact on the comparison across measures.

The coupling direction in our systems is $X→Y$⁠. Accordingly, to be direction-sensitive, a measure should detect a coupling $X→Y$⁠. Moreover, to be direction-specific, the measure should not detect a coupling $Y→X$⁠. We assess direction-sensitivity by the minimal coupling strength value $ε SE$⁠, above which the distribution of values of the measure’s variant used to detect a coupling $X→Y$ differs significantly from the distribution obtained for the uncoupled case. A measure is judged as more direction-sensitive if it has a smaller $ε SE$⁠. Likewise, we assess direction-specificity by the minimal value of $ε SP$⁠, above which the distribution of values of the measure’s variant used to detect a coupling $Y→X$ differs significantly from the distribution at $ε=0$⁠. Since the actual coupling direction is $X→Y$⁠, a measure is more direction-specific if it has a larger $ε SP$⁠. We use the Kolmogorov–Smirnov test at a significance level $α=0.01$ to assess whether two distributions differ. We determine the values of $ε SE$ and $ε SP$ as the first of at least five consecutive values of $ε$ that result in rejections of the Kolmogorov–Smirnov test. The results of the analysis are shown in Fig. 2.

For all measures, the variants used to detect couplings $X→Y$ increase with increasing coupling strength and are larger than the complementary variants used to detect couplings $Y→X$ across wide ranges of coupling strength. This is to be expected of any directional coupling measure. The values $ε SE$ are comparable among five of the six measures. The exception is $ρ$ with slightly smaller values of $ε SE$⁠, which indicates that it is the most direction-sensitive measure. However, it also has the largest variance of the distributions of its values across different realizations, as indicated by broad colored bands in Figs. 2(g)–2(i). This can undermine a better direction-sensitivity in practice as, often, only a small number of realizations of time series pairs is available, and a higher variability of $ρ$ can make the estimation of this measure’s mean value less reliable. The other five measures have not only similar direction-sensitivity but also similar variance across realizations, making them comparable in this aspect.

The value $ε SP$ of any of the measures $c d$⁠, $c r$⁠, or $c e$ is larger than the value $ε SP$ of any of the measures $M$⁠, $L$⁠, or $ρ$⁠. In most cases, $ε SP$ is around ten times larger for the cross-vector measures than for the respective classical counterparts. This means that the $c$ measures are substantially more direction-specific and can be used to reliably infer the coupling direction for a broader range of coupling strengths. The most direction-specific out of all six measures are $c r$ and $c d$⁠. Next is $c e$⁠, and then $ρ$⁠, and the least direction-specific measures are $M$ and $L$⁠.

It is well-known that the coupling direction is generally not well-defined for synchronized dynamics since we no longer have two separate dynamics. Indeed, the values $ε SP$ of measures $c d$ and $c r$ are closely below the generalized synchronization threshold. When the coupling strength is further increased above the threshold, all measures lose the ability to infer the true coupling direction.

Recall that we use coupled non-identical chaotic systems that enter into generalized synchronization. Although they cannot show complete synchronization, the resulting two time series become almost identical at very large $ε$⁠, and four out of six measures reach values close to 1 in all three dynamics. Conversely, the measure $c d$ remains considerably below 1 in all the dynamics, and $c e$ does so for the Duffing dynamics. Furthermore, the values reached by $c d$ and $c e$ differ between the three dynamics, which necessitates further interpretation.

Consider computing the measures $c d$⁠, $c r$⁠, and $c e$ from two identical time series $x( t s)$ and $y( t s)$⁠. The values of the components of $v X → Y$ cannot be smaller than the smallest entry of the distance matrix $D X$⁠. Additionally, a time series estimate obtained from the $k$th conditioned nearest neighbors $x ^ k$ cannot perfectly replicate the original time series $x$⁠. This is because even the closest neighbor to a point is never identical to the point itself. Consequently, the similarity of the estimate to the original time series is limited by the closeness of the nearest neighbors. This example shows that the measures $c d$ and $c e$ are influenced by the distributions of values in the distance matrices. Therefore, comparing $c d , X → Y$ and $c d , Y → X$ (or $c e , X → Y$ and $c e , Y → X$⁠) to determine the coupling direction is less equitable when the distributions of values in $D X$ and $D Y$ are different.

Conversely, the measure $c r$ does not depend on the distributions of distances. It is derived from ranks, which are uniformly distributed regardless of the pair of time series being analyzed. The entries of the cross-rank vector obtained from two identical time series $x( t s)=y( t s)$ are $r j X → Y=j$⁠, resulting in $c r , X → Y= k 2 k 2 + 1$⁠. This holds for all pairs of identical time series. The measure $c r$ is, therefore, not influenced by the distributions of values in the distance matrices. Accordingly, it can determine the coupling direction more accurately than $c d$ and $c e$⁠.

In conclusion, the measures $c r$ and $c d$ have comparable direction-sensitivity and direction-specificity, while $c e$ is slightly less powerful in these aspects. Due to the additional advantage of uniform distribution of ranks, we consider $c r$ the best cross-vector measure. Therefore, we use $c r$ and its classical counterpart $L$ in the subsequent analysis.

It is important to test if a measure stays reliable when the data contain noise, as the measured data often do. To assess the noise robustness of $c r$ and $L$⁠, we add measurement noise to the three deterministic systems from Sec. III A. Specifically, we add Gaussian white noise with zero mean and increasing standard deviation $σ N$⁠. We use noise levels equally spaced on a logarithmic scale starting from $σ N=0.1⋅ σ S$⁠, where $σ S$ is the standard deviation of the noise-free time series. For each noise level, we determine the values $ε SE$ and $ε SP$ of the measures $L$ and $c r$⁠. The results are shown in Fig. 3, which includes the results for the noise-free dynamics already shown in Fig. 2 for comparison.

As the noise level increases, the value $ε SE$ of both measures generally increases for all three systems. The values $ε SE$ of $L$ and $c r$ are nearly identical for noise-free data. However, with increasing noise, $ε SE$ increases faster for $c r$ than for $L$ for all three systems. This shows that $L$ retains better direction-sensitivity than $c r$ when noise is introduced. On the other hand, when the noise level is large, the values $ε SE$ become very similar to $ε SP$⁠. This happens at larger $σ N$ for $c r$ than for $L$⁠, showing that $c r$ retains better direction-specificity. Accordingly, $c r$ allows for correct inference of the coupling direction across a substantially wider range of coupling strength values at all noise levels. Thus, $c r$ exhibits better overall noise robustness than $L$⁠.

In real-world applications, the measured time series are often short. This section, therefore, presents the analysis of the three dynamics from Sec. III A for different time series lengths $N 0$⁠. The time series length is varied linearly on a logarithmic scale. The minimal value is $N 0=200$⁠, which corresponds to 6–10 oscillations in these different dynamics. As the maximal value, we use $N 0= 10 4$⁠, which is the time series length taken in Secs. III A and III B. The embedding parameters $m=20$ and $τ=1$⁠, as well as the Theiler correction window $W=50$⁠, remain unchanged. Note that at the shortest $N 0=200$⁠, this leaves only $N 2= N 0−(m−1)τ−(2W+1)=80$ candidates for the nearest neighbors in total. Accordingly, we deliberately include a very short time series in this analysis, for which hardly any detection of the coupling can be expected.

For each time series length, we determine the values $ε SE$ and $ε SP$ of the measures $L$ and $c r$⁠. Here, care must be taken in choosing the parameters $k 1$ and $k 2$ from Eq. (10), and $k$ from Eq. (3). Recall that these parameters determine the number of nearest neighbors. If a time series is half as long, a chosen point will, on average, have half as many neighbors in a neighborhood of a certain size. Using the same number of nearest neighbors would result in using more distant neighbors. Accordingly, if the time series are half as long, the parameters $k$⁠, $k 1$⁠, and $k 2$ should be half as large in order for the results to be comparable. Therefore, we adjust the parameters $k= k 1$ and $k 2$ such that the ratios $k 1 N 0= 10 − 3$ and $k 2 N 0= 10 − 2$ are kept constant. However, $k 1$ cannot be smaller than $1$⁠, and $k 2$ cannot be smaller than $10$⁠, which ensures the ratio $k 1 k 2=10$ is kept constant for all $N 0$⁠.

The results are shown in Figs. 4(a)–4(c). First, we observe that the values of $ε SE$ of both $L$ and $c r$ decrease with increasing $N 0$⁠. This is expected, as generally longer time series lead to an increase in detection power. The values of $ε SE$ are generally comparable for both measures, indicating that the measures’ direction-sensitivity remains similar between them in short time series. The exception to this are coupled Rössler systems, where $L$ has better direction-sensitivity for very short time series. Conversely, $c r$ continues to be more direction-specific than $L$ for all $N 0$ and all three dynamics, as its values of $ε SP$ remain substantially larger than those of $L$⁠. An interesting observation is that $ε SP$ of $L$ decreases with increasing $N 0$⁠. At first sight, it seems contradictory that the direction-specificity of this measure deteriorates with increasing time series length. It can, however, be understood straightforwardly by the following considerations. Like we have outlined above, the inverse of the asymmetric state similarity conditioning principle (closeness in the driver $X$ maps to closeness in the response $Y$⁠) holds to a much weaker degree than the principle itself (closeness in the response $Y$ maps to closeness in the driver $X$⁠). Although it is weaker, the inverse conditioning still holds and is quantified by $L(Y |X)$⁠. The quantity $ε SP$ is the coupling at which this inverse conditioning is detected by $L(Y |X)$⁠. Accordingly, the deterioration of $L$’s direction-specificity with increasing time series length is explained by an increased sensitivity of $L(Y |X)$ for the inverse conditioning, causing the incorrect detection of a coupling from the response to the driver. Conversely, the values $ε SP$ of $c r$ remain close to the generalized synchronization threshold for all $N 0$⁠. This shows that $c r$ is much less affected by this inverse conditioning, further explaining its better direction-specificity.

To illustrate the importance of adequately setting the nearest neighbor parameters, Figs. 4(d)–4(f) show the results obtained for fixed values $k= k 1=10$ and $k 2=100$⁠.³⁶ We can see that direction-sensitivity is comparable to the results for the adjusted parameters seen in Figs. 4(a)–4(c). However, the direction-specificity values $ε SP$ of $c r$ are much smaller when the parameters are not adjusted to the changing $N 0$⁠. This shows the negative effect of a poor choice of parameters on direction-specificity.

We have shown that the measures $L$ and $c r$ allow one to detect directional couplings in nonlinear deterministic dynamics. However, it is well-known that coupling detection measures can also be influenced by other time series properties. In particular, a strong linear cross-correlation between the time series $x( t s)$ and $y( t s)$ can result in high values of the measures. The problem here is that even in the absence of coupling, one can obtain an arbitrary degree of linear cross-correlation, for example, by mixing the time series originating from independent dynamics. Furthermore, it is also important to assess how these measures behave in stochastic systems, as randomness can potentially affect the reliability of the measures. To assess the effects of linear cross-correlation and randomness on our measures, we analyze linear and nonlinear stochastic processes.

For the two systems (17) and (18), we vary $A$ and $ε$⁠, respectively, on a linear scale and generate $100$ realizations for each value. One realization results in a pair of time series, each consisting of $10 4$ samples. As embedding parameters, we use $τ=1$ and $m=10$⁠. The exclusion window is $W=50$⁠, the $L$ measure’s parameter from Eq. (3) is $k=10$⁠, and the $c r$ measure’s parameters from Eq. (10) are $k 1=10, k 2=100$ for all cases.

For system (17), we evaluate the direction-specificity by determining the minimal $A>0$ at which the results differ significantly from those at $A=0$⁠. We use the same statistics and criteria as for determining $ε SP$⁠, resulting in the analogous quantity $A SP$⁠. We can, however, not determine any quantity analogous to $ε SE$ since system (17) has no coupling and we cannot evaluate the measures’ direction-sensitivity.

Figure 5(a) shows the results of the analysis of system (17). The two time series $x 1( t s)$ and $y 1( t s)$ are always generated by a mixing of independent dynamics (16). In the limit of $A=0$⁠, the time series $x 1( t s)$ and $y 1( t s)$ remain independent. In consequence, the mean values of the directional coupling measures are very close to zero. As $A$ increases and with it the cross-correlation, both $L$ and $c r$ also increase. This illustrates the previously mentioned impact of the linear cross-correlation on coupling detection measures. This can be regarded as a limited specificity of these measures since arbitrary degrees of cross-correlation can be obtained by mixing signals from uncoupled dynamics. At complete mixing, i.e., $A=0.5$⁠, the time series $x 1( t s)$ and $y 1( t s)$ are identical, and accordingly, both measures reach their respective maximal values. Figure 5(a) also shows that $c r$ is less affected by linear cross-correlation than $L$⁠. First, the value of $A SP$ is bigger for $c r$ than for $L$⁠. This shows that for $c r$⁠, a stronger cross-correlation is needed to cause a significant difference between results with and without a cross-correlation. The second aspect is mean values of the measures obtained across multiple realizations of the dynamics for a given value of $A$⁠. Upon increasing $A$⁠, the curve of $c r$ stays close to zero longer and increases slower than the one of $L$⁠. This is of particular relevance for the application to real-world dynamics where often there is no controllable parameter like $A$⁠, and the mean value of the measures across realizations for a fixed setting is all that is available.

Figure 5(b) shows the results of the analysis of system (18). Recall that the coupling direction is $X→Y$⁠. The mean values of $L(X |Y)$ and $c r , X → Y$ both increase with increasing $ε$ and are substantially larger than the mean values of $L(Y |X)$ and $c r , Y → X$⁠, respectively. This shows that both $L$ and $c r$ can be used to determine the coupling direction in stochastic systems. Comparing the values of $ε SE$ reveals that $L$ is more direction-sensitive in this scenario. However, $c r$ is once again found to be more direction-specific than $L$⁠, as indicated by the mean values of $L(Y |X)$ increasing above zero, while the mean values of $c r , Y → X$ remain very close to zero. Furthermore, the value $ε SP$ of $L$ is only slightly larger than its $ε SE$⁠, which means that $L$ detects the correct direction for only a narrow range of coupling strengths. In contrast, the values $c r , Y → X$ obtained for $ε>0$ are significantly different from those at $ε=0$ only at two non-consecutive coupling strength values, so that $ε SP$ remains undefined. Accordingly, also in this setting, $c r$ exhibits strong direction-specificity.

Finally, an important observation can be made by analyzing the cross-correlation values in Fig. 5(b). The mean cross-correlation between $x 2( t s)$ and $y 2( t s)$ stays close to zero regardless of the value of $ε$⁠. However, $L(Y |X)$ still increases with increased $ε$⁠, even though the coupling direction is $X→Y$⁠. This shows that strong couplings from $X$ to $Y$ can cause an increase of the measure $L(Y |X)$⁠, which is supposed to detect a coupling from $Y$ to $X$⁠, even without a linear cross-correlation. This lack of direction-specificity in the absence of linear cross-correlation is not found for $c r$⁠.

In this section, we apply the rank-based measures $c r$ and $L$ to electroencephalographic (EEG) recordings from epilepsy patients to assess the measures’ applicability to real-world data. We use the publicly available Bern–Barcelona EEG database.²⁶ In the initial analysis of this database, the $L$ measure was combined with bivariate surrogates.²⁶ We here extend this analysis by including not only $L$ but also $c r$⁠. We, therefore, at first, briefly describe the database in Sec. IV A, outline the concept of surrogates in Sec. IV B, and describe details of our EEG analysis and results in Sec. IV C.

The database consists of 7500 pairs of signals extracted from intracranial EEG recordings from five epilepsy patients. Each signal pair was recorded simultaneously at the neighboring EEG channels and has a length of $20$ s at a sampling frequency of $512$ Hz. Half of the pairs were recorded from brain areas where the first seizure-related EEG signal changes were detected for individual patients (focal signal pairs). The other half were recorded from brain areas that were not involved in patients’ seizures at their onset (nonfocal signal pairs). Focal and nonfocal signal pairs were picked randomly across patients, time, and channel pairs within the respective brain areas. The location of the seizure onset zone determines whether a given signal pair is focal or nonfocal. However, all signals were taken from recordings of the seizure-free time interval. In other words, neither the 3750 focal nor the 3750 nonfocal signal pairs include actual seizure activity. An exemplary pair of focal signals is shown in Fig. 6(a). For a more detailed description of the signal acquisition, composition, and preprocessing, the reader is referred to the paper.²⁶ 

In Sec. III, we use model systems under controlled conditions. This allows us to directly obtain null distributions of the measures expected for uncoupled and independent dynamics. For the deterministic dynamics and the nonlinear stochastic dynamics, we generate uncoupled dynamics by setting $ε=0$ in their respective Eqs. (13), (14), (15), and (18). Analogously, for the linear mixing of independent autoregressive processes, we generate the non-mixed system by setting $A=0$ in Eq. (17). Our analysis in Sec. III is then based on statistical comparisons of the results obtained for these null distributions against results obtained for nonzero couplings (⁠ $ε>0$⁠) and nonzero mixing (⁠ $A>0$⁠).

When dealing with experimental data, it is typically not possible to generate signals for uncoupled and independent dynamics. This problem can be addressed by the concept of surrogate data.^21,27,28 Surrogates are a very powerful and flexible tool in signal analysis. They allow estimating the null distribution of measures expected under well-defined null hypotheses about the dynamics underlying some measured data. For this purpose, a constrained randomization is applied to the original data, resulting in one realization of the surrogate data. In particular, the randomization constraints are defined such that the surrogate data preserve all properties of the original data that are consistent with the null hypothesis but are otherwise random.

In our case, the data are given by a pair of time series. We construct pairs of surrogate time series that have the amplitude distribution and auto-correlation of the individual time series, as well as the cross-correlation between the time series, in common with the original pair of time series. All other properties of the time series, including possible signatures of the asymmetric state similarity conditioning between nonlinear dynamics, are destroyed by randomization. This corresponds to the null hypothesis $H 0$ that “the dynamics is a stationary bivariate linear stochastic correlated Gaussian process. The measurement functions by which the pair of time series was derived from the dynamics are invertible but potentially nonlinear. The auto-correlation, cross-correlation, mean, and variance of the underlying Gaussian process are such that the measurement results in the auto-correlation, cross-correlation, and amplitude distribution of the observed time series.”²⁶ 

To generate the surrogates, we use the bivariate iterative amplitude adjusted Fourier transform algorithm described in Ref. 28. In line with Ref. 28, our choice to use this particular surrogate algorithm is motivated by a previous study of EEG recordings from epilepsy patients.³⁷ Apart from the surrogates used here, Ref. 37 tested time-shifted surrogates, which represent the null hypothesis that the dynamics are independent and not cross-correlated. No assumption is made about the nature of individual dynamics. Accordingly, time-shifted surrogates maintain the individual signals’ structure and autocorrelation but not the cross-correlation between signals. Reference 37 showed that a combination of $L$ with the surrogates used in the present study allows one for a much better discrimination between focal signals and nonfocal signals than a combination of $L$ with time shifted surrogates. It is, therefore, essential to not only preserve the autocorrelation but also the cross correlation of the signals in our study of the pairs of EEG signals. An exemplary pair of surrogate time series constructed for the EEG signals shown in Fig. 6(a) is displayed in Fig. 6(b).

To test $H 0$⁠, we calculate the measures $c r$ and $L$ for the original pair of time series and a set of 19 independent realizations of pairs of surrogate time series. We reject the null hypothesis if a measure’s value obtained for the original pair of time series is larger than for any of its surrogates, corresponding to a one-sided rank-order test with significance level $α=0.05$⁠.²⁸ 

Unlike the analysis in Ref. 26, we did not downsample the data by a factor of four. We, therefore, took a four times larger time embedding window and autocorrelation exclusion window $W$⁠. Thus, we use embedding parameters $τ=1$⁠, $m=32⋅4=128$⁠, and $W=19⋅4=76$⁠. We use $τ=1$ as we did in Sec. III since it is the temporal distance between the first and the last entry in the embedding vectors $(m−1)τ$ that matters. The algorithm parameters are $k= k 1=10 and k 2=100$⁠. Figures 6(c) and 6(d) show the cross-rank vectors obtained for the exemplary pair of EEG signals [Fig. 6(a)] along with the vectors obtained for $19$ surrogate signals pairs [one of which is shown in Fig. 6(b)]. Figures 6(e)–6(h) show the corresponding $c r$ and $L$ indices.

Recall that the signal pairs of the Bern–Barcelona EEG database were selected randomly from different patients, times, and brain areas. Moreover, the assignment to $X$ and $Y$ was random for each pair. Hence, the directions $X→Y$ and $Y→X$ carry no particular meaning for this real-world data. Therefore, the authors of Ref. 26 did not consider the aspect of direction and used a symmetrized measure to quantify the overall strength of the interdependence: $L S=(L(X |Y)+L(Y |X))/2$⁠. In order to compare our results to the ones of Ref. 26, we adopt this approach and define $c r , S= ( c r , X → Y+ c r , Y → X )/2$⁠.

To study the influence of non-stationarity on the proportions of rejections of $L$⁠, Ref. 26 introduced a surrogate-based test for the temporal variability of the signals’ amplitude distribution, power spectrum, and cross-correlation between signals. The test was designed to be very strict, as a rejection with any one of the three features was considered a rejection of the stationarity test. Using this test, $2000$ and $1750$ focal signal pairs were classified as stationary and non-stationary, respectively. For the nonfocal signal pairs, these numbers were found to be $1674$ and $2076$⁠, respectively. These classification results are publicly available along with the EEG data.³⁸ The last three columns of Table I show the proportions of rejections we obtained exclusively for the subsets of stationary signals along with the resulting contrasts $D$⁠.

For $L S$⁠, the proportions of rejections decrease substantially when excluding the non-stationary signal pairs. This holds especially for nonfocal signals. Consequently, the contrast of $L S$ increases from $D=0.160$ to $D=0.273$⁠. Conversely, for $c r , S$⁠, only the proportion of rejections for the nonfocal signals increases slightly, resulting in a small decrease from $D=0.313$ to $D=0.297$⁠. This shows that $c r , S$ is much less sensitive to non-stationarity. Accordingly, with or without the exclusion of non-stationary signal pairs, the contrast between focal and nonfocal signals is higher for $c r , S$ than for $L S$⁠.

In this paper, we considered the problem of inferring the strength and the direction of couplings from measured bivariate time series. This work is a continuation of Ref. 24, in which we introduced a new state-space-based directional coupling measure $c d$⁠, which is derived from the tails of the cross-distance vectors $v$⁠. Expanding on this concept, we here introduced two novel measures $c r$ and $c e$⁠. They are obtained by quantifying the tails of the cross-rank vectors $r$ and the cross-estimate vectors $e$⁠, respectively.

We first studied the direction-sensitivity and the direction-specificity of the cross-vector measures $c d$⁠, $c r$⁠, and $c e$⁠, the nonlinear interdependence measures $M$²¹ and $L$⁠,²² and the convergent cross mapping measure $ρ$⁠.²³ Our analysis of simulated coupled dynamics shows that all three cross-vector measures allow for considerably more reliable detection of the coupling direction for a wider range of dynamics compared to the other measures. This higher direction-specificity is the main advantage of this cross-vector approach. Furthermore, while the convergent cross mapping measure $ρ$ demonstrates slightly higher direction-sensitivity than the other five measures, it also exhibits the largest variance across different realizations of the dynamics. The other five measures have comparable direction-sensitivity. The measure $c r$ performed at least as well as $c d$ and $c e$ in all aspects. The fact that it is computed from ranks provides additional advantages. It can reach values closer to $1$ and is more robust to asymmetries between $X$ and $Y$ since the distributions of ranks are uniform for all time series. Therefore, $c r$ also allows for more equitable testing with surrogates than $c d$ or $c e$⁠, as the distributions of ranks are identical for both original and surrogate time series. This does not hold for the distributions of distances or cross-correlation values between the original and estimated time series. Hence, we consider $c r$ as the most promising cross-vector measure and, thus, studied its performance further. In particular, we compared it to the established rank-based measure $L$ to ensure a fair evaluation, given that both measures are rank-based.

We investigated the impact of measurement noise on the two chosen measures. As noise levels increase, the measure $c r$ retains better direction-specificity, while $L$ preserves better direction-sensitivity. Overall, $c r$ allows for a correct inference of the coupling direction across a substantially wider range of coupling strength than $L$ for noisy time series. This shows that $c r$ is a better choice when discerning the coupling direction from noisy data. Furthermore, we explored the performance of the rank-based measures on short time series. We found that $L$ slightly outperformed $c r$ in terms of direction-sensitivity in very short time series. However, $c r$ was once again found to be substantially more direction-specific than $L$ for all time series lengths considered. Accordingly, $c r$ is more reliable in inferring the correct coupling direction also for short time series. Additionally, we analyzed the signals generated by the linear mixing of two independent stochastic dynamics, which allowed us to introduce an arbitrary degree of linear cross-correlation to the pair of time series. Our results show that $c r$ is substantially less affected by linear cross-correlation than $L$⁠. We conjecture that this is the key factor contributing to the better direction-specificity of the cross-vector measure. However, this is not the sole factor, as further analysis reveals that $c r$ is significantly more direction-specific than $L$ even for two nonlinearly coupled stochastic dynamics with zero linear cross-correlation. This is an important finding, suggesting that the measure $c r$ is unaffected by multiple properties of the time series that are not directly related to directional coupling.

Finally, we applied the rank-based measures $c r$ and $L$ to the Bern–Barcelona EEG database,²⁶ a publicly available database of EEG signals from epilepsy patients. Previous studies using state-space-based measures have shown that focal signals exhibit a larger degree of nonlinear interdependence compared to nonfocal signals during the seizure-free interval.^18,26,37,39 In this paper, we build upon the study in Ref. 26, which used the measure $L$ combined with bivariate surrogates, by including the measure $c r$⁠. Our analysis shows that $c r$ can better distinguish between focal and nonfocal signals compared to $L$⁠, as shown by a larger contrast $D$⁠. Excluding those signals classified as non-stationary increases the contrast of $L$ close to that of $c r$⁠. Conversely, the results for $c r$ remain nearly unchanged. This is an important finding, as it reveals that $c r$ is much less influenced by non-stationarity in the analyzed time series compared to $L$⁠.

An important point needs to be made about bidirectional couplings. Recall that state-space-based measures were introduced using the asymmetric state similarity conditioning principle, which enables them to quantify the strength of unidirectional couplings in driver-response dynamics. In practice, they can also accurately quantify the influence of both directional couplings $X→Y$ and $Y→X$ when the coupling is bidirectional.²⁴ However, one should be cautious about asymmetries between the two systems $X$ and $Y$ in such scenarios, as these measures quantify the coupling impact rather than the coupling strength parameter from the governing equations.⁴⁰ 

In conclusion, the new measure $c r$ demonstrates substantially improved accuracy in inferring the coupling direction, compared to the established state-space approaches we considered. To specifically infer the direction of coupling, $c r$ has proven to be the most reliable choice across the various settings we examined. Additionally, it shows great potential for applications to real-world data, as illustrated by our study of focal and nonfocal EEG signals from patients with epilepsy. Furthermore, its robustness to non-stationarity reflected in time series eliminates the need to exclude such signals, enabling more effective use of valuable data. In closing, our study offers a new perspective on utilizing the asymmetric state similarity conditioning principle for the detection of directional couplings, which provides an interesting theoretical insight. We hope this advancement will spark further research and inspire additional investigations in this area.

M.B. and P.B. acknowledge the Research Core Funding No. P2-0001 and the Research Project No. J3-4525, both financially supported by the Slovenian Research and Innovation Agency. R.G.A. acknowledges funding from the Spanish Ministry of Science and Innovation and the State Research Agency (Grant No. PID2020-118196GBI00/MICIU/AEI/10.13039/501100011033).

The authors have no conflicts to disclose.

Martin Brešar: Conceptualization (supporting); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (equal); Visualization (lead); Writing – original draft (equal). Ralph G. Andrzejak: Conceptualization (lead); Data curation (lead); Methodology (supporting); Supervision (equal); Validation (equal); Writing – original draft (equal). Pavle Boškoski: Investigation (supporting); Resources (lead); Software (supporting); Supervision (equal); Writing – original draft (supporting).

The simulated data that support the analysis of this article have been generated by the authors and can be fully reproduced from the repository in Ref. 41. The EEG data analyzed in this article are the publicly available Bern–Barcelona EEG database.³⁸