# Directional coupling detection through cross-distance vectors

Martin Brešar<EMAIL_ADDRESS>Jožef Stefan Institute, Jamova cesta 39,
SI-1000 Ljubljana, Slovenia  Jožef Stefan International Postgraduate School,
Jamova cesta 39, SI-1000 Ljubljana, Slovenia  Pavle Boškoski  Jožef Stefan
Institute, Jamova cesta 39, SI-1000 Ljubljana, Slovenia

###### Abstract

Inferring the coupling direction from measured time series of complex systems
is challenging. We propose a new state space based causality measure obtained
from cross-distance vectors for quantifying interaction strength. It is a
model-free noise-robust approach that requires only a few parameters. The
approach is applicable to bivariate time series and is resilient to artefacts
and missing values. The result is two coupling indices that quantify coupling
strength in each direction more accurately than the already established state
space measures. We test the proposed method on different dynamical systems and
analyse numerical stability. As a result, a procedure for optimal parameter
selection is proposed, circumventing the challenge of determining the optimal
embedding parameters. We show it is robust to noise and reliable in shorter
time series. Moreover, we show that it can detect cardiorespiratory
interaction in measured data. A numerically efficient implementation is
available at https://repo.ijs.si/mbresar/cd-vec.

##  I Introduction

Complex systems found in nature can often be considered as many interacting
subsystems. Subsystems are often inherently connected and can not be
considered isolated from each other, which raises the question of how they
interact with each other. In this article, we propose a new method that
belongs to the family of state space distance approaches. It is capable of
detecting and quantifying interactions in a computationally efficient way. It
can be applied to bivariate time series to quantify the coupling strength in
both directions. Furthermore, it is applicable to both linear and nonlinear
coupling.

When considering two subsystems, there are four possibilities for the
direction of their interaction. They can be independent, unidirectionally
coupled (in either direction), or bidirectionally coupled. Another property
often of interest is the nature of coupling, e.g., linear or nonlinear.
Additionally, the coupling can be time-dependent, which poses a new problem in
detecting it. The problem becomes even more complex in the case of more than
two subsystems.

Typically, we measure a subsystem. For example, consider the system of a human
body. We can characterise the cardiac subsystem by measuring the heart’s
electrical activity with an electrocardiogram. By performing such measurements
for each subsystem, we obtain time series. The goal is to infer the direction
and the nature of the interactions between the subsystems from measured time
series.

This article shows that the proposed coupling indices can infer the coupling
direction for various regular and chaotic systems. It can quantify the
coupling strength in both directions for unidirectional or bidirectional
coupling. It turns out to be noise robust and also efficient for shorter time
series. Most importantly, compared to the most prominent state space
approaches, the proposed indices more accurately and reliably detect the
coupling direction. We discuss the possible pitfalls of state space approaches
in detail. We present an in-depth analysis of the behaviour of the proposed
quantities. As a result, we propose a procedure of selecting optimal parameter
values, thus achieving a robust performance. In such a way, we circumvent the
problem of optimal embedding parameters selection. Additionally, we present a
means of dealing with artefacts or missing values. Finally, we test the
effectiveness of the approach on measured cardiorespiratory data.

### Related work

Identifying causal relationships arises in different fields that deal with
complex systems. For this purpose, different methods for detecting coupling
between subsystems are being developed. Some of the most widely used are
Granger causality, information theory, phase dynamics, and state space
methods. Prominent examples of areas where they are applicable are physiology
[1, 2], neuroscience [3], earth system sciences [4], ecology, [5, 6] and
economics [7, 8].

A brief overview of coupling detection methods is given below. For a more
comprehensive overview of the available approaches, the prospective reader is
referred to the overviews [9, 10] and references therein.

#### Granger causality

Coupling detection methods often follow the idea of Granger causality [7]. The
method in the original work is based on fitting a vector autoregressive model.
Based on this idea, numerous other methods were proposed, such as Partial
directed coherence [11], which is a multivariate frequency approach to the
original Granger causality.

#### Information transfer

A widely used family of methods for studying interactions is based on
information theory. Most commonly studied is the transfer entropy [12], which
quantifies information transfer from one subsystem to another. It is
mathematically equivalent to conditional mutual information [13]. It is a
model-free method that can detect both linear and nonlinear coupling. Entropy
methods rely on estimating multidimensional probability distributions, which
can be challenging for shorter time series.

#### Phase dynamics

For oscillatory systems, methods based on the oscillation phase can be used.
Rosenblum and Pikovsky [14] proposed a method for quantifying asymmetry in the
interaction between two oscillating subsystems. Phase transformation can be
used together with information measures to more accurately quantify
information flow between oscillators [15]. Time-frequency approaches can also
be used, such as wavelet phase coherence [16] for detecting linear and wavelet
bispectral analysis [17] for detecting nonlinear interactions.

#### State space methods

In weakly coupled bivariate systems, close states of the driven subsystem are
mapped to close states of the driving subsystem. The opposite effect is much
smaller. Different methods quantifying this effect in both directions have
been proposed [18, 19, 20].

#### Bivariate versus multivariate

An important distinction is between bivariate and multivariate methods. In the
case of more than two interacting subsystems, all of them must be accounted
for when analysing interactions. This necessity led to the multivariate
approach of many of the mentioned approaches. Partial granger causality [21]
is an extension to the Granger causality that excludes the effects of latent
variables. Causation entropy [22] is a measure similar to transfer entropy
that considers multiple variables. Multivariate methods, however, often lead
to the curse of dimensionality, which means the estimation of multivariate
measures becomes increasingly problematic with an increasing dimension of the
measured system.

##  II Cross-distance vectors

Consider a pair of unidirectionally coupled subsystems
$\bm{x}(t)=(x_{1}(t),\dots,x_{n_{x}}(t))$ and
$\bm{y}(t)=(y_{1}(t),\dots,y_{n_{y}}(t))$. If the coupling direction is
$y\rightarrow x$, their time evolution is described by

$\displaystyle\frac{d\bm{x}(t)}{dt}=\bm{f}(\bm{x},t)+\bm{g}(\bm{x},\bm{y})$
(1a) $\displaystyle\frac{d\bm{y}(t)}{dt}=\bm{h}(\bm{y},t),$ (1b)

where $\bm{g}(\bm{x},\bm{y})$ is the coupling function. We consider time-
independent coupling functions, though they generally can depend on time. Our
goal is to define a measure for the magnitude and direction of coupling from
the observed time series. For $\bm{g}(\bm{x},\bm{y})\neq 0$, the subsystems
and, therefore, their trajectories are not independent.

The trajectories of both subsystems $\bm{x}(t)$ and $\bm{y}(t)$ are observed
at equally spaced times $t_{i}=t_{0}+i\Delta t$ for $i\in\mathbb{Z}$. Thus, we
obtain time series of these trajectories $(\bm{x}(t_{i}):i\in\mathbb{Z})$ and
$(\bm{y}(t_{i}):i\in\mathbb{Z})$. We assume that $\Delta t$ is sufficiently
small to capture all the necessary information. The influence of $\Delta t$
and potential downsampling are described in detail in Section IV.5.
Furthermore, we assume that the signals do not contain any trends, no
commensurate frequency components, and other trivial artefacts that can be
easily removed by a simple preprocessing.

Most often, only one dimension of a subsystem is measured. Thus, a single time
series of length $N$ is obtained. Therefore, we consider time series of one-
dimensional values and omit the bold notation.

We split the time series into segments of length $L$. There is a total of
$N-L+1$ segments for each time series. Segments at a time moment $t_{i}$ are
defined as

$\displaystyle\bm{o}^{x}_{i}=(x(t_{i}),x(t_{i+1}),\dots,x(t_{i+L-1}))$ (2)
$\displaystyle\bm{o}^{y}_{i}=(y(t_{i}),y(t_{i+1}),\dots,y(t_{i+L-1})).$ (3)

Furthermore, let us assume a measure of the similarity of two segments, i.e.,
a distance between two vectors

$d(\bm{o}^{x}_{i},\bm{o}^{x}_{j})=\|\bm{o}^{x}_{i}-\bm{o}^{x}_{j}\|.$ (4)

We ask two questions:

  1. 1.

If two segments of the driven subsystem $\bm{o}^{x}_{i}$ and $\bm{o}^{x}_{j}$
at times $t_{i}$ and $t_{j}$ are similar, are the segments at those times of
the driving subsystem $\bm{o}^{y}_{i}$ and $\bm{o}^{y}_{j}$ also similar?

  2. 2.

And inversely: If two segments of the driving subsystem $\bm{o}^{y}_{i}$ and
$\bm{o}^{y}_{j}$ at times $t_{i}$ and $t_{j}$ are similar, are the segments at
those times of the driven subsystem $\bm{o}^{x}_{i}$ and $\bm{o}^{x}_{j}$ also
similar?

To clarify, we ask whether the following statements are true in the coupling
direction $y\rightarrow x$:

$1.\>\>\>\>\>\>\>\>d(\bm{o}_{i}^{x},\bm{o}_{j}^{x})\approx 0\implies
d(\bm{o}_{i}^{y},\bm{o}_{j}^{y})\approx 0$ (5)
$2.\>\>\>\>\>\>\>\>d(\bm{o}_{i}^{y},\bm{o}_{j}^{y})\approx 0\implies
d(\bm{o}_{i}^{x},\bm{o}_{j}^{x})\approx 0$ (6)

It is well known that statement (5) is true and (6) is false for weakly
coupled oscillators [18]. We provide an explanation of why this is the case
for different subsystems in Appendix A. State space methods take advantage of
this property to infer the coupling direction from measured time series. In
what follows, we propose a new measure that quantifies coupling strength in
each direction more accurately than the already established state space
indices.

###  II.1 Cross-distance vectors algorithm

First, we construct segments (2) and (3) and choose a distance measure (4). In
this article, we use

$d(\bm{o}_{i}^{x},\bm{o}_{j}^{x})=\sqrt{\frac{1}{L}\sum_{m=1}^{L}\big{(}x(t_{i+m-1})-x(t_{j+m-1})\big{)}^{2}}.$
(7)

Next, we construct the distance matrix

$\bm{D}^{x},\>D^{x}_{ij}=d(\bm{o}^{x}_{i},\bm{o}^{x}_{j})$ (8)

of size $(N-L+1)\times(N-L+1)$, which contains the distances between all pairs
of segments. We similarly calculate the distance matrix of the other time
series

$\bm{D}^{y},\>D^{y}_{ij}=d(\bm{o}^{y}_{i},\bm{o}^{y}_{j}).$ (9)

The goal is to rearrange the $i$-th row of $\bm{D}^{y}$ in the same order as
if the $i$-th row of $\bm{D}^{x}$ was to be sorted in ascending order. This is
achieved by calculating an index permutation $P^{x}_{i}$ so that the $i$-th
row is sorted in ascending order under that permutation. We do this for every
row and obtain $N-L+1$ permutations.

Now, we sort $i$-th row of $\bm{D}^{y}$ according to the permutation
$P^{x}_{i}$ for all rows and obtain the matrix

$\bm{D}^{y\rightarrow x},\>D^{y\rightarrow
x}_{ij}=D_{i,P_{i}^{x}(j)}^{y}\>\>.$ (10)

This matrix highlights whether segments of $y$ are similar at times at which
segments of $x$ are similar. For example, if the first few elements of $i$-th
row of matrix (10) are the smallest entries in that row, this indicates that
statement (5) is true. To evaluate this for all rows, we finally average out
the rows of the matrix to obtain the cross-distance vector

$\bm{v}^{y\rightarrow x},\>v^{y\rightarrow
x}_{j}=\frac{1}{N-L+1}\sum_{i}D^{y\rightarrow x}_{ij}.$ (11)

Similarly, we calculate the matrix $\bm{D}^{x\rightarrow y}$ and the vector
$\bm{v}^{x\rightarrow y}$ by inverting the roles of $x$ and $y$. The vector
$\bm{v}^{y\rightarrow x}$ can be used to assess the truth of statement (5) and
thus detect coupling in the direction $y\rightarrow x$. Conversely, the vector
$\bm{v}^{x\rightarrow y}$ can be used to assess the truth of statement (6) and
thus detect coupling in the direction $x\rightarrow y$. The complete algorithm
is schematically shown in Figure 1.

It should be noted that the first entry $v_{0}^{y\rightarrow x}$ of cross-
distance vectors (11) will always be precisely zero due to zero diagonals of
the distance matrices. So, $v_{1}^{y\rightarrow x}$ is the first non-zero
entry. In the subsequent analysis, we will always omit the zero values in
graphs.

$\drawmatrix[width=.5,height=1.5]{\mathbf{o}_{i}^{x}}$
$\drawmatrix[diag,bbox/.appendstyle={fill=blue!10},width=1.5,height=1.5]0$
$D_{ij}^{x}$Index Sort$\drawmatrix[width=1.5,height=1.5]{P_{i}^{x}}$
$\drawmatrix[width=.5,height=1.5]{\mathbf{o}_{i}^{y}}$
$\drawmatrix[diag,bbox/.appendstyle={fill=blue!10},width=1.5,height=1.5]0$
$D_{ij}^{y}$Sort
rows(10)$\drawmatrixset{bboxstyle={fill=blue!10}}\drawmatrix[width=.2,height=1.5,bboxwidth=1.5,,bboxheight=1.5]{D_{ij}^{y\to
x}}$ (11)$\bm{v}^{y\to x}$ $\begin{matrix}0\\\ \vdots\\\ 0\\\ \end{matrix}$ 0s
at diagonal0s at diagonal$L$$N-L+1$$N-L+1$$N-L+1$$L$$N-L+1$$N-L+1$(8)(9)
Figure 1: Schematic representation of the algorithm for calculating the cross-
distance vector $\bm{v}^{y\to x}$ in (11).

###  II.2 Coupling index

We can define two coupling indices to quantify coupling strength in each
direction. Ideally, an index is zero when there is no coupling and increases
with increased coupling. Let us again consider unidirectionally coupled
subsystems (1). Mind that the coupling direction is $y\rightarrow x$, i.e.,
$x$ is the driven and $y$ is the driving subsystem. Consider three limits of
the coupling strength: no coupling, weak coupling and strong coupling. The
general behaviour of cross-distance vectors in these limits is shown in Figure
2.

The two subsystems $\bm{x}$ and $\bm{y}$ are independent if there is no
coupling, i.e., $\bm{g}(\bm{x},\bm{y})=0$ in (1a). Therefore, one subsystem’s
sorting permutation $P_{i}^{x}$ is random for the subsystem $y$ and the cross-
distance vector (11) is expected to be roughly constant. In the limit of
infinite time series, both $\bm{v}^{y\rightarrow x}$ and $\bm{v}^{x\rightarrow
y}$ will be constant. There are examples where this is not the case (e.g.
uncoupled subsystems which contain an oscillatory component with the same
frequency), but these are exceptions.

In the limit of strong coupling, the term $\bm{g}(\bm{x},\bm{y})$ in (1a)
causes the subsystems to synchronise. In this case, the cross-distance vectors
generally cannot be used to infer the coupling direction.

If there is weak coupling in the direction $y\rightarrow x$, the values of
$\bm{v}^{y\rightarrow x}$ change significantly from a constant since statement
(5) is true. The beginning of $\bm{v}^{y\rightarrow x}$ decreases (note the
initial tails in Figure 2 (a)). Apart from this initial tail, the bulk of
$\bm{v}^{y\rightarrow x}$ remains roughly constant (though often gains an
increasing trend). Conversely, the other cross-distance vector
$\bm{v}^{x\rightarrow y}$ stays roughly constant without the initial tail
since statement (6) is false. This is shown in Figure 2 (b).

The decrease in the initial values of of $\bm{v}^{y\rightarrow x}$ is what
allows us to detect the coupling direction from the cross-distance vectors,
since the tail is present due to the term $\bm{g}(\bm{x},\bm{y})\neq 0$ in
equation (1a). Mind that if a similar term also exists in (1b) (which means
that coupling is bidirectional), the initial values of both cross-distance
vectors $\bm{v}^{y\rightarrow x}$ and $\bm{v}^{x\rightarrow y}$ exhibit an
initial tail. To obtain the coupling indices, i.e., two values quantifying the
detected coupling strength in each direction, we must somehow quantify this
effect. While this can be done in different ways, we propose a simple but
effective approach.

Figure 2: A generic example of the behaviour of cross-distance vectors
$\bm{v}^{y\rightarrow x}$ and $\bm{v}^{x\rightarrow y}$ in different limits of
the coupling strength. The direction of coupling in the example is
$y\rightarrow x$. The orange arrows represent increasing coupling strength for
the three weak coupling examples. Note the appearance of tails in the
beginning of $\bm{v}^{y\rightarrow x}$, but not in $\bm{v}^{x\rightarrow y}$.

We define the coupling index $c^{y\rightarrow x}$, which quantifies coupling
strength in the direction $y\rightarrow x$, as the normalised difference
between the means of $v^{y\rightarrow x}_{j}$ for small $j$ and for larger $j$

$c^{y\rightarrow x}=\frac{V^{y\rightarrow x}_{k_{1}+1:k_{2}}-V^{y\rightarrow
x}_{1:k_{1}}}{V^{y\rightarrow x}_{1:k_{2}}},$ (12)

where

$V^{y\rightarrow x}_{i:j}=\frac{1}{j-i+1}\sum_{k=i}^{j}v_{k}^{y\rightarrow
x}.$ (13)

The possible range of $k_{1}$ and $k_{2}$ is $1\leq k_{1}<k_{2}\leq N-L$.

To quantify coupling strength in the other direction, we calculate
$c^{x\rightarrow y}$ by swapping $x$ and $y$. Defined in this way, coupling
indices are zero in the absence of coupling and increase with increasing
coupling strength (at least in the limit of weak coupling). Due to
normalisation, indices will roughly range from $0$ to $1$. They can be
slightly negative if there is no coupling due to finite time series, as seen
in most examples in this article. This is, however, not a problem, as a
negative index can be taken as a strong indicator of the absence of coupling.

The purpose of equation (12) is to quantify the prominence of the initial
tail. The parameter $k_{1}$ defines the part of $\bm{v}^{y\rightarrow x}$ that
includes only the initial tail. Conversely, $k_{2}$ defines the subsequent
part, where the rate of increase of $\bm{v}^{y\rightarrow x}$ is significantly
smaller.

One might argue that this selection is arbitrary. Other options, such as
kurtosis, could be explored to quantify the prominence of the initial tail.
For simplicity, we quantify coupling strength with (12) in the subsequent
analysis.

#### Related state space causality measure

From the family of the state space approaches, our index (12) is closest to
$M(Y|X)$ [19]. It will therefore be compared to and used as a benchmark to
test the reliability of the proposed indices $c$. Using our notation, $M(Y|X)$
is defined as

$M(Y|X)=\frac{1}{N-L+1}\sum_{i=1}^{N-L+1}\frac{R_{i}(Y)-R_{i}^{k}(Y|X)}{R_{i}(Y)-R_{i}^{k}(Y)},$
(14)

where

$\displaystyle R_{i}(Y)=\frac{1}{N-L}\sum_{j=1,j\neq i}^{N-L+1}D_{ij}^{y},$
(15a) $\displaystyle
R_{i}^{k}(Y)=\frac{1}{k}\sum_{j=2}^{k+1}D_{ij}^{y\rightarrow y},$ (15b)
$\displaystyle R_{i}^{k}(Y|X)=\frac{1}{k}\sum_{j=2}^{k+1}D_{ij}^{y\rightarrow
x}.$ (15c)

Roughly speaking, $M(Y|X)$ quantifies the mean of the initial tail of
$\bm{v}^{y\rightarrow x}$ compared to the mean of the whole vector, which is
highly influenced by possible trends in the vector. This undesired property is
also present in other state space measures. On the other hand, the indices
(12) quantify only the prominence of the initial tail, ignoring possible
trends and thus resulting in smaller values in the direction without coupling,
which is desired.

###  II.3 An example

Consider the system of analytically solvable unidirectionally coupled harmonic
oscillators.

$\displaystyle\ddot{x}=-\omega_{1}^{2}x+\epsilon(y-x)$ (16)
$\displaystyle\ddot{y}=-\omega_{2}^{2}y$ (17)

The coupling direction is $y\rightarrow x$. The solution to these equations
are oscillations with frequencies $\omega_{2}$ and
$\sqrt{\omega_{1}^{2}+\epsilon}$. For example, given the initial conditions
$x(0)=1,\dot{x}(0)=0,y(0)=1,\dot{y}(0)=0$, the solution is

$\begin{split}x(t)&=\frac{1}{\epsilon+\omega_{1}^{2}-\omega_{2}^{2}}\big{(}\epsilon\cos(\omega_{2}t)+(\omega_{1}^{2}-\omega_{2}^{2})\cos(\sqrt{\epsilon+\omega_{1}^{2}}t)\big{)}\\\
y(t)&=\cos(\omega_{2}t).\end{split}$ (18)

Since we have the analytical solution, we can verify statements (5) and (6)
directly.

Two segments of $x$ will be similar when the phases of both its oscillatory
components will match, one of them being $\cos(\omega_{2}t)$. Since this is
also the component of $y(t)$, its segments at those times will also be
similar. This confirms that the statement (5) is indeed true.

If two segments of $y(t)$ are similar, the phase of only one oscillatory
component of $x(t)$ will match. In contrast, the phase of the other component
can take any value (unless the frequencies of the components of $x(t)$ are
commensurable, which is an exception). This confirms that the statement (6) is
indeed false.

Numerically calculated cross-distance vectors are shown in Figure3. As
expected, the beginning of $\bm{v}^{y\rightarrow x}$ is significantly smaller
than the bulk (Figure3 (a)), and $\bm{v}^{x\rightarrow y}$ stays roughly
constant at nonzero coupling (Figure3 (b)).

Additionally, Figure3 shows that the coupling index of the direction of
coupling $c^{y\rightarrow x}$ increases with increased coupling parameter
$\epsilon$. In contrast, the coupling index in the direction without coupling
$c^{x\rightarrow y}$ remains roughly constant. By simply comparing the values,
we can correctly determine the coupling direction.

We also notice that the end of $\bm{v}^{y\rightarrow x}$ has larger values
than the bulk of the vector, which indicates that a statement similar to (5)
is also true for dissimilarity: if two segments at times $t_{i}$ and $t_{j}$
of the driven subsystem are very dissimilar, the segments at these times of
the driving subsystem are more likely also to be dissimilar. This opposite
effect is, however, not as expressed and does not appear in all systems.

When numerically calculating the cross-distance vectors, we added Gaussian
noise $\mathcal{N}(0,\sigma^{2}=10^{-4})$ to the time series. In noiseless
periodic systems, cross-distance vectors can exhibit periodic oscillations. By
adding small noise, they disappear, and the vectors behave the same as in
chaotic systems. We still used this example due to the availability of the
analytical solution. In practice, periodic systems are of lesser interest and
always contain noise.

In the following section, we will present a numerical analysis that shows that
cross-distance vectors can be used to determine the direction of coupling in
different systems.

Figure 3: The cross-distance vectors for the system of two coupled harmonic
oscillators (18). The system parameters are $\omega_{1}=0.83,\omega_{2}=2.11$,
the time series parameters are $N=10^{4},\Delta t=0.05$, the segment length is
$L=10$, and the coupling index parameters are $k_{1}=100,k_{2}=1000$.

##  III Different systems analysis

The previous section shows that the cross-distance vectors can detect the
coupling direction in the coupled harmonic oscillators system. Here, we will
show that detection is possible not only for regular subsystems but also for
discrete chaotic subsystems and for autonomous and nonautonomous continuous
chaotic subsystems. For an overview on chaotic dynamical systems, the reader
is referred to [23, 24]. We chose systems that have already been analysed with
a coupling detection method. The indices $M(Y|X)$ and $M(X|Y)$ will be shown
next to $c^{y\rightarrow x}$ and $c^{x\rightarrow y}$. At the end of this
section, these two methods will be compared.

###  III.1 Hénon maps

Hénon map is a discrete-time dynamical system. Based on the values of its two
parameters $a$ and $b$ it can be chaotic or converge to a periodic orbit. We
choose the most commonly studied map with parameter values $a=1.4$ and
$b=0.3$, which yield chaotic dynamics. Two unidirectionally coupled maps are
defined by four equations. The driven subsystem is described by

$\displaystyle x_{1}^{\prime}=a-\big{(}\epsilon
x_{1}y_{1}+(1-\epsilon)x_{1}^{2}\big{)}+bx_{2}$ (19) $\displaystyle
x_{2}^{\prime}=x_{1},$ (20)

and the driving one by

$\displaystyle y_{1}^{\prime}=a-y_{1}^{2}+by_{2}$ (21) $\displaystyle
y_{2}^{\prime}=y_{1}$ (22)

The coupling direction is $y\rightarrow x$. This system was analysed in [25].
We choose time series length $N=2\cdot 10^{4}$ and segment length $L=10$. The
cross-distance vectors at a few values of the coupling parameter $\epsilon$
are shown in Figure4. They behave as expected. In the absence of coupling,
they are both constant. For weak coupling, the initial values of
$\bm{v}^{y\rightarrow x}$ decrease with increased coupling, while this does
not happen with $\bm{v}^{x\rightarrow y}$. For strong coupling, both
$\bm{v}^{y\rightarrow x}$ and $\bm{v}^{x\rightarrow y}$ lose the straight
shape, and the coupling direction cannot be inferred from them.

Figure 4: The cross-distance vectors for unidirectionally coupled Hénon maps
at different values of the coupling parameter $\epsilon$. They are obtained
from the time series of the subsystems’ coordinates $x_{1}$ and $y_{1}$ and
with segment length $L=10$. Figure 5: The indices $c^{y\rightarrow x}$,
$c^{x\rightarrow y}$, $M(Y|X)$, and $M(X|Y)$ (a), and index ratios (b) for
unidirectionally coupled Hénon maps at different values of the coupling
parameter $\epsilon$. The coupling indices $c$ were calculated by (12) with
$k_{1}=10,k_{2}=100$, and the indices $M$ were calculated by (14) with $k=10$.

The dependence of the coupling indices $c^{x\rightarrow y}$ and
$c^{y\rightarrow x}$ on the coupling parameter $\epsilon$ is shown in Figure5
(a). In the absence of coupling, both coupling indices are close to zero. With
increased coupling, $c^{y\rightarrow x}$ increases for small values of
$\epsilon$. The other coupling index $c^{x\rightarrow y}$ is close to zero
until the synchronisation threshold at $\epsilon\approx 0.75$. Therefore, the
cross-distance vectors can be used to infer the coupling direction in
unsynchronised Hénon maps correctly.

###  III.2 Rössler systems

Rössler system is a system of three nonlinear ordinary differential equations,
which define a chaotic continuous-time dynamical system. Two coupled Rössler
subsystems are defined by six equations. We chose the same parameters for the
subsystems as in [13]. The driven subsystem is thus described by

$\displaystyle\dot{x}_{1}=-0.985x_{2}-x_{3}+\epsilon(y_{1}-x_{1})$ (23)
$\displaystyle\dot{x}_{2}=0.985x_{1}+0.15x_{2}$ (24)
$\displaystyle\dot{x}_{3}=0.2+x_{3}(x_{1}-10).$ (25)

and the driving subsystem by

$\displaystyle\dot{y}_{1}=-1.015y_{2}-y_{3}$ (26)
$\displaystyle\dot{y}_{2}=1.015y_{1}+0.15y_{2}$ (27)
$\displaystyle\dot{y}_{3}=0.2+y_{3}(y_{1}-10).$ (28)

Coupling is unidirectional in the direction $y\rightarrow x$. It should be
noted that the first subsystem is in a regular regime for $\epsilon=0$, and in
a chaotic regime for $\epsilon>0$, while the second subsystem is always
chaotic. Therefore, we omit analysis at zero coupling.

Time series parameters are $N=2\cdot 10^{4},\Delta t=0.5$. The system
integration was done with the Runge-Kutta 4 integrator with time step 0.01. We
choose the time series of the subsystems’ coordinates $x_{1}$ and $y_{1}$ and
segment length $L=20$ for the calculation of the cross-distance vectors.

The behaviour of the cross-distance vectors for a few values of the coupling
parameter $\epsilon$ is shown in Figure6. It is similar to the examples seen
so far. For weak coupling, the initial values of $\bm{v}^{y\rightarrow x}$
decrease from the bulk, while $\bm{v}^{x\rightarrow y}$ stays roughly
constant. For strong coupling, both subsystems become identical, and so do
both cross-distance vectors.

The dependence of the coupling indices $c^{x\rightarrow y}$ and
$c^{y\rightarrow x}$ on the coupling parameter $\epsilon$ is shown in Figure7
(a). Similar to the previous example, $c^{y\rightarrow x}$ increases with
increased coupling strength while $c^{x\rightarrow y}$ stays close to zero. At
$\epsilon\approx 0.15$, synchronisation occurs, and the coupling direction can
no longer be inferred from the coupling indices.

This analysis shows that the cross-distance vectors can again be used to infer
the coupling direction in unsynchronised Rössler subsystems correctly.

Figure 6: The cross-distance vectors $\bm{v}^{y\rightarrow x}$ and
$\bm{v}^{x\rightarrow y}$ for unidirectionally coupled Rössler subsystems at
different values of the coupling parameter $\epsilon$. They are obtained from
the time series of the subsystems’ coordinates $x_{1}$ and $y_{1}$ and with
segment length $L=20$. Figure 7: The indices $c^{y\rightarrow x}$,
$c^{x\rightarrow y}$, $M(Y|X)$, and $M(X|Y)$ (a), and index ratios (b) for
unidirectionally coupled Rössler subsystems at different values of the
coupling parameter $\epsilon$. The coupling indices $c$ were calculated by
(12) with $k_{1}=10,k_{2}=100$, and the indices $M$ were calculated by (14)
with $k=10$.

###  III.3 Duffing systems

Duffing system is a periodically forced nonlinear oscillator with damping. It
is a nonautonomous continuous-time dynamical system. It can exhibit chaotic or
periodic dynamics based on the values of its parameters. We chose the
parameters of the two subsystems the same as in [26], resulting in coupled
chaotic subsystems. They are described by

$\displaystyle\ddot{x}+0.2\dot{x}-x+x^{3}=0.3\cos(t)+\epsilon_{1}(y-x)$ (29)
$\displaystyle\ddot{y}+0.3\dot{y}-y+y^{3}=0.5\cos(1.2t)+\epsilon_{2}(x-y).$
(30)

We consider a unidirectional case with $\epsilon_{2}=0$ and a bidirectional
case with $\epsilon_{2}=0.1$.

Time series parameters are $N=5\cdot 10^{4}$ and $\Delta t=0.5$. The system
integration was done with the Runge-Kutta 4 integrator with a time step of
$0.01$. The time series of the first coordinates (positions) and segment
length $L=20$ were chosen for calculating the cross-distance vectors. Since
the behaviour of the cross-distance vectors is very similar to the above
examples, here we only consider the coupling index dependences.

The dependence of the coupling indices $c^{x\rightarrow y}$ and
$c^{y\rightarrow x}$ on the coupling parameter $\epsilon_{1}$ is shown in
Figure8 (a) (unidirectional case) and Figure9 (a) (bidirectional case). In the
unidirectional case, $c^{y\rightarrow x}$ again increases with increased
coupling strength while $c^{x\rightarrow y}$ stays close to zero.
Synchronisation occurs at around $\epsilon_{1}\approx 0.75$. Therefore, cross-
distance vectors can also be used to infer the coupling direction in
nonautonomous Duffing subsystems correctly.

In the bidirectional case, $c^{x\rightarrow y}$ always has a positive value,
while $c^{y\rightarrow x}$ increases with increased $\epsilon_{1}$ similar to
the unidirectional case. At around $\epsilon_{1}=\epsilon_{2}=0.1$, the values
of the coupling indices are the same. Therefore, cross-distance vectors can
also be used to quantify coupling strength in each direction in
bidirectionally coupled subsystems.

Figure 8: The indices $c^{y\rightarrow x}$, $c^{x\rightarrow y}$, $M(Y|X)$,
and $M(X|Y)$ (a), and index ratios (b) for unidirectionally coupled Duffing
subsystems at different values of the coupling parameter $\epsilon_{1}$. The
coupling indices $c$ were calculated by (12) with $k_{1}=10,k_{2}=100$, and
the indices $M$ were calculated by (14) with $k=10$.  Figure 9: The indices
$c^{y\rightarrow x}$, $c^{x\rightarrow y}$, $M(Y|X)$, and $M(X|Y)$ (a), and
index ratios (b) for bidirectionally coupled Duffing subsystems at different
values of the coupling parameter $\epsilon_{1}$ and at $\epsilon_{2}=0.1$. The
sudden increases in all the indices, seen in four values of $\epsilon_{1}$,
are due to system bifurcations. The coupling indices $c$ were calculated by
(12) with $k_{1}=10,k_{2}=100$, and the indices $M$ were calculated by (14)
with $k=10$.

###  III.4 Comparison to established indices

We have discussed the behaviour of the dependences of the coupling indices on
the coupling parameter for various systems. The comparison between the
coupling indices $c$ and the established indices $M$ is shown in Figure5 (a),
Figure7 (a), Figure8 (a) and Figure9 (a). Both index variants were obtained
from the same segments. The choice of the parameter $k$ in (14) was done in a
similar manner as the choice of $k_{1}$ in (12). It turns out that the optimal
value of $k$ is similar to the optimal value of $k_{1}$. This is not
surprising since $k_{1}$ is chosen such that it captures the initial tail,
much like $k$. In all of the examples in this article, $k_{1}=k=10$.

In the unidirectional cases, both indices behave similarly in the sense that
in the direction $y\rightarrow x$, they increase with increased coupling,
while the index in the other direction stays close to zero.

An important difference is between $c^{x\rightarrow y}$ and $M(X|Y)$. Ideally,
these indices should be zero when the coupling direction is $y\rightarrow x$.
The new index $c^{x\rightarrow y}$ is significantly smaller than $M(X|Y)$.
Most importantly, the ratio $c^{y\rightarrow x}/c^{x\rightarrow y}$ is
significantly larger than the ratio $M(Y|X)/M(X|Y)$, which is shown in Figure5
(b), Figure7 (b), and Figure8 (b). As explained in Section II.2, the summation
in (15a) goes up to $N-L+1$, which captures the increasing trend seen in the
cross-distance vectors and leads to larger values of $M(X|Y)$.

Ideally, the index ratios should be infinite. By comparing the coupling
indices, the coupling direction is more accurately determined with the $c$
indices.

In the bidirectional case, the main difference between the two methods is in
the indices of the direction $x\rightarrow y$. Ideally, they should be
constant since $\epsilon_{2}=0.1$ is constant. The coupling parameter
$c^{x\rightarrow y}$ does vary slightly, but less than $M(X|Y)$.
Interestingly, bifurcations have a small impact on the coupling indices $c$
and a large impact on $M$.

##  IV Numerical stability

We have shown that the cross-distance vectors can detect the coupling
direction from measured time series. In this section, we will discuss the
numerical stability and parameter selection of this method. All the analysis
will be done on the test system of unidirectionally coupled Duffing subsystems
with the coupling parameters $\epsilon_{1}=0.1$ and $\epsilon_{2}=0$, where
the system integration is done with the Runge-Kutta 4 integrator and then
sampled at $\Delta t=0.5$. As in the previous section, the position
coordinates’ time series will be used to obtain the cross-distance vectors. We
chose this system because it consists of continuous chaotic subsystems, a
property commonly found in real systems. However, analysis of all systems in
Section III gives similar results. The analysis in this section focuses solely
on the cross-distance vectors and the new index $c$. A detailed comparison of
numerical properties of the indices $c$ and $M$ is done in Appendix B.

###  IV.1 Time series length

An important property is the convergence of the cross-distance vectors with
the length of the time series $N$. In Figure10, the dependence of the cross-
distance vectors on the length of the time series $N$ is shown at constant
$L=20$.

The beginning of the cross-distance vector $\bm{v}^{y\rightarrow x}$ is
significantly smaller than the bulk, i.e., the vector has an initial tail. By
increasing $N$, the first point (represented by the black line) lowers even
further from the bulk, which is desired. This happens because any entry of the
cross-distance vector $\bm{v}^{y\rightarrow x}$ cannot be smaller than the
smallest entry of the distance matrix (9), since the vector entry itself is an
average of these distances. When dealing with a finite number of segments
(finite $N$), the smallest distance between a pair of segments will be a
finite value. With increased $N$, the smallest distance (most likely)
decreases due to a bigger number of segments. Thus, the lowest possible value
of the cross-distance vector also decreases, and so do the initial values of
$\bm{v}^{y\rightarrow x}$. This means that coupling is easier to detect for
longer time series. In the limit, $N\rightarrow\infty$, the initial values of
$\bm{v}^{y\rightarrow x}$ would reach zero.

In the opposite direction with no coupling, the effect is opposite. The
initial values of $\bm{v}^{x\rightarrow y}$ converge towards the bulk of the
vector, i.e., there is no initial tail, which is also desired. This means that
at large enough $N$, no coupling is detected in the direction with no
coupling.

Figure 10:  The dependence of the cross-distance vectors $\bm{v}^{y\rightarrow
x}$ and $\bm{v}^{x\rightarrow y}$ on the length of the time series $N$. The
blue lines represent the cross-distance vectors at a certain $N$ (these are
the same lines as in, for example, Figure2). The black lines represent the
values of the first points of the cross-distance vectors. For a better
visibility, six cross-distance vectors are highlighted and have $N$ that is
written in the figure ticks. The scale on the $N$ axis is logarithmic. The $i$
axis of each plotted cross-distance vector was normalised to $i/N$ (a value
between 0 and 1) for a simpler comparison.

Another behaviour we notice from Figure10 is the change of the whole shape of
the cross-distance vectors. Their smoothness increases with increased $N$, and
they seem to converge to a particular shape (which depends on the system).

To sum up, by increasing the length of the time series $N$, the reliability of
this method increases. For large enough $N$, the coupling will not be detected
in the direction without coupling, and it will be detected in the direction of
coupling. Interestingly, the coupling direction is reliably inferred in short
time series that contain only about 40 oscillations.

###  IV.2 Algorithm parameter dependence

We use the same test system for this analysis. The time series length is
$N=2\cdot 10^{4}$. The only parameter of the cross-distance vectors algorithm
is the segment length $L$. Figure11 (a) shows the dependence of the first
points and of the mean of the cross-distance vectors on the segment length
$L$.

The value of the first point $v_{1}^{y\rightarrow x}$ generally decreases with
increasing $L$ up to $L\approx 200$. Since the mean values of both cross-
distance vectors increase only by a little with increasing $L$, this indicates
that the significant change is in the initial tail. In the other direction,
$v_{1}^{x\rightarrow y}$ stays close to the mean
$v_{\text{mean}}^{x\rightarrow y}$ until it starts to decrease at around
$L\approx 50$. This is reflected in the coupling indices $c^{y\rightarrow x}$
and $c^{x\rightarrow y}$, the dependence of which is shown in Figure11 (b). It
tells us two important properties.

First, we notice that $c^{y\rightarrow x}$ does not detect coupling for $L=1$.
This is because, at $L=1$, the subsystems are not well reconstructed with the
segments. The Takens’ embedding theorem [27] gives a minimum dimension of
delay embedding vectors needed for reconstructing a system’s attractor. It is
$2n+1$, where $n$ is the system dimension (though often, less than that is
needed). For Duffing subsystems in our example, $L=1$ does not reconstruct the
attractor, which results in falsely not detecting coupling.

Second, we can see that for $L>1$, $c^{y\rightarrow x}$ is significantly
larger than $c^{x\rightarrow y}$ and generally increases with increased $L$,
except for a local maximum at $L=10$. The other index $c^{x\rightarrow y}$
stays close to zero until $L\approx 50$, where it starts to increase visibly.
This means that a too large $L$ results in falsely detecting coupling. The
inherent similarity of neighbouring segments is the reason for false coupling
detection in the direction $x\rightarrow y$ at $L>50$.

The similarity between segments will generally decrease if we increase segment
length $L$. This becomes obvious when we notice that with increased $L$, more
dimensions of segments must match to maintain high similarity. However, the
neighbouring segments $\bm{o}^{x}_{i}$ and $\bm{o}^{x}_{i+1}$ are
autocorrelated and will, therefore, always be very similar for any $L$. This
can be seen from the definition of the distance measure (7).

Therefore, when $L$ becomes large enough, the closest segment to
$\bm{o}_{i}^{x}$ will most likely be $\bm{o}_{i+1}^{x}$ (or
$\bm{o}_{i-1}^{x}$). The same will hold for segments of the other subsystem,
i.e., the closest to $\bm{o}_{i}^{y}$ will most likely be $\bm{o}_{i+1}^{y}$.
This makes it seem as if both statements (5) and (6) are technically correct
(but only for time autocorrelated segments), regardless of whether there is
coupling. This is the reason for the increase of $c^{x\rightarrow y}$ at
$L\approx 50$ in our example. It shows that $L$ cannot be too large as it can
result in falsely detecting coupling. A more detailed analysis of this effect
is done in Appendix C. Guidelines for tuning this parameter are given in
Section IV.5.

The detected coupling direction turns out to be robust with respect to the
choice of $k_{1}$ and $k_{2}$ when calculating the coupling indices. The
general rule is that $k_{1}$ should contain the initial tail, which contains
the information about coupling. The choice of $k_{2}$ is significantly less
important. In this article, we chose $10k_{1}<k_{2}<100k_{1}$, though smaller
and larger values give similar results. The point is that $k_{2}$ must be much
smaller than $N$. An analysis of the influence of $k_{1}$ and $k_{2}$ values
on the indices $c$ is done in Appendix D.

Figure 11: The dependence of the first points and the averages of the cross-
distance vectors (a) and the coupling indices (b) on the segment length $L$.
Coupling indices were obtained with equation (12) with $k_{1}=10$ and
$k_{2}=100$.

###  IV.3 Noise dependence

Let us consider the noise robustness of this method. We use time series of the
same test system. The time series length is $N=5\cdot 10^{4}$. We add Gaussian
noise to each point of the time series $x(t_{i})\rightarrow
x(t_{i})+\xi_{i},\xi_{i}\sim\mathcal{N}(0,\sigma^{2})$, where $\sigma$ is the
standard deviation of the noise.

In Figure12, the dependence of the coupling indices $c^{y\rightarrow x}$ and
$c^{x\rightarrow y}$ on $\sigma$ is shown. As one might expect, the
reliability of the detected coupling direction decreases with increased
$\sigma$. For $L=10$, the index $c^{y\rightarrow x}$ decreases to zero at
around $\sigma\approx 0.4$, and for $L=30$ at around $\sigma\approx 0.7$,
which is more than a third of the subsystems’ amplitudes. The other index
$c^{x\rightarrow y}$ is close to zero for all values of $\sigma$.

Noise robustness can be further increased by increasing $L$. However, one must
be careful not to increase it to the point of false detection, explained in
Section IV.2. This method is, therefore, quite robust to noise.

Figure 12: The noise dependence of the coupling indices $c^{y\rightarrow x}$
and $c^{x\rightarrow y}$ with segment lengths $L=10$ (a) and $L=30$ (b).
Coupling index parameters are $k_{1}=10$ and $k_{2}=100$.

###  IV.4 Robustness to artefacts and missing data

Measurements can contain artefacts, such as spikes, or have missing values for
a time period. Coupling detection in such data can be problematic. Coupling
indices are resilient to such imperfections in time series.

The solution is to delete the distance matrix elements whose values were
obtained from segments with artefacts or missing values. This is done by
deleting the corresponding rows and columns, decreasing the matrix size. They
must be deleted from both distance matrices, even if the artefact is only
present in one time series. This allows us to ignore any unwanted sections of
either time series, making this approach very flexible. The only requirement
is manually choosing the time series points to ignore.

To give an example, consider the time series of the same test system. The time
series length is $N=5\cdot 10^{4}$. We change the values of $x$ between points
$5000$ and $9000$ and the values of $y$ between $30000$ and $36000$ to random
values to simulate artefacts. Therefore, when computing the distance matrices,
we ignore $4000+6000=10000$ rows and columns (which will result in cross-
distance vectors with $10000$ fewer points). Let us compare cross-distance
vectors, $c$ indices, and $M$ indices obtained from clean time series and time
series containing artefacts.

The results are shown in Figure 13. The cross-distance vectors obtained from
clean data and data containing artefacts are nearly indistinguishable. The $c$
and the $M$ indices are also nearly the same. The difference in $c$ is seen in
the third decimal. The difference in $M$ is larger but still small. This shows
that state space approaches are resilient to artefacts and missing values.

Figure 13: Cross-distance vectors obtained from clean data (C) and from data
with artefacts (A). $\bm{v}^{y\rightarrow x}$ is shown in (a) and
$\bm{v}^{x\rightarrow y}$ is shown in (b). The $j$ axis of each plotted cross-
distance vector was normalised to $j/N$ (a value between 0 and 1) for a
simpler comparison. Note that the vectors obtained from clean and from data
with artefacts are nearly indistinguishable. The index values obtained from
clean data are $c^{y\rightarrow x}=0.268$, $c^{x\rightarrow y}=0.000$,
$M(Y|X)=0.407$, $M(X|Y)=0.052$. The index values obtained from data with
artefacts are $c^{y\rightarrow x}=0.263$, $c^{x\rightarrow y}=-0.001$,
$M(Y|X)=0.385$, $M(X|Y)=0.049$.

###  IV.5 Practical computational details

In this section, we will briefly discuss the practical aspects of the
algorithm and provide a suggestion that can prove useful.

#### Tuning the L parameter

For this purpose, we can use the large similarity of neighbouring segments
(explained in Section IV.2). For a chosen $L$, we can check whether the
closest segment is often one of the neighbours, i.e., whether the distance
matrices (8) and (9) have very small subdiagonals. If at least one of them
does, we must lower $L$ until subdiagonals have similar values to the rest of
the matrices since small subdiagonals lead to false detection. This gives the
maximal possible value $L_{max}$. In practice, we recommend a value close to
the maximal $L_{max}$, for example, $L_{max}/2$, since larger $L$ generally
gives more accurate and noise-robust results.

#### Neighbouring segments

In some cases, especially in time series with a large sampling rate (small
time step $\Delta t$), lowering $L$ cannot adequately raise the values of the
subdiagonals. If the time step is very small, the neighbouring segments will
be autocorrelated and will therefore always be very close such that
$d(\bm{o}_{i}^{x},\bm{o}_{i+1}^{x})\rightarrow 0$. This can lead to false
coupling detection, as explained in IV.2. We propose two solutions.

The first solution is to downsample the signal. The neighbouring segments
become less similar by increasing the time step $\Delta t$. We must downsample
the time series to the point where neighbouring segments are no longer the
most similar, i.e., when the subdiagonal values of the distance matrices are
no longer small compared to the rest of the values.

The second solution is manually decreasing the dimensions of the distance
matrices. This is needed in a case where downsampling is not an option. One
such example is stiff subsystems, i.e., subsystems which contain small and
large frequency components. In such cases, downsampling can erase the high-
frequency component (and thus erase possibly crucial information) but still
not raise the values of the subdiagonals. For this purpose, we suggest an
alternative approach. We can only calculate every $M$-th value of the full
distance matrices, i.e., we construct them from every $M$-th segment. Choosing
$M$ large enough ignores the (small) subdiagonals that appear in full
matrices.

Both approaches should be done with a fixed $L$, for example, $L=5$. Once the
subdiagonals are adequately raised, the chosen approach can be repeated for a
larger $L$. If successful, the larger $L$ should be taken for better accuracy
and noise robustness.

Both approaches ignore small subdiagonals that would appear in full matrices,
as well as other values that can be redundant. For example, a full matrix
contains elements $d(\bm{o}_{i}^{x},\bm{o}_{j}^{x})$ and
$d(\bm{o}_{i}^{x},\bm{o}_{j+1}^{x})$, which have similar values if
$d(\bm{o}_{j}^{x},\bm{o}_{j+1}^{x})\approx 0$.

The downsampling approach is computationally more efficient, as it effectively
decreases segment length $L$. Therefore, if possible, it should be chosen over
manually decreasing dimensions.

#### Time and memory limitations

In a time series of length $N$, the computational complexity of the algorithm
is $\mathcal{O}(N^{2}\log(N))$ due to the sorting of all rows of the distance
matrices. As the previous paragraph explains, downsampling the original time
series is recommended since it significantly decreases the execution time.

Furthermore, the matrices $\bm{D}^{x}$ and $\bm{D}^{y\rightarrow x}$ can, in
practice, be too large to store in computer memory. To avoid this, the rows of
the matrices (8) and (9) can be calculated individually to obtain a single row
of $\bm{D}^{y\rightarrow x}$. One summation in (11) is done for each row, and
the cross-distance vector is obtained by only storing a few sets of data of
size $N-L+1$ at once.

#### GPU implementation

Since the algorithm is based on matrix operations, GPU devices can
significantly decrease the execution time. For this purpose, we are providing
a GPU-based implementation written in Python using the JAX library [28]. A
time-efficient and memory-efficient implementation is available at [29].

##  V Application to physiological signals

In this section, we will apply the cross-distance vectors algorithm to time
series of physiological measurements. The goal is to show the applicability of
cross-distance vectors to real-life systems. We will analyse coupling in the
human cardiorespiratory system. The cardiac subsystem is characterised by ECG
(electrocardiogram), and the respiratory subsystem by respiration curves
obtained with a respiratory belt sensor. These subsystems are inherently
bidirectionally coupled, i.e., both subsystems depend on each other. A survey
of the mechanisms responsible for this dependence is given in [30].
Specifically, we will investigate coupling in the direction from the
respiratory to the cardiac subsystem, which we label with $R\rightarrow C$.

The signals used in the analysis are shown in Figure14. Thirty minutes of ECG
and respiration curves were measured on two subjects. In such a way, we obtain
two pairs of signals, marked $R_{1}$ and $C_{1}$ for the first subject and
$R_{2}$ and $C_{2}$ for the second subject.

Figure 14: A part of the ECG and respiration signals used for analysis. Label
$C$ means ECG, and label $R$ means respiration. $R_{1}$ and $C_{1}$ were
measured simultaneously on one subject, and $R_{2}$ and $C_{2}$ simultaneously
on another subject. Time series parameters are $N=172800$ and $\Delta
t=\frac{1}{96}s$, corresponding to a measurement time of $30$ minutes.

The most dominant frequency component in ECG signals is around 1Hz. The R peak
is, however, very short, which means this subsystem is stiff. Therefore we do
not downsample the time series but rather manually decrease the distance
matrix sizes by a factor $M$, as explained in Section IV.5. In this case, we
chose $M=20$, which means we are comparing segments that are apart by a
multiple of $M\Delta t=20/96s\approx 0.2s$.

Figure15 (a) and (b) represent cross-distance vectors for $R_{1}\rightarrow
C_{1}$ and $R_{2}\rightarrow C_{2}$. In these cases, ECG and respiratory
signals belong to the same person and were measured simultaneously. Therefore,
we expect to detect coupling. Indeed, the cross-distance vectors have an
initial tail, suggesting coupling is present in the underlying subsystems.

Figure15 (c) and (d) represent cross-distance vectors for $R_{2}\rightarrow
C_{1}$ and $R_{1}\rightarrow C_{2}$. Since the signals belong to different
subjects that are inherently independent, we do not expect to detect coupling.
Indeed, the cross-distance vectors do not have an initial tail, suggesting no
coupling in the underlying subsystems.

Let us compare the $c$ indices, obtained with parameters
$k_{1}=10,k_{2}=1000$. The indices $c^{R_{1}\rightarrow C_{1}}=0.105$ and
$c^{R_{2}\rightarrow C_{2}}=0.176$ are nearly two orders of magnitude larger
than $c^{R_{2}\rightarrow C_{1}}=-0.0004$ and $c^{R_{1}\rightarrow
C_{2}}=-0.006$. This means that coupling indices have significantly larger
values when the coupling is present in the underlying subsystems than when the
subsystems are independent. Coupling indices $c^{R_{2}\rightarrow C_{1}}$ and
$c^{R_{1}\rightarrow C_{2}}$ even have a small negative value, which is a
strong indication of the absence of coupling.

For comparison, the $M$ indices were also computed with parameter $k=10$. They
were computed from the same distance matrices as the $c$ indices, i.e.,
matrices obtained by comparing every 20th segment with $L=20$. The obtained
values are $M(R_{1}|C_{1})=0.150$, $M(R_{2}|C_{2})=0.254$,
$M(R_{2}|C_{1})=-0.003$, and $M(R_{1}|C_{2})=-0.009$. The $M$ indices also
correctly detect the existence of coupling in the first two cases and have a
small negative value in the second two cases.

When considering the coupling indices $c$ in an application, it should be
noted that they require the time series to belong to a dynamical system that
can be reconstructed with time-delay embedding. Otherwise, approaches such as
Granger causality or information transfer are more suitable. However, contrary
to the indices $c$, these approaches may not be reliable in short time series.
To sum up, the new coupling indices $c$ excel in causality detection from
bivariate time series generated by a dynamical systems, especially in short
signals.

Figure 15: The first 1000 points of the cross-distance vectors in the
direction from the respiratory to the cardiac subsystem. In (a) and (b), the
analysed signals belong to the same subject, while in (c) and (d), they are
independent. The segment length is $L=20$. Dimensions of distance matrices
were decreased by a factor $20$.

##  VI Conclusion

In this article, we have defined cross-distance vectors as a means of
inferring the direction of coupling from bivariate time series and provided an
algorithm for calculating them. Cross-distance vectors provide information
about coupling by calculating two coupling indices, which quantify coupling
strength in both directions.

The new coupling indices can infer the coupling direction in various coupled
dynamical systems. Comparing the new coupling indices to the conventional
indices based on state space distance shows more accurate results with the new
approach. Analysis of numerical stability has shown that the reliability of
the new indices increases with the length of the time series but is also
reliable in short time series, containing only about 40 oscillations. The
performance of the coupling indices was also evaluated in detecting
cardiorespiratory interaction in measured data.

The algorithm has a simple implementation that requires the choice of only a
maximum of three parameters. The selection of the optimal parameters’ values
can be made systematically, which results in robust performance and
circumvents the challenge of determining optimal embedding parameters. A
numerically efficient implementation is available at [29].

A logical continuation is an extension of the new method for inferring
connections in a network of multiple subsystems. Should the multivariate
extension be as effective as the bivariate variant, it could become an
essential tool in analysing complex multivariate problems such as brain
connectivity.

##  Appendix A Mapping of close states

In what follows, we will provide an argument to illustrate why close states in
the driven subsystem are mapped to close states in the driving subsystem for
various subsystems with different properties.

One might assume that the statement (5) is false and the statement (6) is
true. It may seem that the times when the driving subsystem is self-similar
will appear as times of large self-similarity in the driven subsystem via
coupling. Hence (6) is true. Since this logic does not apply in the other
direction, one might expect that the statement (5) is false. This is, however,
not the case. In a special case under two assumptions, we provide an analytic
argument. For a more general case, we will provide a heuristic argument.

#### Special case

Consider that two segments of the driven subsystem $\bm{o}_{i}^{x}$ and
$\bm{o}_{j}^{x}$ are identical (this is only possible for subsystems in a
regular dynamical regime)

$\displaystyle\bm{x}(t_{i+k})=\bm{x}(t_{j+k}),$ (31) $\displaystyle
k=0,1,\dots,L-1.$

If the subsystem coordinates match over a time period, their time evolution
must also match. The first assumption is that $\bm{f}$ is time-independent,
which gives us the following $L$ equations

$\begin{split}&\bm{f}(\bm{x}(t_{i+k}))+\bm{g}(\bm{x}(t_{i+k}),\bm{y}(t_{i+k}))=\\\
&\bm{f}(\bm{x}(t_{j+k}))+\bm{g}(\bm{x}(t_{j+k}),\bm{y}(t_{j+k})),\\\
&k=0,1,\dots,L-1.\end{split}$ (32)

Additionally, from (31) follows

$\displaystyle\bm{f}(\bm{x}(t_{i+k}))=\bm{f}(\bm{x}(t_{j+k}))$ (33)
$\displaystyle k=0,1,\dots,L-1.$

Combining (32) and (33) gives us

$\displaystyle\bm{g}(\bm{x}(t_{i+k}),\bm{y}(t_{i+k}))=\bm{g}(\bm{x}(t_{j+k}),\bm{y}(t_{j+k})),$
(34) $\displaystyle k=0,1,\dots,L-1.$

Since (31) holds, we can consider the first argument of $\bm{g}$ as constant
and define $k$ functions

$\displaystyle\tilde{\bm{g}}(\bm{y}(t_{i+k}))=\bm{g}(\bm{x}(t_{i+k}),\bm{y}(t_{i+k})),$
(35) $\displaystyle k=0,1,\dots,L-1.$

The second assumption is that functions $\tilde{\bm{g}}(\bm{y}(t_{i+k}))$ are
injective. If this is the case, from (34) follows

$\displaystyle\bm{y}(t_{i+k})=\bm{y}(t_{j+k}),$ (36) $\displaystyle
k=0,1,\dots,L-1.$

This means that statement (5) is true. This argument does not work in the
other direction with swapped $\bm{x}$ and $\bm{y}$, since $\bm{y}$ does not
depend on $\bm{x}$, i.e., there is no coupling function in the time evolution
of $\bm{y}$. Therefore, the statement (6) is false.

#### General case

We argue that if the subsystems are sufficiently nice, this also holds without
the two assumptions if segment length $L$ is large enough. To obtain (34) from
(32) with time dependent $\bm{f}(\bm{x},t)$, we argue that it would seem
unlikely for the sum of $\bm{f}$ and $\bm{g}$ to match over a time period, if
they do not match individually. The analysis in Section III.3 shows that this
indeed holds for coupled periodically forced Duffing subsystems. To obtain
(36) from (34) if $\tilde{\bm{g}}(\bm{y}(t_{i+k}))$ are not injective, we
would similarly argue that if $\bm{g}$ matches over a long enough time period,
so must its arguments.

We can use the same arguments for chaotic subsystems by demanding the distance
of segments to be less than $\delta$ instead of them being identical. Under
the admissibility conditions of the Poincaré recurrence theorem, the choice of
$\delta$ can be arbitrarily small. The arguments hold for chaotic subsystems
by swapping all the equalities in equations (31)–(36) with arbitrarily small
proximity.

##  Appendix B Detailed comparison of $c$ and $M$

Here the goal is to compare the numerical stability of $c$ and $M$ indices.
For that purpose, the analysis from Section IV is repeated for $M$, and both
results are presented jointly. All the analysis is done on the same test
system of coupled Duffing oscillators as in Section IV. Also, all the
parameters are the same as in Section IV. Coupling is unidirectional with
direction $y\rightarrow x$.

###  B.1 Dependence on time series length

The analysis in Section IV.1 was done on the cross-distance vectors. For
results to be comparable to $M$, we compute the $c$ indices. As explained in
Section III.4, the indices are most fairly compared when their parameters are
$k=k_{1}$. In our case, $k=k_{1}=10$ and $k_{2}=100$ are taken, the same as in
Section III.3. While these are not optimal parameters for every $N$, we keep
them fixed for simplicity. The segment length is $L=20$, as in Section IV.1.
The comparison is done in Figure 16.

Figure 16: The dependence of the indices $c^{y\rightarrow x}$,
$c^{x\rightarrow y}$, $M(Y|X)$, $M(X|Y)$ (a) and the index ratios
$c^{y\rightarrow x}/c^{x\rightarrow y}$, $M(Y|X)/M(X|Y)$ (b) on the length of
the time series $N$. The $M$ index parameter is $k=10$ and the $c$ index
parameters are $k_{1}=10,k_{2}=100$.

The general behaviour is similar for both indices, as seen in Figure 16(a). At
small $N$, all indices are large. The main difference between $M$ and $c$ is
seen at large $N$ where $c^{x\rightarrow y}$ converges to zero, which is
desirable, while $M(X|Y)$ seems to converge to a finite positive value. As
explained in Section II.2, this is because $c$ ignores possible trends in the
cross-distance vectors, while $M$ does not.

The index ratios are shown in Figure 16(b). Ideally, the ratios are infinite
since the coupling is unidirectional. For nearly all values of $N$ the ratio
$c^{y\rightarrow x}/c^{x\rightarrow y}$ is larger than the ratio
$M(Y|X)/M(X|Y)$, especially for large $N$. For small $N$, the difference is
smaller, but the $c$ ratio is generally still larger. This shows that $c$
better determines the coupling direction for both short and long time series.

###  B.2 Dependence on $L$

Let us compare the $L$ dependence from Figure 11 for indices $c$ and $M$. The
chosen time series length is $N=2\cdot 10^{4}$, same as in Section IV.2.

The dependence of the indices $c$ and $M$ on the segment length $L$ is shown
in Figure 17(a). The dependence is very similar for both index variants. The
main difference is that both $M$ indices are slightly larger (since they do
not ignore trends seen in cross-distance vectors). At $L=1$, the $M$ indices
both have a very similar positive value, while the $c$ indices are both much
closer to zero. Interestingly, the bias of $M(X|Y)$ due to trends is nearly
constant for segment length up to $L\approx 50$, where $L$ becomes too large
(as explained in Section IV.2).

The index ratios $c^{y\rightarrow x}/c^{x\rightarrow y}$ and $M(Y|X)/M(X|Y)$
are shown in Figure 17(b). The ratio of the $c$ indices is significantly
larger than that of the $M$ indices for all values of $L$, except for $L=1$.
This means that the $c$ indices better determine the coupling direction
regardless of the choice of segment length $L$.

Figure 17: The dependence of the $c$ and $M$ indices (a) and the index ratios
$c^{y\rightarrow x}/c^{x\rightarrow y}$ and $M(Y|X)/M(X|Y)$ (b) on the segment
length $L$. The $M$ index parameter is $k=10$ and the $c$ index parameters are
$k_{1}=10,k_{2}=100$.

###  B.3 Dependence on noise

Let us compare the noise robustness of the $c$ and the $M$ indices similarly
as in Section IV.3. The chosen time series length is again $N=5\cdot 10^{4}$
and the chosen segment length is $L=10$, the same as in Figure 12(a).

The dependence of the $c$ and $M$ indices on the standard deviation of noise
$\sigma$ is shown in Figure 18(a). It turns out that the dependence is very
similar for both indices. They are both robust to noise. As seen in Figure
18(b), the index ratio is larger for $c$ than for $M$ up to around
$\sigma\approx 0.5$, at which point both ratios become close to $1$.

Figure 18: The dependence of the $c$ and $M$ indices (a) and the index ratios
$c^{y\rightarrow x}/c^{x\rightarrow y}$ and $M(Y|X)/M(X|Y)$ (b) on the
standard deviation of noise $\sigma$. The $M$ index parameter is $k=10$ and
the $c$ index parameters are $k_{1}=10,k_{2}=100$. The segment length is
$L=10$.

##  Appendix C The $N$ and $L$ dependence of cross-distance vectors

In order to get a complete picture of the behaviour of the cross-distance
vectors, we compute them for numerous values of pairs $(N,L)$. We again use
the test system from Section IV. We are interested in the cross-distance
vectors’ first (nonzero) point. The results are shown in Figure 19. Figure 19
can be understood as plotting the black lines in Figure 10, calculated for
different $L$. Equivalently, Figure 19(a) can be seen as plotting
$v_{1}^{y\rightarrow x}$ from Figure 11(a) and Figure 19(b) as
$v_{1}^{x\rightarrow y}$ from Figure 11(a), both calculated for different $N$.

Figure 19(a) shows $v_{1}^{y\rightarrow x}$, which corresponds to the
direction of coupling. One can roughly identify three areas in the $(N,L)$
grid that are highlighted using two red lines. Mind that the areas are not
strictly defined but assist in explaining the figure. Area 1 has very small
values. There, $L$ is too large at a given $N$, which results in detecting
coupling regardless of the underlying dynamics (area of falsely detecting
coupling). Area 3 has large values. There $L=1$, which is too small (at any
$N$) to reconstruct the underlying state space, which results in not detecting
coupling regardless of the underlying dynamics (area of falsely not detecting
coupling). Area 2 has moderate values. There, $L$ is large enough to
reconstruct the underlying state space but small enough to avoid false
detection.

The white line highlights the approximate area where the optimal $L$ at a
given $N$ is. We understand the optimal $L$ as the value in area 2 where
$v_{1}^{y\rightarrow x}$ is the smallest at a given $N$. The optimal $L$
increases with increased $N$.

Figure 19(b) shows $v_{1}^{x\rightarrow y}$, which corresponds to the
direction without coupling. In this case, the three areas and the optimal $L$
are not as obvious as in the other direction. We only clearly see where $L$
becomes too large, which leads to false coupling detection.

It should be noted that in both figures, there is a white area in the upper
left corner. This areas represents the impossible pairs $(N,L)$, i.e., $L\geq
N$.

Figure 19: The first point of the cross-distance vectors $v_{1}^{y\rightarrow
x}$ (a) and $v_{1}^{x\rightarrow y}$ (b) for different values of the time
series length $N$ and the segment length $L$. In (a), the two red lines
separate the grid into three areas, and the white line approximately
represents optimal $L$ at a given $N$. The white area in the upper left corner
represents the impossible pairs $(N,L)$ where $L\geq N$.

##  Appendix D The choice of parameters $k_{1}$ and $k_{2}$

In order to analyse the influence of the parameters $k_{1}$ and $k_{2}$ on the
coupling indices $c$, we consider the same test system used in Section IV.
Specifically, we will plot the $c$ indices from Figure 8(a) (blue and orange
lines), obtained with different $k_{1}$ and $k_{2}$.

The influence of the change of $k_{1}$ is shown in Figure 20. In this case, we
fixed the parameter $k_{2}=150$, which is close to $k_{2}$ most often used in
the article. We can make some important observations. First, we notice that
for small values of $k_{1}$, the values of $c^{y\rightarrow x}$ increase
quickly with increased $\epsilon$, which allows for more reliable detection at
small coupling. Furthermore, the values of $c^{y\rightarrow x}$ reach higher
values, while the values of $c^{x\rightarrow y}$ remain unchanged, which is
also desirable. However, there is a negative side to using such small values
of $k_{1}$. The variance of both indices is not negligible, as seen by a
larger spread of $c^{x\rightarrow y}$ around zero and a negative value of
$c^{y\rightarrow x}$ at $\epsilon=0$. On the other hand, increasing $k_{1}$
decreases the variance of both indices. However, it also reduces the rate of
increase of $c^{y\rightarrow x}$ and introduces a bias at moderate values of
$\epsilon$ (around $0.5$). Therefore, the choice of $k_{1}$ offers a tradeoff
between bias and variance of the indices $c$.

The influence of the change of $k_{2}$ is shown in Figure 21. Similar to the
influence of $k_{1}$, the variance of both indices is larger at smaller
$k_{2}$. Furthermore, larger $k_{2}$ increases $c^{x\rightarrow y}$ at
moderate values of $\epsilon$, increasing the false positive error. Most
importantly, this false positive error increases drastically when $k_{2}$ is
close to the time series length $N$ ($k_{2}=49000$ in our case). One can
notice a strong similarity between the case of $k_{2}=49000$ and the $M$
indices in Figure 8(a) (green and red lines). Based on this analysis, we draw
the conclusion that the safe choice for $k_{2}$ is approximately
$10k_{1}<k_{2}<100k_{1}$.

Figure 20: The dependence of the coupling indices $c$ on the coupling
parameter $\epsilon$ for different values of the parameter $k_{1}$ and fixed
$k_{2}=150$. Figure 21: The dependence of the coupling indices $c$ on the
coupling parameter $\epsilon$ for different values of the parameter $k_{2}$
and fixed $k_{1}=10$.

###### Acknowledgements.

The authors acknowledge the project J3-4525 and the research core funding No.
P2-0001 that were financially supported by the Slovenian Research Agency.

## Author declarations

The experiment from Section V was approved by the Ethical Committee of the
University Medical Centre Ljubljana, Slovenia. All subjects provided written
informed consent.

## Data availability

The data that support the analysis of this article have been generated by the
authors and can be fully reproduced from the repository
https://repo.ijs.si/mbresar/cd-vec [29].

## References

  * Runge  _et al._ [2015] J. Runge, M. Riedl, A. Müller, H. Stepan, J. Kurths, and N. Wessel, Quantifying the causal strength of multivariate cardiovascular couplings with momentary information transfer, Physiological Measurement 36, 813 (2015).
  * Schulz  _et al._ [2013] S. Schulz, F.-C. Adochiei, I.-R. Edu, R. Schroeder, H. Costin, K.-J. Bär, and A. Voss, Cardiovascular and cardiorespiratory coupling analyses: a review, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 371, 20120191 (2013).
  * Wibral  _et al._ [2014] M. Wibral, R. Vicente, and J. T. Lizier,  _Directed Information Measures in Neuroscience_ (Springer Berlin Heidelberg, 2014).
  * Runge  _et al._ [2019] J. Runge, S. Bathiany, E. Bollt, G. Camps-Valls, D. Coumou, E. Deyle, C. Glymour, M. Kretschmer, M. D. Mahecha, J. Muñoz-Marí, E. H. van Nes, J. Peters, R. Quax, M. Reichstein, M. Scheffer, B. Schölkopf, P. Spirtes, G. Sugihara, J. Sun, K. Zhang, and J. Zscheischler, Inferring causation from time series in earth system sciences, Nature Communications 10, 10.1038/s41467-019-10105-3 (2019).
  * Sugihara  _et al._ [2012] G. Sugihara, R. May, H. Ye, C. hao Hsieh, E. Deyle, M. Fogarty, and S. Munch, Detecting causality in complex ecosystems, Science 338, 496 (2012).
  * Ye  _et al._ [2015] H. Ye, E. R. Deyle, L. J. Gilarranz, and G. Sugihara, Distinguishing time-delayed causal interactions using convergent cross mapping, Scientific Reports 5, 10.1038/srep14750 (2015).
  * Granger [1969] C. W. J. Granger, Investigating causal relations by econometric models and cross-spectral methods, Econometrica 37, 424 (1969).
  * Hoover [2008] K. D. Hoover, Causality in economics and econometrics, in  _The New Palgrave Dictionary of Economics_ (Palgrave Macmillan UK, 2008) pp. 1–13.
  * Papana  _et al._ [2021] A. Papana, E. Siggiridou, and D. Kugiumtzis, Detecting direct causality in multivariate time series: A comparative study, Communications in Nonlinear Science and Numerical Simulation 99, 105797 (2021).
  * Clemson  _et al._ [2016] P. Clemson, G. Lancaster, and A. Stefanovska, Reconstructing time-dependent dynamics, Proceedings of the IEEE 104, 223 (2016).
  * Baccalá and Sameshima [2001] L. A. Baccalá and K. Sameshima, Partial directed coherence: a new concept in neural structure determination, Biological Cybernetics 84, 463 (2001).
  * Schreiber [2000] T. Schreiber, Measuring information transfer, Physical Review Letters 85, 461 (2000).
  * Paluš and Vejmelka [2007] M. Paluš and M. Vejmelka, Directionality of coupling from bivariate time series: How to avoid false causalities and missed connections, Physical Review E 75, 056211 (2007).
  * Rosenblum and Pikovsky [2001] M. G. Rosenblum and A. S. Pikovsky, Detecting direction of coupling in interacting oscillators, Physical Review E 64, 045202 (2001).
  * Paluš and Stefanovska [2003] M. Paluš and A. Stefanovska, Direction of coupling from phases of interacting oscillators: An information-theoretic approach, Physical Review E 67, 055201 (2003).
  * Bandrivskyy  _et al._ [2004] A. Bandrivskyy, A. Bernjak, P. McClintock, and A. Stefanovska, Wavelet phase coherence analysis: Application to skin temperature and blood flow, Cardiovascular Engineering 4, 89 (2004).
  * Jamšek  _et al._ [2007] J. Jamšek, A. Stefanovska, and P. V. E. McClintock, Wavelet bispectral analysis for the study of interactions among oscillators whose basic frequencies are significantly time variable, Physical Review E 76, 046221 (2007).
  * Arnhold  _et al._ [1999] J. Arnhold, P. Grassberger, K. Lehnertz, and C. Elger, A robust method for detecting interdependences: application to intracranially recorded EEG, Physica D: Nonlinear Phenomena 134, 419 (1999).
  * Andrzejak  _et al._ [2003] R. G. Andrzejak, A. Kraskov, H. Stögbauer, F. Mormann, and T. Kreuz, Bivariate surrogate techniques: Necessity, strengths, and caveats, Physical Review E 68, 066202 (2003).
  * Chicharro and Andrzejak [2009] D. Chicharro and R. G. Andrzejak, Reliable detection of directional couplings using rank statistics, Physical Review E 80, 026217 (2009).
  * Guo  _et al._ [2008] S. Guo, A. K. Seth, K. M. Kendrick, C. Zhou, and J. Feng, Partial granger causality—eliminating exogenous inputs and latent variables, Journal of Neuroscience Methods 172, 79 (2008).
  * Sun and Bollt [2014] J. Sun and E. M. Bollt, Causation entropy identifies indirect influences, dominance of neighbors and anticipatory couplings, Physica D: Nonlinear Phenomena 267, 49 (2014).
  * Strogatz [2015] S. H. Strogatz,  _Nonlinear Dynamics and Chaos_ (Taylor & Francis Inc, 2015).
  * Ott [2002] E. Ott,  _Chaos in Dynamical Systems_ (Cambridge University Press, 2002).
  * Paluš  _et al._ [2001] M. Paluš, V. Komárek, Z. Hrnčíř, and K. Štěrbová, Synchronization as adjustment of information rates: Detection from bivariate time series, Physical Review E 63, 046211 (2001).
  * Brešar  _et al._ [2021] M. Brešar, P. Boškoski, and M. Horvat, Detection of coupling in duffing oscillator systems, Chaos: An Interdisciplinary Journal of Nonlinear Science 31, 063130 (2021).
  * Takens [1981] F. Takens, Detecting strange attractors in turbulence, in  _Lecture Notes in Mathematics_ (Springer Berlin Heidelberg, 1981) pp. 366–381.
  * [28] Jax library, https://github.com/google/jax.
  * [29] Python implementation of the cross-distance vectors algorithm, https://repo.ijs.si/mbresar/cd-vec.
  * Dick  _et al._ [2014] T. E. Dick, Y.-H. Hsieh, R. R. Dhingra, D. M. Baekey, R. F. Galán, E. Wehrwein, and K. F. Morris, Cardiorespiratory coupling, in  _Progress in Brain Research_ (Elsevier, 2014) pp. 191–205.