# Supplementary Text for ‘Unifying Pairwise Interactions in Complex Dynamics’

Oliver M. Cliff  School of Physics, The University of Sydney, Camperdown NSW
2006, Australia.  Centre for Complex Systems, The University of Sydney,
Camperdown NSW 2006, Australia.  Annie G. Bryant  School of Physics, The
University of Sydney, Camperdown NSW 2006, Australia.  Centre for Complex
Systems, The University of Sydney, Camperdown NSW 2006, Australia.  Joseph T.
Lizier  Centre for Complex Systems, The University of Sydney, Camperdown NSW
2006, Australia.  Faculty of Engineering, The University of Sydney, Camperdown
NSW 2006, Australia.  Naotsugu Tsuchiya  Turner Institute for Brain and Mental
Health & School of Psychological Sciences, Faculty of Medicine, Nursing, and
Health Sciences, Monash University, Melbourne, Victoria, Australia.  Center
for Information and Neural Networks (CiNet), National Institute of Information
and Communications Technology (NICT), Suita-shi, Osaka 565-0871, Japan.
Advanced Telecommunications Research Computational Neuroscience Laboratories,
2-2-2 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-0288, Japan.  Ben D. Fulcher
School of Physics, The University of Sydney, Camperdown NSW 2006, Australia.
Centre for Complex Systems, The University of Sydney, Camperdown NSW 2006,
Australia.

This supplementary document includes detailed information about the statistics
for pairwise interactions (SPIs) in Sec. S1 and the data included in our
multivariate time series (MTS) library in Sec. S2. We also provide additional
information about the empirical similarity index in Sec. S3.1.

##  S1 Library of statistics for pairwise interactions (SPIs)

This section introduces our library of 237 SPIs that we have collected for
evaluating pairwise interactions in multivariate time series (MTS). Each
method is a real-valued measurement, $s(x,y)\in\mathbb{R}$, of some type of
dependency between two time series, $x$ and $y$, and as such are referred to
as statistics for pairwise interactions (SPIs). We only include methods that
can be computed on continuous-valued $M$-variate time series (see Sec. S2),
which can be used to compute an $M$-by-$M$ square matrix of pairwise
interactions (referred to as an MPI, see Sec. S3.1). Prior to evaluating each
statistic, we standardize ($z$-score) the MTS along the time axis. We include
methods that both operate directly on pairs of time series (i.e., bivariate
time series, like Kendall’s $\tau$ [1]) and methods that operate directly on
the full MTS (like precision matrices); both are used to generate an MPI from
an MTS.

Most SPIs in our library are based on existing implementations, but in some
cases, especially where the output of a method was not a single real number,
we implemented new real-valued statistics. Moreover, many of the algorithms
include a number of free parameters that we set either using optimization
procedures (often available within implementations) or fix to a small number
of sensible predefined settings. The combination of both the different
parameter configurations and the different summary statistics taken from each
method gives each SPI in the library a unique identifier (as a string) that
will be used throughout this supplement. Consider the SPI identifier
‘xcorr_mean_sig-True’ (cf. Sec. S1.1.5); here, xcorr refers to the method of
cross-correlation between $x$ and $y$, which itself does not provide a single
statistic but rather a correlogram (i.e., a series of correlation coefficients
for each time delay between the two signals). However, the two additional
modifiers that are separated by underscores in the identifier, mean and sig-
True, gives a scalar value. The first modifier, mean, indicates that we are
taking the average across lags of the cross-correlation function. The second
modifier, sig-True, indicates that we will only take the mean over
statistically significant lags. Therefore, the combination of the method
(xcorr) with the modifiers (mean and sig-True) gives a scalar value from the
bivariate time series $(x,y)$. By using different parameters and modifiers of
45 distinct methodologies, we obtain the 237 SPIs used throughout this work.

In order to simplify our presentation of the methods (and their modifiers)
that are used to generate these SPIs, we have organized each method into one
of the following broad literature categories: ‘basic’, ‘distance similarity’,
‘causal’, ‘information theory’, ‘spectral’, and ‘miscellaneous’. Each of these
categories is included as a subsection below. Note that this grouping is for
convenience and is not unambiguous—as we find throughout this work, many
existing statistics fit multiple descriptions and could thus be organized
validly into multiple categories. Finally, for each SPI, we include a small
number of keywords to indicate whether the interactions they measure are:
undirected or directed; linear or nonlinear; signed or unsigned; bivariate or
multivariate; and contemporaneous, time-dependent, frequency-dependent, or
time-frequency dependent. These keywords help users to understand key
assumptions of each of the SPIs. Note that our aim for the pyspi library is
for it to be an evolving library that is refined and improved over time, with
this section capturing an early snapshot of it.

###  S1.1 Basic statistics (21 SPIs)

In this section, we detail the 21 SPIs that we have categorized as ‘basic
statistics’.

####  S1.1.1 Covariance (cov)

Keywords: undirected, linear, signed, multivariate, contemporaneous.

The covariance matrix is estimated for a wide variety of statistical
procedures. Due to the $z$-scoring of the time series, the correlation and
covariance matrices are equivalent and thus the covariance statistic is within
$[-1,1]$. We use v0.24.1 of scikit-learn [2] to compute the covariance matrix
via a number of estimators: the standard maximum likelihood estimate (MLE)
(denoted by the modifier EmpiricalCovariance); the elliptic envelope
(EllipticEnvelope) [3] and the minimum covariance determinant (MinCovDet)
methods for outlier removal; the Lasso technique, which uses an
$\ell_{1}$-regularization to sparsify the covariance matrix GraphicalLasso (as
well as a method with the regularization method chosen through cross
validation with five splits, GraphicalLassoCV); a basic shrinkage covariance
estimator with a fixed shrinkage coefficient of 0.1 (ShrunkCovariance) [2];
the Ledoit-Wolf method for optimizing the shrinkage coefficient (LedoitWolf)
[4]; and the oracle approximating shrinkage, an improved method for optimizing
the shrinkage coefficient if the data are Gaussian (OAS) [5]. By using all
estimators, we obtain 8 SPIs:

  1.   2. 1.

cov_EmpiricalCovariance

  3. 2.

cov_EllipticEnvelope

  4. 3.

cov_MinCovDet

  5. 4.

cov_GraphicalLasso

  6. 5.

cov_GraphicalLassoCV

  7. 6.

cov_ShrunkCovariance

  8. 7.

cov_LedoitWolf

  9. 8.

cov_OAS

####  S1.1.2 Precision (prec)

Keywords: undirected, linear, signed, multivariate, contemporaneous.

The precision matrix is the matrix inverse of the covariance matrix, and can
be used to quantify the association between each pair of time series while
controlling for concomitant effects of all other time series. For normalized
time-series data, the precision matrix is equivalent to the partial
correlation between each pairwise time series, conditioned on all other time
series, and is within $[-1,1]$. The precision matrix is computed via the same
module as the covariance matrix (in scikit-learn [2]), and has the same
estimators. By using all estimators, we obtain 8 SPIs:

  1.   2. 9.

prec_EmpiricalCovariance

  3. 10.

prec_EllipticEnvelope

  4. 11.

prec_MinCovDet

  5. 12.

prec_GraphicalLasso

  6. 13.

prec_GraphicalLassoCV

  7. 14.

prec_LedoitWolf

  8. 15.

prec_OAS

  9. 16.

prec_ShrunkCovariance

####  S1.1.3 Spearman’s rank-correlation coefficient (spearmanr)

Keywords: undirected, nonlinear, signed, bivariate, contemporaneous.

Spearman’s $\rho$ is a nonparametric measure of rank correlation between
variables. The use of ordinal (ranked) variables allows the statistic to
capture non-linear (but monotonic) relationships between random variables. The
method is implemented via function spearmanr in SciPy [6], and has a value in
$[-1,1]$.

  1. 17.

spearmanr

####  S1.1.4 Kendall’s rank-correlation coefficient (kendalltau)

Keywords: undirected, nonlinear, signed, bivariate, contemporaneous.

Kendall’s $\tau$ [1] assesses the association of ordinal variables, similar to
Spearman’s $\rho$, but has certain differences, such as becoming more
mathematically tractable in the event of ties [7]. The method is implemented
via function kendalltau in SciPy [6], and has a value in $[-1,1]$.

  1. 18.

kendalltau

####  S1.1.5 Cross correlation (xcorr)

Keywords: undirected, linear, signed/unsigned, bivariate, time-dependent.

The cross-correlation function is defined as the Pearson correlation between
two time series for all lags [8], giving values in $[-1,1]$ for each lag. To
estimate the cross-correlation function, we use SciPy [6], which outputs a
correlogram, i.e., the correlation from the MLE of the cross-covariance at a
given lag, normalised by the autocovariance. The cross-correlation is computed
with fewer observations at larger lags and so it is common to truncate the
function at a given level [9, 10], which we do by only using the first $T/4$
lags, where $T$ is the number of observations. Moreover, a correlation below
$1.96/\sqrt{T}$ is considered statistically insignificant [8], thus we
optionally cut-off the lags at this level (the modifier sig-True means we only
use the statistically significant values, and sig-False means we use all
values). We take the two summary statistics of the correlogram: the maximum
over the considered lags (denoted by modifier max), and the average over the
considered lags (mean). By using each approach outlined above, we obtain three
SPIs:

  1.   2. 19.

xcorr_max_sig-True

  3. 20.

xcorr_mean_sig-True

  4. 21.

xcorr_mean_sig-False

###  S1.2 Distance similarity (26 SPIs)

In this section, we detail SPIs that we have categorized as ‘distance-based
similarity’ measures in that they aim to establish statistical similarity or
independence based on the pairwise distance between bivariate observations.
All of the methods presented in this section are implemented using one of the
following toolboxes:

  * •

For distance correlation (and related) statistics, as well as reproducing
kernel Hilbert space (RKHS)-based statistics, we use the Hypothesis Testing in
Python (hyppo) package [11].

  * •

For dynamic time warping (and related) statistics, we use the tslearn package
[12].

####  S1.2.1 Distance correlation (dcorr)

Keywords: undirected, nonlinear, unsigned, bivariate, contemporaneous.

Distance correlation is used to infer the independence between two random
variables via pairwise distance metrics and hypothesis tests [13]. The sample
distance correlation is computed by summing over the entry-wise product of
Euclidean distance matrices. This statistic is biased, however an unbiased
estimator can be obtained by first double-centering the distance matrices
[14]. Although any pairwise distance metric can be used, here we only use
Euclidean distance because distance correlation computed with this metric has
been shown to be universally consistent (asymptotically converging to the true
value). We compute both the biased and unbiased statistics from the hyppo
package, yielding two SPIs:

  1.   2. 22.

dcorr

  3. 23.

dcorr_biased

####  S1.2.2 Cross distance correlation (dcorrx)

Keywords: undirected, nonlinear, unsigned, bivariate, time-dependent.

The cross-distance correlation [11] quantifies the independence between two
univariate time series based on distance correlation [13] (see above, Sec.
S1.2.1). This measure is the average of lagged distance correlations between
the past of time series $x$ to the future of time series $y$, up to a given
lag. We include both a low-order (lag-1, maxlag-1) and a high-order (lag-10,
maxlag-10) assumption of the number of relevant lags, yielding two SPIs:

  1.   2. 24.

dcorrx_maxlag-1

  3. 25.

dcorrx_maxlag-10

####  S1.2.3 Multiscale graph correlation (mgc)

Keywords: undirected, nonlinear, unsigned, bivariate, contemporaneous.

Multiscale graph correlation (MGC) [14] is a generalization of distance
correlation [13] (see Sec. S1.2.1) that is designed to overcome its
limitations in inferring nonlinear relationships such as circles and
parabolas. Specifically, the algorithm truncates the Euclidean distance
matrices (the procedure at the core of distance correlation) at an optimal
threshold. MGC also includes a smoothing function that is intended to remove
any bias introduced by disconnected components in the graph introduced by
truncating the distance matrix. Only this unbiased estimator is included here
(as in hyppo):

  1. 26.

mgc

####  S1.2.4 Cross multiscale graph correlation (mgcx)

Keywords: undirected, nonlinear, unsigned, bivariate, time-dependent.

The cross-multiscale graph correlation [11] is defined similarly to dcorrx
(Sec. S1.2.2) but uses lagged MGCs instead of lagged distance correlations
(see Sec. S1.2.3). By taking lag-1 (maxlag-1) and lag-10 (maxlag-10), we get
two SPIs:

  1.   2. 27.

mgcx_maxlag-1

  3. 28.

mgcx_maxlag-10

####  S1.2.5 Hilbert-Schmidt Independence Criterion (hsic)

Keywords: undirected, nonlinear, unsigned, bivariate, contemporaneous.

The Hibert-Schmidt Independence Criterion (HSIC) [15] is an RKHS-based
statistic that quantifies a statistical dependence between random variables
via a sample kernel matrix. As with distance correlation, HSIC yields a biased
statistic, where the unbiased (and consistent) estimator can be derived by
double centering the kernel distance matrix [15]. Both biased and unbiased
estimators are computed using hyppo, yielding two SPIs:

  1.   2. 29.

hsic

  3. 30.

hsic_biased

####  S1.2.6 Heller-Heller-Gorfine Independence Criterion (hhg)

Keywords: undirected, nonlinear, unsigned, bivariate, contemporaneous.

The Heller-Heller-Gorfine (HHG) method [16] yields an RKHS-based statistic
that uses the ranks of random variables to obtain sample kernel matrices,
rather than their distances (cf. HSIC in Sec. S1.2.5). This SPI is computed
via the hyppo package:

  1. 31.

hhg

####  S1.2.7 Dynamic time warping (dtw)

Keywords: undirected, nonlinear, unsigned, bivariate, time-dependent.

Dynamic time warping (DTW) [17] extends the ideas of measuring the pairwise
Euclidean distance between time series by allowing for potentially dilated
time series of variable size. Specifically, DTW finds the minimum distance
between two time series through alignment (shifting and dilating of the
sequences). This algorithm (and many outlined below) also include the Sakoe-
Chiba band [17] and the Itakura parallelogram [18] global constraints on the
alignments to prevent pathological warpings. We compute this statistic using
the tslearn package for the three global constraints (i.e., no constraints (no
modifier), the band (modifier sakoe_chiba), and the parallelogram (itakura)),
yielding three SPIs:

  1.   2. 32.

dtw

  3. 33.

dtw_constraint-itakura

  4. 34.

dtw_constraint-sakoe_chiba

####  S1.2.8 Longest common subsequence (lcss)

Keywords: undirected, nonlinear, unsigned, bivariate, time-dependent.

The longest common subsequence [19] generalizes ideas from DTW by measuring
the similarity between continuous subsections of time series, rather than the
entire time series themselves, subject to distance thresholds and alignment
constraints. We use the default threshold from the tslearn package
($\epsilon=1$), and include the three global constraints discussed above:

  1.   2. 35.

lcss

  3. 36.

lcss_constraint-itakura

  4. 37.

lcss_constraint-sakoe_chiba

####  S1.2.9 Soft dynamic time warping (softdtw)

Keywords: undirected, nonlinear, unsigned, bivariate, time-dependent.

Soft dynamic time warping [20] uses a smoothed formulation of DTW to optimize
the minimal-cost alignment as a differentiable loss function. This method is
computed via the tslearn package (with the default hyperparameter,
$\gamma=1$), and includes the same global constraints as DTW:

  1.   2. 38.

softdtw

  3. 39.

softdtw_constraint-itakura

  4. 40.

softdtw_constraint-sakoe_chiba

####  S1.2.10 Barycenter (bary)

Keywords: undirected, nonlinear, unsigned, bivariate, time-dependent.

A barycenter (or Fréchet mean) [21, 22] is a (univariate) time series that
minimizes the sum of squared distances between MTS. The tslearn package
provides functions to obtain barycenters by minimizing the sum-of-squares of
the following distance metrics: (unwarped) Euclidean distance (euclidean);
alignment via expectation maximization (dtw); alignment via a differential
loss function (softdtw); and alignment via a subgradient descent algorithm
(sgddtw). For each pair of (bivariate) time series in an $M$-variate MTS, we
compute the raw barycenters then summarize them by taking their mean (modifier
mean) and maximum (max), yielding seven SPIs:

  1.   2. 41.

bary_euclidean_max

  3. 42.

bary_dtw_mean

  4. 43.

bary_dtw_max

  5. 44.

bary_softdtw_mean

  6. 45.

bary_softdtw_max

  7. 46.

bary_sgddtw_mean

  8. 47.

bary_sgddtw_max

###  S1.3 Causal inference (10 SPIs)

In this section, we detail the 10 SPIs that we have categorized as ‘causal
inference’. These statistics aim to establish directed independence from
bivariate observations, typically by making assumptions about the underlying
model. We use the following two packages:

  * •

For convergent cross-mapping, we use the Empirical Dynamic Modeling (pyEDM)
package [23].

  * •

For all other SPIs, we use v0.5.23 of the Causal Discovery Toolbox (cdt) [24].

####  S1.3.1 Additive noise model (anm)

Keywords: directed, nonlinear, unsigned, bivariate, contemporaneous.

Additive noise models [25] are used for hypothesis testing directed nonlinear
dependence (or causality) of $x\to y$ by making the assumption that the effect
variable, $y$, is a function of a cause variable, $x$, plus a noise term (that
is independent of the cause). In this framework we use the statistic from cdt
as our SPI, which is computed by first predicting $y$ from $x$ via a Gaussian
process (with a radial basis function kernel), and then computing the
normalized HSIC test statistic from the residuals (see Sec. S1.2.5 above)
yielding the SPI:

  1. 48.

anm

####  S1.3.2 Information-geometric causal inference (igci)

Keywords: directed, nonlinear, unsigned, bivariate, contemporaneous.

Information-geometric causal inference [26] is a method for inferring causal
influence from $x\to y$ for deterministic systems with invertible functions.
The statistic is computed using cdt as the difference in differential
entropies where the probability density is computed via nearest-neighbor
estimators.

  1. 49.

igci

####  S1.3.3 Conditional distribution similarity fit (cds)

Keywords: directed, nonlinear, unsigned, bivariate, contemporaneous.

The conditional distribution similarity fit [27] is the standard deviation of
the conditional probability distribution of $y$ given $x$, where the
distributions are estimated by discretizing the values.

  1. 50.

cds

####  S1.3.4 Regression error-based causal inference (reci)

Keywords: directed, nonlinear, unsigned, bivariate, contemporaneous.

The regression error-based causal inference method [28] is an estimate of the
causal effect of $x\to y$ by quantifying the error in a regression of $y$ on
$x$ with a monomial (power product) model. In the bivariate case, this
statistic is the MSE of the linear regression of the cubic (plus constant) of
$x$ with $y$, giving one SPI:

  1. 51.

reci

####  S1.3.5 Convergent cross-mapping (ccm)

Keywords: directed, nonlinear, unsigned, bivariate, time-dependent.

The idea behind convergent cross-mapping (CCM) [29] is that there is a causal
influence from time series $x\to y$ if the Takens time-delay embedding [30] of
$y$ can be used to predict the observations of $x$. The algorithm quantifies
the prediction error (in terms of Pearson’s $\rho$) of time series $x$ from
the delay embedding of time series $y$ for increasing library sizes (i.e.,
time series length being used in the predictions). If, as the library size
increases, the correlation converges and is higher in one direction than the
other, there is an inferred causal link. The results of CCM are typically
represented as two curves [29]; one for each causal direction ($x\to y$ and
$y\to x$) with the library size on the horizontal axis and the prediction
quality (correlation) on the vertical axis.

We use the pyEDM package [23] to compute CCM, which requires an embedding
dimension to be set for the delay embedding of each time series. We use both
fixed embedding dimensions (with dimension 1 and 10, indicated by modifiers
E-1 and E-10, respectively) and an inferred embedding dimension from
univariate phase-space reconstruction methods (modifier E-None). Following the
Supplementary Materials of the original CCM paper [29] and the documentation
in the pyEDM package, we infer the embedding dimension to be the maximum of
the two univariate delay embeddings that best predicted each time series.
Given a fixed or inferred embedding dimension, we have an upper and lower
bound on the minimum and maximum library size that can be used for computing
CCM. In this work we use 21 uniformly sampled library sizes between this
minimum and maximum to generate the CCM curves. Once the curve (prediction
quality as a function of library size) is obtained, we take summary statistics
of the mean (modifier mean) and the maximum (max) across the curves. We do not
explicitly measure convergence of the algorithm as a function of library size,
consistent with common practice in the literature (e.g., [31]), but note that
this differs from the original theory; no automatic algorithm or heuristic was
originally proposed for quantifying convergence [29]. The combination of
statistics and embedding dimension parameters yields six SPIs:

  1.   2. 52.

ccm_E-1_mean

  3. 53.

ccm_E-1_max

  4. 54.

ccm_E-10_mean

  5. 55.

ccm_E-10_max

  6. 56.

ccm_E-None_mean

  7. 57.

ccm_E-None_max

###  S1.4 Information theory (37 SPIs)

The pairwise measures that we employ from information theory are either
intended to operate on serially independent observations (e.g., joint entropy
and mutual information) or bivariate time series (e.g., transfer entropy and
stochastic interaction). We primarily use v1.5 of the Java Information
Dynamics Toolkit [32] in this section, which allows us to compute differential
entropy, mutual information, and transfer entropy, in order to construct many
information-theoretic measures. A density estimation is required to compute
information-theoretic measures [33], and in this work we use four different
estimators (see references in [32]):

  * •

The Gaussian-distribution model (denoted by modifier gaussian) assumes a
linear-Gaussian multivariate, where the measure is derived from the cross-
covariance matrix;

  * •

Kernel estimation (kernel) uses a box kernel method with a specific kernel
width (default width of 0.5 standard deviations (denoted by the modifier
W-0.5));

  * •

The Kozachenko-Leonenko technique (kozachenko) is a nearest-neighbour approach
that is suitable when measures can only be constructed using entropy or joint
entropy (not mutual information); or

  * •

The Kraskov-Stögbauer-Grassberger (KSG) technique (ksg) combines nearest-
neighbor estimators for mutual information based measures (default of four
nearest neighbors, indicated by modifier NN-4). The KSG estimator is
effectively a combination of multiple Kozachenko estimators that includes
techniques to remove a bias that is introduced by taking the difference
between two differential entropy estimates.

An alternative approach to using the continuous estimators above would be to
discretise each observation (e.g., through binning) and use a discrete
estimator; however, discrete estimators are known to be heavily dependent on
the discretization size [34].

Density estimates for mutual information (and related measures) are biased by
the autocorrelation present in individuals signals [10, 35]. A common solution
to reduce bias in nonlinear estimators (like the KSG technique) is to use
dynamic correlation exclusion (also known as a Theiler window) [36, 35], which
excludes any data points within a given time window from the density estimate
and avoids oversampling. The window size should be large enough to render
observations included in the density estimate uncorrelated (sometimes called
the “autocorrelation time” [36]); here, we set the window as the product of
the autocorrelation functions of both time series (a heuristic for the
autocorrelation time based on Bartlett’s formula [10, 37]). The DCE modifier
indicates that the Theiler window is used for the KSG estimator.

####  S1.4.1 Joint entropy (je)

Keywords: undirected, nonlinear, unsigned, bivariate, contemporaneous.

The joint entropy [33] quantifies the uncertainty over the paired
observations:

  1.   2. 58.

je_gaussian

  3. 59.

je_kozachenko

  4. 60.

je_kernel_W-0.5

####  S1.4.2 Conditional entropy (ce)

Keywords: undirected, nonlinear, unsigned, bivariate, contemporaneous.

Conditional entropy [33] quantifies the uncertainty over the observations in
$y$ in the context of simultaneously observing $x$:

  1.   2. 61.

ce_gaussian

  3. 62.

ce_kozachenko

  4. 63.

ce_kernel_W-0.5

####  S1.4.3 Mutual information (mi)

Keywords: undirected, nonlinear, unsigned, bivariate, contemporaneous.

Mutual information [33] is an undirected measure of the (potentially
nonlinear) dependence between paired observations of $x$ and $y$.

  1.   2. 64.

mi_gaussian

  3. 65.

mi_kraskov_NN-4

  4. 66.

mi_kraskov_NN-4_DCE

  5. 67.

mi_kernel_W-0.25

####  S1.4.4 Time-lagged mutual information (tlmi)

Keywords: undirected, nonlinear, unsigned, bivariate, time-dependent.

Time-lagged mutual information [35] is an undirected measure of the dependence
between time series $x$ and a time-lagged instance of the series $y$. We
include statistics for only lag-one mutual information, giving four
statistics:

  1.   2. 68.

tlmi_gaussian

  3. 69.

tlmi_kraskov_NN-4

  4. 70.

tlmi_kraskov_NN-4_DCE

  5. 71.

tlmi_kernel_W-0.25

####  S1.4.5 Transfer entropy (te)

Keywords: directed, nonlinear, unsigned, bivariate, time-dependent.

Transfer entropy [35] is a measure of information transfer from a source time
series $x$ to a target time series $y$, based on the Takens time-delay
embedding. Delay embeddings capture the relevant history of a time series that
can be used as a predictor of its future and are constructed from an embedding
length and a time delay. The embedding lengths are denoted by modifier’s l and
k for the source, $x$, and target, $y$, respectively; the time delay is
denoted lt and kt for time series $x$ and $y$. The embedding parameters can be
obtained in a number of ways, so long as their product is (significantly) less
than the number of observations. We compute transfer entropy for both fixed
and optimized embedding parameters. The fixed parameters are minimal values,
i.e., typically with a fixed target embedding length of 1 (denoted by k-1) or
2 (denoted by k-2). The optimized parameters are inferred by choosing the
embedding parameters (up to a maximum embedding length of 10, denoted by
k-max-10, and maximum time delay of 2, tau-max-2) that maximize a univariate
information-theoretic measure known as active information storage (see the
MAX_CORR_AIS method of [32] for details). Note that there is no inclusion of
the gaussian estimator for transfer entropy, since this is equivalent to
Granger causality (see below) [38]. We also include a symbolic estimator for
transfer entropy (denoted by the modifier symbolic), with details in [39].

  1.   2. 72.

te_kraskov_NN-4_k-max-10_tau-max-4

  3. 73.

te_kraskov_NN-4_DCE_k-max-10_tau-max-4

  4. 74.

te_kraskov_NN-4_DCE_k-1_kt-1_l-1_lt-1

  5. 75.

te_kraskov_NN-4_DCE_k-2_kt-1_l-1_lt-1

  6. 76.

te_kraskov_NN-4_k-1_kt-1_l-1_lt-1

  7. 77.

te_kernel_W-0.25_k-1

  8. 78.

te_symbolic_k-1_kt-1_l-1_lt-1

  9. 79.

te_symbolic_k-10_kt-1_l-1_lt-1

####  S1.4.6 Granger causality (gc)

Keywords: directed, linear, unsigned, bivariate, time-dependent.

Granger causality [40] is obtained by assessing directed dependence of $x\to
y$ as the predictive power of a bivariate autoregressive model (comprising $x$
and $y$) over the univariate autoregressive model (with $y$ only). The
statistic included in our framework is the $\log$-ratio of residual variance
for the two models, which is a variant of Granger causality popularized by
Geweke [41] later found to be equivalent to transfer entropy with a Gaussian
estimator [38]. We compute the time delay and embedding dimension in the same
way as transfer entropy (see Sec. S1.4.5), allowing for both fixed and
optimized embedding lengths (k and l) and a time delays (kt and lt), where the
optimization procedure is identical to transfer entropy. This gives two SPIs:

  1.   2. 80.

gc_gaussian_k-max-10_tau-max-2

  3. 81.

gc_gaussian_k-1_kt-1_l-1_lt-1

####  S1.4.7 Causally conditioned entropy (cce)

Keywords: directed, nonlinear, unsigned, bivariate, time-dependent.

Causally conditioned entropy [42, 43, 44] aims to measure the uncertainty
remaining in time series $y$ in the context of the entire (causal) past of
both time series $x$ and $y$. The measure is computed as a sum of conditional
entropies (of $y$ given both the past of $x$ and the past of $y$) with
increasing history lengths. The standard assumption is that we consider the
entire past of both time series (i.e., there are $T-1$ conditional entropies
in the sum, from $1$ to $T-1$); however, for computational reasons, we
restrict the history length to $10$. This corresponds to the assumption that
the joint process is, at maximum, a 10th-order Markov chain.

  1.   2. 82.

cce_gaussian

  3. 83.

cce_kozachenko

  4. 84.

cce_kernel_W-0.5

####  S1.4.8 Directed information (di)

Keywords: directed, nonlinear, unsigned, bivariate, time-dependent.

Directed information [45, 42] is a measure of information flow from a source
time series $x$ to a target time series $y$ that is related to transfer
entropy, but has no time lag between source and target [46, App. C]. The
directed information can be computed as a difference between the conditional
entropy of $y$ given its own past, and causally conditioned entropy [42, 43,
44]. For the same reasons as for causally conditioned entropy, we restrict its
computation up to a history length of 10.

  1.   2. 85.

di_gaussian

  3. 86.

di_kozachenko

  4. 87.

di_kernel_W-0.5

####  S1.4.9 Stochastic interaction (si)

Keywords: undirected, nonlinear, unsigned, bivariate, time-dependent.

Stochastic interaction [47, 48] is a measure of integrated information between
two processes in the context of their own past. It is quantified by the
difference between the joint entropy of the bivariate process and the
individual entropies of each univariate process. Both entropies are measured
in context of (conditioned on) the history of the processes; we restrict this
history to be only one-step, assuming a first-order Markov process.

  1.   2. 88.

si_gaussian_k-1

  3. 89.

si_kozachenko_k-1

  4. 90.

si_kernel_W-0.5_k-1

####  S1.4.10 Integrated information (phi)

Keywords: undirected, nonlinear, unsigned, bivariate, time-dependent.

Integrated information was proposed to capture some aspects of consciousness
as part of Integrated Information Theory (IIT) [49, 50]. We implement two
proxy measures for IIT 2.0 [51] from the PhiToolbox [52]:

  * •

$\Phi^{*}$ is a proxy of integrated information [53]. It is an undirected
measure that uses the concept of mismatched decoding in information theory
[54] and can be considered as the amount of loss of information due to the
disconnection between two variables.

  * •

$\Phi_{G}$ is a measure of integrated information derived from information
geometry [55]. It is an undirected measure that quantifies the divergence
between the actual probability distribution of a system and an approximated
probability distribution where influences among elements are statistically
disconnected. Here we implement the Gaussian-distribution model [53].

Both measures are optionally normalized (divided) by entropy (denoted by the
modifier norm), yielding four SPIs:

  1.   2. 91.

phi_star_t-1_norm-0

  3. 92.

phi_star_t-1_norm-1

  4. 93.

phi_Geo_t-1_norm-0

  5. 94.

phi_Geo_t-1_norm-1

###  S1.5 Spectral (126 SPIs)

Spectral SPIs are computed in the frequency or time-frequency domain, using
either Fourier or wavelet transformations to derive spectral matrices. Unless
otherwise stated, the frequency-domain (i.e., Fourier-based) measures are
computed using the Spectral Connectivity Toolbox [56]. Each measure is based
on the cross- and power-spectral densities at a given frequency and sampling
rate, which is estimated via the multitaper technique. All directed measures
(excluding parametric spectral Granger causality) are nonparametrically
estimated through spectral matrix factorization [57, 58]. Moreover, the phase
slope index is both computed in the frequency domain (via Fourier
transformation) and the time-frequency domain (via Morlet transformation). For
more detail on the calculation, limitations, or interpretation of each of the
following measures, see the review by Bastos and Schoffelen [59].

Most methods in this literature category return a discrete set of values
across a given frequency range. Specifically, there is one value per frequency
bin, $f$, for a given sampling frequency, $f_{s}$. In order to compare
processes with very different timescales in this work (e.g., economic time
series sampled daily, or neural activity sampled at millisecond scale), we
consider the time step, $\Delta t$ (s), between successive measurements of a
time series to be rescaled by a timescale, $t_{s}$, appropriate for the
process of interest, yielding a dimensionless time step,
$\Delta\tilde{t}=\Delta t/t_{s}$. Accordingly, we assume a sampling frequency
$f_{s}=1$ throughout (denoted by modifier fs-1). We use 125 uniformly sampled
bins across the entire frequency range, $f\in[f_{0},f_{s}/2]$, where $f_{s}/2$
is the Nyquist frequency and $f_{0}=4/T$ is chosen as the minimal frequency
for computational reasons. In order to obtain SPIs, we take the mean (denoted
by modifier mean) and the maximum (max) values of each measure over three
ranges of the spectrum with corresponding modifiers, outlined below.

  * •

fmin-0_fmax-0-5: the full-frequency range, $f\in[0,f_{s}/2]$;

  * •

fmin-0_fmax-0-25: the lower frequencies, $f\in[0,f_{s}/4]$; and

  * •

fmin-0-25_fmax-0-5: the higher frequencies, $f\in[f_{s}/4,f_{s}/2]$.

####  S1.5.1 Coherence magnitude (cohmag)

Keywords: undirected, linear, unsigned, bivariate, frequency-dependent.

The coherence (also known as coherence magnitude and ordinary coherence [60]
or coherence coefficient [59]) is an undirected frequency-dependent measure of
linear association between time series $x$ and $y$. Mathematically, it is the
frequency domain equivalent of the squared time-domain cross correlation [59].
We compute mean (denoted by mean) and maximum (max) summary statistics of the
coherence for the three frequency ranges described above, giving six SPIs:

  1. 95.

cohmag_multitaper_mean_fs-1_fmin-0_fmax-0-5

  2. 96.

cohmag_multitaper_mean_fs-1_fmin-0_fmax-0-25

  3. 97.

cohmag_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  4. 98.

cohmag_multitaper_max_fs-1_fmin-0_fmax-0-5

  5. 99.

cohmag_multitaper_max_fs-1_fmin-0_fmax-0-25

  6. 100.

cohmag_multitaper_max_fs-1_fmin-0-25_fmax-0-5

####  S1.5.2 Coherence phase (phase)

Keywords: undirected, linear, unsigned, bivariate, frequency-dependent.

By omitting the magnitude operator of the coherence (Sec. S1.5.1), we obtain
the complex-valued coherency, where the phase-difference angle (which we refer
to as the coherence phase) has been used to infer a time-delayed dependence
between two signals [59]. We compute the coherence phase for the two summary
statistics and three frequency ranges described above, yielding six SPIs:

  1. 101.

phase_multitaper_mean_fs-1_fmin-0_fmax-0-5

  2. 102.

phase_multitaper_mean_fs-1_fmin-0_fmax-0-25

  3. 103.

phase_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  4. 104.

phase_multitaper_max_fs-1_fmin-0_fmax-0-5

  5. 105.

phase_multitaper_max_fs-1_fmin-0_fmax-0-25

  6. 106.

phase_multitaper_max_fs-1_fmin-0-25_fmax-0-5

####  S1.5.3 Group delay (gd)

Keywords: directed, linear, unsigned, bivariate, frequency-dependent.

The group delay [61] infers a directed, average time-delay between two signals
by measuring the slope of the phase differences (coherence phase, Sec. S1.5.2)
as a function of the frequency (obtained through linear regression). The slope
is only computed for statistically significant coherence values and the time
delay is obtained by a simple rescaling of the slope by $2\pi$. We output the
(rescaled) time-delay statistic for three frequency splits:

  1. 107.

gd_multitaper_delay_fs-1_fmin-0_fmax-0-5

  2. 108.

gd_multitaper_delay_fs-1_fmin-0_fmax-0-25

  3. 109.

gd_multitaper_delay_fs-1_fmin-0-25_fmax-0-5

####  S1.5.4 Phase slope index (psi)

Keywords: directed, linear/nonlinear, unsigned, bivariate, frequency-dependent
and time-frequency dependent.

The phase slope index (PSI) [62] is a directed measure of information flow
computed using the complex-valued coherency (see Sec. S1.5.2). Specifically,
the phase slope index evaluates the consistency of the changes in phase
differences across a pre-specified frequency range, weighted by the coherence.
Due to its availability in v0.23.0 of the MEG and EEG Analysis python package
(MNE) [63], this is the only spectral measure that is computed in the time-
frequency domain (denoted by modifier wavelet) in addition to the frequency
domain (denoted by multitaper). The time-frequency decomposition is given by a
Morlet wavelet, with spectral densities indexed by a time-frequency tuple,
given a sampling rate and number of cycles. To obtain a statistic in the time-
frequency domain, we take the average (denoted by the mean modifier) or the
maximum (max) PSI across all time points. We use the same pre-specified
frequency ranges as outlined at the beginning of this section (i.e., low, high
or full), giving three frequency-domain measures and six time-frequency domain
measures due to the two summary statistics, i.e., nine SPIs in total:

  1. 110.

psi_multitaper_mean_fs-1_fmin-0_fmax-0-5

  2. 111.

psi_multitaper_mean_fs-1_fmin-0_fmax-0-25

  3. 112.

psi_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  4. 113.

psi_wavelet_mean_fs-1_fmin-0_fmax-0-5_mean

  5. 114.

psi_wavelet_mean_fs-1_fmin-0_fmax-0-25_mean

  6. 115.

psi_wavelet_mean_fs-1_fmin-0-25_fmax-0-5_mean

  7. 116.

psi_wavelet_max_fs-1_fmin-0_fmax-0-5_max

  8. 117.

psi_wavelet_max_fs-1_fmin-0_fmax-0-25_max

  9. 118.

psi_wavelet_max_fs-1_fmin-0-25_fmax-0-5_max

####  S1.5.5 Imaginary coherence (icoh)

Keywords: undirected, linear, unsigned, bivariate, frequency-dependent.

The imaginary part of the complex-valued coherency (referred to as imaginary
coherence [59]) is argued as a means to obtain the coherence exclusively
caused by a time delay (i.e., by removing instantaneous interactions that are
present in the real axis) [64]. We compute the imaginary coherence for the
three typical frequency ranges and two summary statistics, giving six SPIs:

  1. 119.

icoh_multitaper_mean_fs-1_fmin-0_fmax-0-5

  2. 120.

icoh_multitaper_mean_fs-1_fmin-0_fmax-0-25

  3. 121.

icoh_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  4. 122.

icoh_multitaper_max_fs-1_fmin-0_fmax-0-5

  5. 123.

icoh_multitaper_max_fs-1_fmin-0_fmax-0-25

  6. 124.

icoh_multitaper_max_fs-1_fmin-0-25_fmax-0-5

####  S1.5.6 Phase locking value (plv)

Keywords: undirected, linear, unsigned, bivariate, frequency-dependent.

The phase locking value (PLV) [65] is computed using the same formula as
coherence, however applied to the amplitude-normalized Fourier transformed
signals (i.e., normalized by individual tapers), which is argued to make it a
more robust measure of phase synchronization than the coherence measures
presented earlier [59]. We compute the PLV for the three typical frequency
ranges and two summary statistics, giving six SPIs:

  1. 125.

plv_multitaper_mean_fs-1_fmin-0_fmax-0-5

  2. 126.

plv_multitaper_mean_fs-1_fmin-0_fmax-0-25

  3. 127.

plv_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  4. 128.

plv_multitaper_max_fs-1_fmin-0_fmax-0-5

  5. 129.

plv_multitaper_max_fs-1_fmin-0_fmax-0-25

  6. 130.

plv_multitaper_max_fs-1_fmin-0-25_fmax-0-5

####  S1.5.7 Pairwise phase consistency

Keywords: undirected, linear, unsigned, bivariate, frequency-dependent.

The pairwise phase consistency (PPC) [66] measures phase synchronization by
quantifying the distribution of the (per taper) phase differences, which means
that, unlike PLV, it is not biased by the sample size [59]. We compute the PPC
for the three typical frequency ranges and two summary statistics, giving six
SPIs:

  1. 131.

ppc_multitaper_mean_fs-1_fmin-0_fmax-0-5

  2. 132.

ppc_multitaper_mean_fs-1_fmin-0_fmax-0-25

  3. 133.

ppc_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  4. 134.

ppc_multitaper_max_fs-1_fmin-0_fmax-0-5

  5. 135.

ppc_multitaper_max_fs-1_fmin-0_fmax-0-25

  6. 136.

ppc_multitaper_max_fs-1_fmin-0-25_fmax-0-5

####  S1.5.8 Phase lag index (pli)

Keywords: undirected, linear, unsigned, bivariate, frequency-dependent.

The phase lag index (PLI) [67] measures phase synchronization by averaging the
sign of the (per taper) phase difference [59]. In addition to computing the
phase lag index (pli) for the frequency domain, we include weighted variations
that make the measure more robust against electrophysical artifacts (namely,
field spread), noise, and sample-size bias by weighting components based on
the imaginary coherence [68]: the weighted phase lag index (wpli), the
debiased squared phase lag index (dspli), and the debiased squared weighted
phase lag index (dswpli). In total, this gives 24 SPIs:

  1. 137.

pli_multitaper_mean_fs-1_fmin-0_fmax-0-5

  2. 138.

pli_multitaper_mean_fs-1_fmin-0_fmax-0-25

  3. 139.

pli_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  4. 140.

pli_multitaper_max_fs-1_fmin-0_fmax-0-5

  5. 141.

pli_multitaper_max_fs-1_fmin-0_fmax-0-25

  6. 142.

pli_multitaper_max_fs-1_fmin-0-25_fmax-0-5

  7. 143.

wpli_multitaper_mean_fs-1_fmin-0_fmax-0-5

  8. 144.

wpli_multitaper_mean_fs-1_fmin-0_fmax-0-25

  9. 145.

wpli_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  10. 146.

wpli_multitaper_max_fs-1_fmin-0_fmax-0-5

  11. 147.

wpli_multitaper_max_fs-1_fmin-0_fmax-0-25

  12. 148.

wpli_multitaper_max_fs-1_fmin-0-25_fmax-0-5

  13. 149.

dspli_multitaper_mean_fs-1_fmin-0_fmax-0-5

  14. 150.

dspli_multitaper_mean_fs-1_fmin-0_fmax-0-25

  15. 151.

dspli_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  16. 152.

dspli_multitaper_max_fs-1_fmin-0_fmax-0-5

  17. 153.

dspli_multitaper_max_fs-1_fmin-0_fmax-0-25

  18. 154.

dspli_multitaper_max_fs-1_fmin-0-25_fmax-0-5

  19. 155.

dswpli_multitaper_mean_fs-1_fmin-0_fmax-0-5

  20. 156.

dswpli_multitaper_mean_fs-1_fmin-0_fmax-0-25

  21. 157.

dswpli_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  22. 158.

dswpli_multitaper_max_fs-1_fmin-0_fmax-0-5

  23. 159.

dswpli_multitaper_max_fs-1_fmin-0_fmax-0-25

  24. 160.

dswpli_multitaper_max_fs-1_fmin-0-25_fmax-0-5

####  S1.5.9 Directed transfer function (dtf)

Keywords: directed, linear, unsigned, bivariate, frequency-dependent.

The cross-spectral density matrix can be decomposed into a noise covariance
matrix and a spectral transfer matrix, from which we obtain the directed
transfer function (DTF) [69], which quantifies the inflow from $x\to y$
(according to the transfer matrix) normalized by the total inflow from all
other signals into $y$ (the row-wise sum of the spectral transfer matrix). In
addition to the DTF, we include direct-DTF (ddtf) [70], which extends DTF by
conditioning out the influence of other signals (using the partialized cross-
spectrum), giving 12 SPIs:

  1. 161.

dtf_multitaper_mean_fs-1_fmin-0_fmax-0-5

  2. 162.

dtf_multitaper_mean_fs-1_fmin-0_fmax-0-25

  3. 163.

dtf_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  4. 164.

dtf_multitaper_max_fs-1_fmin-0_fmax-0-5

  5. 165.

dtf_multitaper_max_fs-1_fmin-0_fmax-0-25

  6. 166.

dtf_multitaper_max_fs-1_fmin-0-25_fmax-0-5

  7. 167.

ddtf_multitaper_mean_fs-1_fmin-0_fmax-0-5

  8. 168.

ddtf_multitaper_mean_fs-1_fmin-0_fmax-0-25

  9. 169.

ddtf_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  10. 170.

ddtf_multitaper_max_fs-1_fmin-0_fmax-0-5

  11. 171.

ddtf_multitaper_max_fs-1_fmin-0_fmax-0-25

  12. 172.

ddtf_multitaper_max_fs-1_fmin-0-25_fmax-0-5

####  S1.5.10 Directed coherence (dcoh)

Keywords: directed, linear, unsigned, bivariate, frequency-dependent.

The directed coherence [71] is obtained from the inflow of $x\to y$ using the
spectral transfer matrix (as per DTF, Sec. S1.5.9), normalized by their noise
covariance. We obtain six SPIs:

  1. 173.

dcoh_multitaper_mean_fs-1_fmin-0_fmax-0-5

  2. 174.

dcoh_multitaper_mean_fs-1_fmin-0_fmax-0-25

  3. 175.

dcoh_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  4. 176.

dcoh_multitaper_max_fs-1_fmin-0_fmax-0-5

  5. 177.

dcoh_multitaper_max_fs-1_fmin-0_fmax-0-25

  6. 178.

dcoh_multitaper_max_fs-1_fmin-0-25_fmax-0-5

####  S1.5.11 Partial directed coherence (pdcoh)

Keywords: directed, linear, unsigned, bivariate, frequency-dependent.

The partial directed coherence [72] from $x\to y$ is the inflow (as per DTF,
Sec. S1.5.9) normalized by the total outflow from all other signals into $y$
(the column-wise sum of the spectral transfer matrix). In addition to the
partial directed coherence, we include the generalized partial directed
coherence (gpdoh) [73], which scales the relative strength of inflow from
$x\to y$ by their noise covariance, yielding 12 SPIs:

  1. 179.

pdcoh_multitaper_mean_fs-1_fmin-0_fmax-0-5

  2. 180.

pdcoh_multitaper_mean_fs-1_fmin-0_fmax-0-25

  3. 181.

pdcoh_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  4. 182.

pdcoh_multitaper_max_fs-1_fmin-0_fmax-0-5

  5. 183.

pdcoh_multitaper_max_fs-1_fmin-0_fmax-0-25

  6. 184.

pdcoh_multitaper_max_fs-1_fmin-0-25_fmax-0-5

  7. 185.

gpdcoh_multitaper_mean_fs-1_fmin-0_fmax-0-5

  8. 186.

gpdcoh_multitaper_mean_fs-1_fmin-0_fmax-0-25

  9. 187.

gpdcoh_multitaper_mean_fs-1_fmin-0-25_fmax-0-5

  10. 188.

gpdcoh_multitaper_max_fs-1_fmin-0_fmax-0-5

  11. 189.

gpdcoh_multitaper_max_fs-1_fmin-0_fmax-0-25

  12. 190.

gpdcoh_multitaper_max_fs-1_fmin-0-25_fmax-0-5

####  S1.5.12 Spectral Granger causality (sgc)

Keywords: directed, linear, unsigned, bivariate, frequency-dependent.

Spectral Granger causality [40] is the frequency domain equivalent to Granger
causality, computed via the spectral transfer matrix and noise covariance that
are estimated using either a parametric (VAR model) approach [41] or
nonparametric (spectral factorization) approach [57, 58]. We implement the
nonparametric form from the Spectral Connectivity Toolbox [56] and the
parametric implementation from v0.9 of NiTime [74]. The VAR model is estimated
via least squares using the same order for both processes, where the order is
either inferred (denoted by modifier order-None) or fixed. We infer a VAR
model order from the Bayesian information criterion (with a maximum lag of 50,
implemented in the NiTime toolkit), as well as choosing fixed orders of one
(denoted by order-1) and 20 (order-20). The combination of parametric and
nonparametric approaches, with each autoregressive order and summary
statistic, gives 24 SPIs:

  1. 191.

sgc_nonparametric_mean_fs-1_fmin-0_fmax-0-5

  2. 192.

sgc_nonparametric_mean_fs-1_fmin-0_fmax-0-25

  3. 193.

sgc_nonparametric_mean_fs-1_fmin-0-25_fmax-0-5

  4. 194.

sgc_nonparametric_max_fs-1_fmin-0_fmax-0-5

  5. 195.

sgc_nonparametric_max_fs-1_fmin-0_fmax-0-25

  6. 196.

sgc_nonparametric_max_fs-1_fmin-0-25_fmax-0-5

  7. 197.

sgc_parametric_mean_fs-1_fmin-0_fmax-0-5_order-None

  8. 198.

sgc_parametric_mean_fs-1_fmin-0_fmax-0-25_order-None

  9. 199.

sgc_parametric_mean_fs-1_fmin-0-25_fmax-0-5_order-None

  10. 200.

sgc_parametric_mean_fs-1_fmin-0_fmax-0-5_order-1

  11. 201.

sgc_parametric_mean_fs-1_fmin-0_fmax-0-25_order-1

  12. 202.

sgc_parametric_mean_fs-1_fmin-0-25_fmax-0-5_order-1

  13. 203.

sgc_parametric_mean_fs-1_fmin-0_fmax-0-5_order-20

  14. 204.

sgc_parametric_mean_fs-1_fmin-0_fmax-0-25_order-20

  15. 205.

sgc_parametric_mean_fs-1_fmin-0-25_fmax-0-5_order-20

  16. 206.

sgc_parametric_max_fs-1_fmin-0_fmax-0-5_order-None

  17. 207.

sgc_parametric_max_fs-1_fmin-0_fmax-0-25_order-None

  18. 208.

sgc_parametric_max_fs-1_fmin-0-25_fmax-0-5_order-None

  19. 209.

sgc_parametric_max_fs-1_fmin-0_fmax-0-5_order-1

  20. 210.

sgc_parametric_max_fs-1_fmin-0_fmax-0-25_order-1

  21. 211.

sgc_parametric_max_fs-1_fmin-0-25_fmax-0-5_order-1

  22. 212.

sgc_parametric_max_fs-1_fmin-0_fmax-0-5_order-20

  23. 213.

sgc_parametric_max_fs-1_fmin-0_fmax-0-25_order-20

  24. 214.

sgc_parametric_max_fs-1_fmin-0-25_fmax-0-5_order-20

###  S1.6 Miscellaneous (17 SPIs)

A small number of methods do not fit squarely into any category listed above,
and so we place them in a ‘miscellaneous’ category. Here, we discuss the use
of linear and nonlinear model fits, for which we use scikit-learn [2],
cointegration, for which we use statsmodels [75], and envelope correlation,
for which we use MNE [63].

####  S1.6.1 Linear model fit (lmfit)

Keywords: directed, linear, unsigned, bivariate, contemporaneous.

Linear regression is commonly used for establishing independence through model
fits (e.g., see additive noise models in Sec. S1.3.1). As such, we use a
number of linear models and record the mean squared error (MSE) of a
regression of $y$ on $x$. The following models (with the default parameters)
from scikit-learn [2] are included: stochastic gradient descent regression
with a squared loss (ordinary least squares fit) function (denoted by modifier
SGDRegressor); Ridge regression, which uses $\ell_{2}$-norm regularization
(Ride); the Elastic-Net model [76, 77], which uses both $\ell_{1}$ and
$\ell_{2}$-norm regularization (ElasticNet); and the Bayesian Ridge regressor
uses a gamma distribution prior (with $\lambda_{1}=\lambda_{2}=10^{-6}$) for
the $\ell_{2}$-norm regularizor in Ridge regression (BayesianRidge). This
yields four SPIs:

  1.   2. 215.

lmfit_SGDRegressor

  3. 216.

lmfit_Ridge

  4. 217.

lmfit_ElasticNet

  5. 218.

lmfit_BayesianRidge

####  S1.6.2 Gaussian process model fit (gpfit)

Keywords: directed, nonlinear, unsigned, bivariate, contemporaneous.

Similar to the linear model fits, we also use Gaussian process [78] model fits
as a nonparametric measure of influence of $x$ on $y$. Here, we use a
combination of kernels (with parameters chosen from default settings) from the
scikit-learn package [2] and compute the MSE of their fit. The fits are
computed for the dot-product kernel with inhomogenity parameter $\sigma_{0}=1$
(DotProduct) and the radial basis function (RBF) kernel with length scale
$l=1$ (RBF). Each of these kernels are separately combined with the constant
kernel (with a constant of 1.0) and the white kernel (with a noise level of
1), giving two SPIs:

  1.   2. 219.

gpfit_DotProduct

  3. 220.

gpfit_RBF

####  S1.6.3 Cointegration (coint)

Keywords: undirected, linear, unsigned, bivariate, time-dependent.

If two time series are individually integrated but some linear combination of
them has a lower order of integration, then the series are said to be
‘cointegrated’ [79]. We implement statistics for quantifying the cointegration
of bivariate time series from two tests included in v0.12.0 of statsmodels
[75]: the Augmented Engle-Granger (AEG, modifier aeg) [80] two-step test and
the Johansen test (johansen) [81].

  * •

For the AEG test, we use a lag that is inferred via either AIC (modifier aic)
or BIC (bic), with a maximum lag of 10 and obtain a $t$-statistic (tstat) of
the unit-root test on the residuals. The time series is first detrended by
assuming either a constant (c) or a constant and linear trend (ct).

  * •

For the Johansen test, we output both the maximum eigenvalue (modifier
max_eig_stat) and the trace (trace_stat) of the vector error correction model.
Similar to the AEG test, we also assume a constant (order-0) or constant and
linear trend (order-1) and fixed autoregressive lags of 1 (ardiff-1) and 10
(ardiff-10).

The combination of all tests and parameters gives 11 SPIs:

  1. 221.

coint_johansen_max_eig_stat_order-0_ardiff-10

  2. 222.

coint_johansen_trace_stat_order-0_ardiff-10

  3. 223.

coint_johansen_max_eig_stat_order-0_ardiff-1

  4. 224.

coint_johansen_trace_stat_order-0_ardiff-1

  5. 225.

coint_johansen_max_eig_stat_order-1_ardiff-10

  6. 226.

coint_johansen_trace_stat_order-1_ardiff-10

  7. 227.

coint_johansen_max_eig_stat_order-1_ardiff-1

  8. 228.

coint_johansen_trace_stat_order-1_ardiff-1

  9. 229.

coint_aeg_tstat_trend-c_autolag-aic_maxlag-10

  10. 230.

coint_aeg_tstat_trend-ct_autolag-aic_maxlag-10

  11. 231.

coint_aeg_tstat_trend-ct_autolag-bic_maxlag-10

####  S1.6.4 Power envelope correlation (pec)

Keywords: undirected, linear, unsigned, bivariate, time-dependent.

The envelope correlation [82, 83] is the correlation between the two amplitude
envelopes of $x$ and $y$. Power envelope correlation is computed using MNE,
where we use six combinations of parameters, including whether the method was
orthogonalized (modifier orth), the envelopes were squared and logs were taken
prior to correlation (log), or the absolute of the correlation coefficient was
used (abs):

  1.   2. 232.

pec

  3. 233.

pec_orth

  4. 234.

pec_log

  5. 235.

pec_orth_log

  6. 236.

pec_orth_abs

  7. 237.

pec_orth_log_abs

##  S2 Library of multivariate time series

This section summarizes the MTS assembled in our database. The database (or
‘library’ of MTS) includes synthetic time series, which we have generated, and
real-world time series, which we have obtained from various publicly available
resources.

MTS are described by $M$-vector processes, $\bm{Z}=(Z_{1},\ldots,Z_{M})$,
where each $Z_{i}$ is a univariate stochastic process [79]. Here we assume a
constant sampling period, $\Delta t$, such that each univariate time series in
the tuple, $\bm{Z}$, can be represented by a discrete-time stochastic process,
i.e., $Z_{i}=(Z_{i1},Z_{i2},\ldots)$. The time series obtained from recording
$T$ real-valued observations of these processes is referred to as a sample
path (or realization) and is represented by lower-case variables, i.e.,
$\bm{z}=(z_{1},\ldots,z_{M})$, where each $z_{i}=(z_{i1},\ldots,z_{iT})$ and
$z_{it}\in\mathbb{R}$. Following conventional statistical notation, we use
$\bm{Z}$ to denote a stochastic process, with $\bm{z}$ its realization and
$\bm{z}_{0}=(z_{i0},\ldots,z_{iM})$ the initial conditions. The length of a
realization, $\bm{z}=(\bm{z}_{1},\ldots,\bm{z}_{T})$, is denoted $T$ and the
number of variables in a multivariate observation,
$\bm{z}_{1}=(z_{11},z_{21},\ldots,z_{M1})^{\prime}$, is denoted $M$.

Time-series lengths, $T$, were chosen to be large enough to obtain reasonable
statistics but small enough to allow the entire library to be processed within
a feasible time limit: all datasets have a length of greater than $T=100$
samples and less than $T=2000$ samples. The number of processes, $M$, was
chosen to give enough data for correlating the SPIs (and, again, small enough
to ensure we could process the library). That is, for a given MTS, each SPI
outputs $M(M-1)$ pairwise interactions that can be used to correlate this SPI
with every other SPI (see details in Sec. S3.1), and we need enough data to
compute such a correlation. The smallest multivariate system we consider has
$M=5$ processes (yielding $20$ pairwise interactions/samples for directed
SPIs, and $10$ for undirected SPIs) and the largest is $M=40$ (yielding $1560$
pairwise interactions/samples for directed SPIs, and $780$ for undirected
SPIs).

Finally, certain high-dimensional dynamics that are recorded from real-world
systems or generated from synthetic models can exhibit a low intrinsic
dimensionality, meaning that the dynamics can be reproduced with a smaller
number of variables than originally present in the signal. For instance, fully
synchronized oscillators appear as a set of processes that are all copies of
one another, and so there is no more information in the entire system than the
dynamics of one oscillator. Datasets with low intrinsic dimensionality do not
allow us to meaningfully capture similarities between different SPIs because
the number of effective samples in computing their similarity (as outlined
above) is significantly lower than the original number of pairwise
interactions (due to spatiotemporal correlation) [10, 37]. In analyzing the
datasets for the purposes of inferring an empirical organization, we therefore
discarded any MTS that requires less than three principal components to
explain 95% of the variance, ensuring that the effective sample size is large
enough to estimate the correlation, and that we only include systems that
exhibit complex and non-trivial dynamics.

###  S2.1 Noise models

We generated 96 datasets from uncorrelated noise (independently sampled from a
distribution) and correlated noise (e.g., Brownian motion) models. These
models simulate dynamics for each process independently, and thus involve no
interactions between processes.

####  S2.1.1 Uncorrelated noise

We generate a number of datasets using different distributions from the NumPy
package. We perform this generation procedure for each of the following
standard distributions:

  * •

Standard Normal distribution (mean $\mu=0$ and standard deviation $\sigma=1$)

  * •

Standard Cauchy distribution (location $\mu=0$ and scale $\gamma=1$).

  * •

Standard Exponential distribution (rate parameter $\gamma=1$).

  * •

Standard Gamma distribution (shape $k=1$ and scale $\theta=1$).

  * •

Standard $t$-distribution (degrees of freedom $\mu=2$).

Each dataset was generated by sampling $T=[100,500,1000,2000]$ observations
from $M=[5,10,20]$ independent distributions, giving 12 instances each (60
total).

####  S2.1.2 Correlated noise

We used the Brain Dynamics Toolbox [84] to include stochastic differential
equations for the following correlated noise models:

  * •

Ornstein–Uhlenbeck process: generated using $M$ independent simulations of the
differential equation: $dz_{it}=\theta(\mu-z_{it})\,dt+\sigma\,dW_{it}$.
Parameters chosen were $\mu=0.5$, $\sigma=0.1$, $dt=T/10$, and with randomly
generated initial conditions (with $z_{i0}\in[0,1]$) for each simulation.

  * •

Arithmetic Brownian motion: generated using $M$ independent simulations of the
differential equation: $dz_{it}=\mu\,dt+\sigma\,dW_{t}$. Parameters chosen
were $\mu=1$, $\sigma=1$, $dt=T/100$, and $\bm{z}_{0}=\bm{1}_{M}$ the vector
of ones for each simulation.

  * •

Geometric Brownian motion: generated using $M$ independent simulations of the
differential equation: $dz_{it}=\mu Z_{it}\,dt+\sigma z_{it}\,dW_{it}$.
Parameters were chosen as $\mu=0.5$, $\sigma=0.5$, $dt=T/100$, and
$\bm{z}_{0}=\bm{1}_{M}$ for each simulation.

Each dataset is generated by sampling $T=[100,500,1000]$ observations from
$M=[9,16,25]$ independent processes, giving 12 instances each (36 total).

###  S2.2 State-space models

We generated 144 datasets from state space models, also known as vector
autoregressive moving average (VARMA) models [79]. These dynamics are
governed, in general, by an autoregressive component (up to order $p$) and a
moving average component (up to order $q$):

$\bm{z}_{t}=\bm{a}_{t}+\sum_{i=1}^{p}\bm{\Phi}_{i}\bm{z}_{t-i}+\sum_{i=1}^{q}\bm{\Theta}_{i}\bm{a}_{t-i}\,,$
(1)

where the innovations $\bm{a}_{t}\sim\mathcal{N}(0,\bm{\Sigma})$ are normally
distributed. We randomly generate autoregressive matrices, $\bm{\Phi}_{i}$ and
$\bm{\Theta}_{i}$, for simulating these dynamics. Since it can be difficult to
automate the generation of high-order VARMA models, we restrict our attention
to first-order VAR/VMA. For each type of system (VAR or VMA), we have eight
configurations, given by varying each parameter below:

  * •

Coupling density: either independent or randomly sampled coupling with a
density of 10%.

  * •

Autocorrelation: either low autocorrelation (coefficients of the matrices,
$\bm{\Phi}$ and $\bm{\Theta}$, have magnitude in $[0,0.5]$) or high
autocorrelation (coefficients have magnitude in $[0.5,1]$). The sign of the
coefficients is randomly chosen to be positive or negative.

  * •

Symmetry: either symmetric or asymmetric coefficient matrices.

Each set of these configurations are generated with $T=[100,500,1000,2000]$
observations and $M=[5,10,20]$ processes. Altogether, this gives 72
simulations from VAR models and 72 simulations from VMA models.

###  S2.3 Distributed dynamical systems

In general, we consider three types of dynamical systems (with continuous
states): coupled maps (discrete space and time), system of ordinary
differential equations (discrete space, continuous time), and partial
differential equations (continuous space and time).

####  S2.3.1 Coupled map lattice

We generated 108 datasets from coupled map lattices, which are a popular model
for studying spatiotemporal chaos. Here, we use the model presented in [85],
where each time series is generated from the following equation

$z_{it}=(1-\epsilon)f(z_{it})+\frac{\epsilon}{2}\left[f(z_{i+1,t})+f(z_{i-1,t})\right]\,,$
(2)

applied to the map, $f(z)=1-\alpha z^{2}$. We generate a number of well-
studied parameter choices with periodic boundary conditions (see [85] for
details), where the initial state is sampled from a uniform distribution,
$z_{i0}\in[-0.5,1]$. The qualitative classes of those parameters are listed
below:

  * •

Frozen Chaos: generated with $\alpha=1.45$, $\epsilon=0.2$ (with a transient
period of 10,000 time points excluded).

  * •

Pattern Selection: generated with $\alpha=1.71$, $\epsilon=0.4$ (with a
transient period of 10,000 time points excluded).

  * •

Chaotic Brownian Motion of a Defect: generated with $\alpha=1.85$,
$\epsilon=0.1$ (with a transient period of 10,000 time points excluded).

  * •

Defect Turbulence: generated with $\alpha=1.895$, $\epsilon=0.1$ (with a
transient period of 10,000 time points excluded).

  * •

Spatiotemporal Intermittency (Type I): generated with $\alpha=1.7522$,
$\epsilon=0.00115$, observed every 12 time steps (with no transient period).

  * •

Spatiotemporal Intermittency (Type II): generated with $\alpha=1.75$,
$\epsilon=0.3$ (with no transient period).

  * •

Spatiotemporal Chaos: generated with $\alpha=2.0$, $\epsilon=0.3$ (with a
transient period of 1,000 time points excluded).

  * •

Traveling wave: generated with $\alpha=1.67$, $\epsilon=0.5$, observed every
2,000 time points (with a transient period of 10,000 time points excluded).

  * •

Chaotic traveling wave: generated with $\alpha=1.69$, $\epsilon=0.5$, observed
every 5,000 time points (with a transient period of 10,000 time points
excluded).

Each configuration was simulated with $M=[5,10,20]$ and
$T=[100,500,1000,2000]$, and low-dimensional datasets were removed.

####  S2.3.2 Systems of ordinary differential equations (ODEs)

Coupled oscillators are a type of coupled (nonlinear) ordinary differential
equation (ODE) and popular distributed dynamical systems model exhibiting
spaiotemporal phenomena, like synchrony and chaos. Unless stated otherwise,
the ODEs were simulated with $M=[5,10,20]$ processes and
$T=[100,500,1000,2000]$ observations. Note that some of these datasets were
considered low-dimensional (according to our approximation of their intrinsic
dimension) and so were not included in the final analysis.

  * •

The Kuramoto model [86] is generated using the pyclustering library [87],
which simulates the following equation:

$\frac{d\theta_{i}}{dt}=\omega_{i}+\frac{K}{M}\sum_{j=1}^{M}a_{ij}\sin(\theta_{j}-\theta_{i}),$
(3)

where $\theta_{i}$ is the phase of oscillator $i$ at time $t$, $\omega_{i}$ is
its natural frequency, $K$ is the coupling constant, and $a_{ij}$ is a boolean
indicating whether oscillator $i$ is connected to oscillator $j$. In our
library, we convert the phase, $\theta_{i}$, into a magnitude, i.e.,
$z_{i}=\sin{\theta_{i}}$, to avoid discontinuities. We use three different
coupling schemes for $a_{ij}$: ‘all-to-all’, where every oscillator is
connected to every other oscillator; ‘bidirectional list’, where each
oscillator is connected with two neighbors in a list format; and ‘grid four’,
where each oscillator is connected to four of its neighbors in a grid format.
For each coupling scheme, we use six different coupling strengths:
$K=[-20,-4,-1,1,4,20]$, which affect the synchronization of the network (see
[86] for details). The combination of configurations, network size ($M$), and
observation length ($T$), yielded 216 datasets, of which 130 were sufficiently
high-dimensional to include in our analysis.

  * •

The Kuramoto-Sakaguchi model [86] model extends the Kuramoto model to include
phase shifts:

$\frac{d\theta_{i}}{dt}=\omega_{i}+\frac{K}{M}\sum_{j=1}^{M}a_{ij}\sin(\theta_{j}-\theta_{i}-A_{ij})\,,$
(4)

where $A_{ij}$ is the phase frustration. Similar to the Kuramoto model, we use
the filter $z_{i}=\sin{\theta_{i}}$, and three different coupling schemes: an
all-to-all network; a circular network; and a random network. The phase
frustration is set to the network adjacency (i.e., $A_{ij}=a_{ij}$). These
models were simulated with $M=[9,16,25]$ processes and $T=[100,500,1000]$
observations via the Brain Dynamics Toolkit [84].

  * •

The Stuart–Landau and Kuramoto model [86] combine the Stuart–Landau
equations—representing a nonlinear oscillator near the Hopf bifurcation—with
the Kuramoto model for synchronization. We use the pyclustering [87]
implementation, which simulates the equation:

$\frac{d\theta_{i}}{dt}=(i\omega_{i}+\rho^{2}_{i}-|\theta_{i}|^{2})\theta_{i}+\frac{K}{M}\sum_{j=1}^{N}a_{ij}(\theta_{j}-\theta_{i})\,,$
(5)

where $\rho_{i}$ is the radius, and $i\omega_{i}$ is the imaginary unit, $i$,
multiplied by the natural frequency of that oscillator, $\omega_{i}$. We use
the following (six) configurations:

    * –

Two-way split radius (half the oscillators have $\rho_{i}=1$ and the other
half have $\rho_{i}=2$), fixed natural frequency $\omega_{i}=1$, with coupling
strength $K=1$, and an all-to-all network structure.

    * –

Random radius $\rho_{i}\in[0,1]$, $K=1$, fixed natural frequency
$\omega_{i}=1$, and an all-to-all network structure.

    * –

Fixed radius $\rho_{i}=1$ for all oscillators, $K=1$, random natural frequency
$\omega_{i}\in[0,1]$, all-to-all network structure.

    * –

Fixed radius $\rho_{i}=1$, $K=1$, fixed $\omega_{i}=1$, with a bidirectional
list structure.

    * –

Fixed radius $\rho_{i}=1$, $K=1$, two-way split natural frequency (half with
$\omega_{i}=1$, half with $\omega_{i}=2$), and an all-to-all network
structure.

    * –

Random radius $\rho_{i}\in[0,1]$, random natural frequency
$\omega_{i}\in[0,1]$, $K=0.1$, and an all-to-all network structure.

  * •

The Hodgkin–Huxley oscillatory network [88] is a collection of interacting
artificial neurons where the central element is based on the Hodgkin–Huxley
model and the peripheral neurons take in some stimulus. We include four
configurations that have been simulated using the pyclustering toolbox [87],
where we change only a small number of parameters that affect the
synchronization of the peripheral neurons and central elements. We use the
following (four) configurations:

    * –

Two-way split stimulus with splits 25 and 47. We use the default parameters
(see [87]), except with w1, w2, and w3 set to zero (intended to keep the
oscillators unsynchonized).

    * –

Two-way split stimulus with splits 25 and 27. We use the default parameters,
except with w1 set to 0.1, and w2/w3 set to zero (intended to partially
synchronize the oscillators).

    * –

Two-way split stimulus with splits 25 and 47. We use the default parameters,
except with w1 set to 0.1, w2 set to 5.0, and w3 set to zero (intended to
synchronize the oscillators).

    * –

Three-way split stimulus with splits 25, 47, and 67. We use the default
parameters, except with deltah set to 400.

  * •

A Morris–Lecar model is an approximation to the original Hodgin–Huxley model
for neuronal dynamics. We simulate a network of these neurons, with
$M=[9,16,25]$ processes and $T=[100,500,1000]$ observations, via the Brain
Dynamics Toolkit [84, 89] (with default toolkit settings).

  * •

The local-excitatory, global-inhibitory oscillatory network (LEGION) [90, 91]
is a class of synchronization models where each oscillator corresponds to a
standard relaxation oscillator with two time scales. We use the pyclustering
implementation [87], where the relaxation oscillator is based on the van der
Pol model. LEGION takes as input a number of parameters relating to
biophysical processes, a stimulus for the neural network, and a network
structure. We use the following six configurations:

    * –

Stimulus of 1 for all oscillators, bidirectional list network structure, and
otherwise default parameters.

    * –

Three-way uniform-split stimulus (one-third with stimulus 0, next third with
stimulus 1, remainder with stimulus 0), bidirectional list network structure,
and otherwise default parameters.

    * –

Three-way uniform-split stimulus, bidirectional list network structure,
parameter Wt set to 4, parameter fi set to 10, and otherwise default
parameters.

    * –

Three-way random-split stimulus (same as uniform-split stimulus but the split
locations are chosen randomly), bidirectional list network structure,
parameter Wt set to 4, parameter fi set to 10, and otherwise default
parameters.

    * –

Three-way random-split stimulus, bidirectional list network structure,
parameter Wt set to 4, parameter fi set to 0.8, and otherwise default
parameters.

    * –

Three-way random-split stimulus, bidirectional list network structure,
parameter teta_x set to -1.1, and otherwise default parameters.

  * •

The hysteresis neural network [92], implemented in the pyclustering toolkit,
comprises a network of neurons defined according to the following equations:

$\displaystyle\frac{dz_{i}}{dt}$
$\displaystyle=-z_{i}+\sum_{j=1}^{M}a_{ij}u_{j}\,,$ (6) $\displaystyle u_{j}$
$\displaystyle=h(z_{j})=\begin{cases}+1&\quad\textrm{for }z_{j}>-1\,,\\\
-1&\quad\textrm{for }z_{j}<+1\,,\end{cases}$ (7)

where $a_{ij}$ is the coupling between oscillators $i$ and $j$, $u_{i}$ are
outputs, and $h(z_{j})$ is a bipolar piecewise linear hysteresis. We generate
dynamics using five configurations, each of which have an all-to-all network
configuration and a fixed self-coupling $a_{ii}=-4$, while varying the
neighbor coupling ($a_{ij}$ where $j\neq i$), the initial states $\bm{z}_{0}$,
and the initial outputs $\bm{u}_{0}$:

    * –

$a_{ij}=1$, a two-way split for the initial states (half of $\bm{z}_{0}$ set
to 1, remainder set to 0), and all-positive outputs, $u_{i0}=1$ for all $i$.

    * –

$a_{ij}=-1$, a two-way split for the initial states (half of $\bm{z}_{0}$ set
to 1, remainder set to 0), and all-positive outputs, $u_{i0}=1$ for all $i$.

    * –

$a_{ij}=-1$, random initial states, $\bm{z}_{0}\in[0,1]$, all-negative
outputs, $u_{i0}=-1$ for all $i$.

    * –

$a_{ij}=1$, linearly increasing initial states, $z_{i0}=1/(M-i)$, all-positive
outputs, $u_{i0}=1$ for all $i$.

    * –

$a_{ij}=-1$, linearly increasing initial states, $z_{i0}=1/(M-i)$, all-
negative outputs, $u_{i0}=-1$ for all $i$.

  * •

A system of $M$ coupled van der Pol oscillators is described by the equation:

$\displaystyle\frac{dz_{i}}{dt}$ $\displaystyle=u_{i}\,,$ (8)
$\displaystyle\frac{du_{i}}{dt}$
$\displaystyle=\mu(1-z_{i}^{2})u_{i}-z_{i}-\nu\sum_{j=1}^{M}a_{ij}\,z_{i}\,,$
(9)

where $\mu$ and $\nu$ are scalar parameters, and $a_{ij}$ is the coupling
strength between oscillator $i$ and $j$ (zero if there is no relationship). We
use the Brain Dynamics Toolkit [84, 89] to simulate the dynamics of the
network of oscillators, with: (i) a ring network structure, and (ii) uniformly
random $a_{ij}\in[0,1]$.

  * •

The Wilson–Cowan model describes the mean firing rates of reciprocally coupled
excitatory and inhibitory neurons. We include simulations from a network of
such oscillators from the Brain Dynamics Toolkit [84, 89], which solves the
differential equations for excitatory, $u$, and inhibitory, $v$, neurons
according to the following equation:

$\displaystyle\frac{du_{i}}{dt}$
$\displaystyle=-u_{i}+F\left(w_{uu}u_{i}-w_{uv}v_{i}+J_{u}+\sum_{j=1}^{M}a_{ij}u_{i}\right)\,,$
(10) $\displaystyle\frac{dv_{i}}{dt}$
$\displaystyle=-v_{i}+F\left(w_{ie}u_{i}-w_{ii}v_{i}+J_{i}\right)\,,$ (11)

where $a_{ij}$ is the coupling strength, $w_{uv}$ is the weight of the
connection to $u$ from $v$, $J_{u}$ and $J_{v}$ are injection currents, and
$F(v)=1/(1+\exp(-v))$ is a sigmoid function (note we removed any variables
from Eqs (10) and (11) that were not used). For each of the simulations, we
sample the mean firing rate of the excitatory neurons only: $z_{i}=u_{i}$.
Moreover, we obtain the connectivity structure ($a_{ij}$) from the CoCoMac
database (see [84]), and use the following parameters (obtained from code
examples in [84]): $w_{uu}=10$, $w_{vv}=-1$, $w_{uv}=10$, $w_{vu}=10$,
$J_{v}=-2.5$, and $J_{u}=-8.5$.

####  S2.3.3 Partial differential equations

Finally, partial differential equations represent spatiotemporally continuous
systems that can be sampled discretely in space and time. To obtain dynamics
of this form, we simulated wave equations using the Brain Dynamics Toolkit
[84, 89]. The toolkit simulates these systems for both:

  * •

One spatial dimension using the equations,

$\frac{\partial^{2}z}{\partial t^{2}}=c^{2}\,\frac{\partial^{2}z}{\partial
u^{2}}\,,$ (12)

where $c$ is a scalar value representing the wave propagation speed and $u$ is
the spatial dimension; and

  * •

Two spatial dimensions using the equations,

$\frac{\partial^{2}z}{\partial t^{2}}=c^{2}\left(\frac{\partial^{2}z}{\partial
u^{2}}+\frac{\partial^{2}z}{\partial v^{2}}\right)\,,$ (13)

where $u$ and $v$ are the two spatial dimensions.

For each system, we use $c=10$ (the default setting) and sample uniformly
across space to obtain the MTS, $\bm{z}$. Moreover, we use the default initial
conditions, $\bm{z}_{0}$, which are proportional to a (one- or two-
dimensional) Gaussian function with standard deviation based on the number of
processes, $\sigma=M/20$. The system is then solved using $M=[9,16,25]$
processes and $T=[100,500,1000]$ observations with periodic boundary
conditions.

###  S2.4 Real-world data

We obtained various real-world datasets from publicly available sources that
cover geophysical, medical, physiological, financial, and other datatypes.
Many of these datasets were downloaded from the University of East Anglia
(UEA) multivariate time-series classification database [93] (see this article
for their references), where every task has a number of different classes, and
often many instances of each class. We attempt to minimize biasing towards
certain tasks by using only ten instances of each class when including MTS
from the UEA database.

####  S2.4.1 Medicine

  * •

The UEA ArticularyWordRecognition datasets are electromagnetic articulographs
of subjects speaking 25 words. Nine sensors were used in data collection with
a sampling rate of 200 Hz, giving $T=144$ observations across $M=9$ sensors.
We include 10 instances from each of the 25 classes of this dataset, giving a
total of 250 MTS for this datatype.

  * •

The UEA FaceDetection datasets are 306-channel MEG recordings of subjects
being shown either a face or a scrambled image. The data are recorded over 1.5
s and then down-sampled to 250 Hz and high-pass filtered at 1 Hz, giving
$T=62$ observations per 306 channels. We subsample the channels by randomly
choosing a subsequence of channels, with $M\in[10,30]$. We include 10
instances from each of the two classes of this dataset, giving a total of 20
MTS for this datatype.

  * •

The UEA FingerMovements datasets are 28-channel EEG recordings of subjects
using a keyboard, intended for the classification task of predicting whether
the subject was about to use their left or right hand. The data were
downsampled to 100 Hz, yielding $T=50$ observations for each of the $M=28$
channels. We include 10 instances from each of the two classes of this
dataset, giving a total of 20 MTS for this datatype.

  * •

The UEA HandMovementDirection datasets are 10-channel MEG recordings of
subjects whilst moving a joystick towards one of four predetermined targets on
a monitor. The trials have $T=500$ observations of brain activity resampled at
400 Hz across each of the $M=10$ channels. We include 10 instances from each
of the four classes of this dataset, giving a total of 40 MTS for this
datatype.

  * •

The UEA Heartbeat datasets are spectrograms of various locations (mainly the
aortic, pulmonic, tricuspid, and mitral areas) of normal and abnormal
heartbeats. The spectrogram is recorded over 5 s with a window size of 0.061 s
and an overlap of 70%, where each process of the MTS is a frequency band from
the spectrogram. Originally containing 61 processes representing different
frequency bands, we subsampled the dataset to have between 15–25 processes. We
include 10 instances from each of the two classes of this dataset, giving a
total of 20 MTS for this datatype.

  * •

The UEA MotorImagery datasets are recordings from an $8\times 8$ ECoG platinum
electrode grid, recorded at 1000 Hz while subjects were imagining moving
either the left small finger or the tongue. We subsampled the 64 electrodes by
randomly choosing a subsequence of electrodes, with $M\in[10,30]$. We include
10 instances from each of the two classes of this dataset, giving a total of
20 MTS for this datatype.

  * •

The SelfRegulationSCP1 dataset records the 6-channel EEG signals of a subject,
with the intention to classify if they are moving a cursor up and down a
computer screen. The signals were recorded at a sampling rate of 256 Hz for
3.5 s, resulting in $T=896$ observations per each of the $M=6$ channels. We
include 10 instances from each of the two classes of this dataset, giving a
total of 20 MTS for this datatype.

  * •

The UEA SelfRegulationSCP2 dataset records the 7-channel EEG signals of an
artificially respirated ALS patient, with the intention to classify if they
are moving a cursor up and down a computer screen The signals were recorded at
a sampling rate of 256 Hz for 4.5 s, resulting in $T=1152$ observations per
each of the $M=7$ channels. We include 10 instances from each of the two
classes of this dataset, giving a total of 20 MTS for this datatype.

  * •

We obtained resting-state fMRI data of 100 subjects from the Human Connectome
Project [94]. The fMRI data consists of blood-oxygenated level-dependent
(BOLD) signals from 333 parcels in the brain, which were minimally
preprocessed and contain 1200 observations (see details in [95]). From this
data, we randomly chose 10 subjects from which we used a subsequence of $M=10$
parcels as our MTS.

  * •

We obtained $N$-Back fMRI data of 100 subjects from the Human Connectome
Project [94], which was minimally preprocessed and parcellated to give one
333-dimensional MTS with $T=580$ observations per subject (see details in
[95]). From this data, we randomly chose 10 subjects from which we used a
subsequence of $M=10$ parcels as our MTS.

  * •

We obtained resting-state fMRI data from the UCLA Consortium for
Neuropsychiatric Phenomics on both healthy individuals (130 subjects), and
individuals with neuropsychiatric disorders, including schizophrenia (50
subjects), bipolar disorder (49 subjects), and attention deficit/hyperactivity
disorder (ADHD, 43 subjects) [96]. Each of the recordings had $T=152$
observations and, from each of the four subject types (healthy, schizophrenia,
biopolar disorder, and ADHD), we randomly selected 10 instances, and sampled
$M=15$ (consecutively numbered) brain regions from each instance.

  * •

fMRI data was obtained for 100 anesthetized wild-type mice measured at rest
[97]. The data were parcellated using the Allen Common Coordinate Framework
(CCF v3), and, similar to the human fMRI data, we randomly chose ten subjects,
from which we subsampled MTS with $M=10$ parcels and $T=2000$ observations.

####  S2.4.2 Physiology

  * •

The UEA BasicMotions dataset contains 3D accelerometer and 3D gyroscope
recordings (at 10 Hz for 10 s) of students performing four activities We
include 10 instances from each of the four classes of this dataset, giving a
total of 40 MTS for this datatype.

  * •

The UEA Cricket dataset [98] are 3D accelerometer readings (at 184 Hz) of hand
signals of cricket umpires performing twelve signals We include 10 instances
from each of the twelve classes of this dataset, giving a total of 120 MTS for
this datatype.

  * •

The UEA NATOPS data contains 3D Euclidean coordinate data from sensors
recorded at eight locations on the hands, elbows, wrists and thumbs, while the
subject was performing six different gestures We include 10 instances from
each of the six classes of this dataset, giving a total of 60 MTS for this
datatype.

  * •

The UEA RacketSports data contains 3D accelerometer and 3D gyroscope
recordings (at 10 Hz for 3 s) of students playing a forehand/backhand in
squash, or a clear/smash in badminton We include 10 instances from each of the
four classes of this dataset, giving a total of 40 MTS for this datatype.

####  S2.4.3 Finance

We obtained financial data on the daily open prices for various stocks and
foreign exchange rates using the Pandas datareader [99]. For each datatype, we
downloaded time series for a random number of days, with $T\in[500,1500]$,
between January 1 2010 to January 1 2020.

  * •

We downloaded five MTS of daily open prices from Yahoo Finance, each of which
comprised $M\in[7,15]$ randomly sampled ticker symbols that made up the
Dow–Jones Industrial Average index on August 31, 2020.

  * •

We downloaded five MTS of daily open prices from Yahoo Finance, each of which
comprised $M\in[7,15]$ randomly sampled ticker symbols that made up the SNP500
on August 31, 2020.

  * •

We downloaded five MTS of daily open prices from the St. Louis FED (FRED),
each of which comprised seven of the major foreign exchange rates: JPY/USD,
USD/EUR, USD/CNY, USD/GBP, CAD/AUD, USD/AUD, CHF/USD, BRL/USD, and DKK/USD.

####  S2.4.4 Geophysics

  * •

A number of hybrid (real-world and model-generated) MTS were downloaded from
the CauseMe challenge [100]. For these datasets, daily components were
simulated using a random CMIP5 model, with random climate variables, that
represent different climate subprocesses (thus, we refer to this dataset as a
hybrid of real-world and simulated data). Each of the datasets listed below
are named after the experiment names found on the website, with each
originally comprising $M=40$ processes and $T=250$ observations. Here, we
include 10 instances of each datatype, subsampling the processes to
$M\in[10,20]$:

    * –

CLIM is a climate-type dataset (featuring autocorrelation and time delays)
that is aggregated to a monthly time resolution with dependencies up to a
couple of months.

    * –

CLIMnoise dataset is similar to CLIM but additionally contaminated with
observational noise.

    * –

CLIMnonstat dataset is similar to CLIM but additionally features
nonstationarity.

    * –

WEATH is a weather-type dataset (featuring autocorrelation, nonlinearity, and
time delays) that has a 5-day time resolution with dependencies of up to a
couple of months.

    * –

WEATHnoise is similar to WEATH but additionally contaminated with
observational noise.

    * –

WEATHsub is similar to WEATH but is time-subsampled, i.e., only every third
time step is observed.

  * •

Ten seismograms of recent earthquakes were downloaded from the IRIS Wilber 3
System [101]. We downloaded seismic data from several (up to 25) stations that
were closest to a given earthquake epicenter, with data recorded from the time
of the earthquake until five minutes later. We randomly subsampled these
datasets so that each of the ten MTS represented aligned seismic data from
15–25 stations, with between 1000–1500 observations in each dataset.

####  S2.4.5 Others

  * •

The UEA LSST dataset is from the Photometric LSST Astronomical Time Series
Classification Challenge to classify astronomical time-series data. The
dataset contains simulations of the light curves of 14 astronomical objects
across six different filters (passbands). We include 10 instances from each of
the thirteen classes of this dataset, giving a total of 130 MTS for this
datatype.

  * •

The UEA EigenWorms dataset is intended to classify C. Elegans as either wild
types or (one of four) mutant types based on time series recordings of their
movements (which is translated via a projection onto a base shape known as
eigenworms). We truncate the dataset to have $T=3000$ observations. We include
10 instances from each of the five classes of this dataset, giving a total of
50 MTS for this datatype.

  * •

The UEA PEMS-SF data describes the occupancy rate (between 0 and 1) of
different car lanes of San Francisco bay area freeways from Jan 1, 2008 to Mar
30, 2009. Each time series represents a single day sampled every 10 minutes
(giving 144 observations) and the goal was to classify the correct day of the
week. We subsample the channels by randomly choosing a subsequence of lanes
(processes) with a random length $M\in[10,30]$. We include 10 instances from
each of the seven classes of this dataset, giving a total of 70 MTS for this
datatype.

  * •

We obtained county-level COVID-19 cases in the USA, including records from
3006 counties over $T=700$ days of the pandemic [102]. From this database, we
generated 10 instances of both the incidence (number of daily cases) and the
cumulative incidence (cumulative sum of daily cases), where, in each instance,
we randomly chose $M\in[7,15]$ counties.

##  S3 Supplementary Methods and Results

###  S3.1 The empirical similarity index

Some example MPIs are shown in Fig. S1A for SPIs based on transfer entropy
(SPI 73), coherence magnitude (SPI 95), Kendall’s $\tau$ (SPI 18), and
convergent cross-mapping (SPI 56) applied to three different MTS. Different
SPIs can yield quite different structures in their corresponding MPIs, and the
level of similarity can vary across different datasets. For example, in Fig.
S1A, Kendall’s $\tau$ and convergent cross-mapping exhibit markedly different
behavior for the currency pair but qualitatively similar behavior on the
coupled map lattice and human fMRI data.

Figure S1:  We use the empirical similarity index, $R$, to measure the
similarity between two SPIs, based on the Spearman’s correlation coefficient
between their pairwise matrices, averaged across all MTS. A. From our MTS
library, here we show three selected examples: (i) exchange rates; (ii)
coupled maps; and (iii) a functional Magnetic Resonance Imaging (fMRI)
dataset. For each dataset, we also show how pairwise dependencies between all
elements of the system can be captured as a matrix of pairwise interactions
(MPI). The MPIs are normalized and shaded in grayscale, with white being the
lowest value, and black being the highest. We show four selected examples from
our library of SPIs: (i) Kendall’s $\tau$; (ii) convergent cross-mapping;
(iii) transfer entropy; and _(iv)_ coherence magnitude. B. On a given dataset,
we quantified the similarity between each pair of SPIs as the absolute
Spearman’s correlation coefficient of (non-diagonal) elements of the resulting
pairwise matrices, as $|r_{d}|$, and averaged this quantity across all MTS to
compute our empirical similarity score, $R$. Probability distributions of
$|r_{d}|$ are shown for two SPI pairs: (i) Kendall’s $\tau$ and convergent
cross-mapping, and (ii) transfer entropy and coherence magnitude. The
resulting $R$, the average of this distribution, is shown as a vertical dashed
line. For each of these SPI pairs we also plot an example ranked scatter plot
for the fMRI dataset (yielding an $|r_{d}|$ value annotated red in the
distribution plot). C. The cumulative distribution function for the empirical
similarity index across all pairs of SPIs, with annotations for: (i) transfer
entropy and coherence magnitude ($R=0.17$); and (ii) Kendall’s $\tau$ and
convergent cross-mapping ($R=0.49$).

The 1053 MTS have an average of approximately 11.89 processes, yielding 195
112 pairwise interactions in total with which to compute the empirical
similarity index per SPI pair. Note, however, that undirected statistics have
repeated entries in the upper and lower diagonals of the MPIs, and so the
number of effective samples (in computing the correlation coefficient) is
halved, which may induce a higher variance for undirected SPIs across all
datasets when compared to directed SPIs [10].

The calculation of $r_{d}$ is depicted as ranked scatter plots in Fig. S1B for
the similarity on human fMRI data between: (i) Kendall’s $\tau$ [1] and
convergent cross-mapping (CCM) [29], $r_{d}=0.57$; and (ii) transfer entropy
[35] and coherence, $r_{d}=-0.34$. For both of these pairs of SPIs, the
distribution of $|r_{d}|$ across all 1053 MTS is also shown in Fig. S1B: the
empirical similarity index between two SPIs, $R$, is then the average of this
distribution, as $\langle|r_{d}|\rangle_{d}$. While it was important to have a
single summary of this similarity for a pair of SPIs, as the mean, it is
important to note that some pairs of SPIs have quite a wide distribution of
scores, $|r_{d}|$, across datasets, $d$, indicating that they can give highly
correlated outputs on some MTS, but not on others. An example of a wide
distribution is shown in Fig. S1B, where the distribution of Kendall’s $\tau$
and CCM is approximately uniformly distributed across MTS, ranging from some
MTS for which the two SPIs are uncorrelated, $|r_{d}|\approx 0$, to other MTS
for which there is perfect correlation, $|r_{d}|\approx 1$.

When using this approach to compare each pair of SPIs across the MTS library,
we found a wide range of similarity indices, from $R\approx 0$ to $R\approx
1$, illustrated by the cumulative distribution function in Fig. S1C. Each
index was then used to construct the dendrogram (illustrated in Fig. S2) that
gave us our findings in the main text. While we were able to use our index to
capture important empirical relationships using hierarchical clustering, it is
difficult to interpret the raw index value, $R$, in terms of a statistical
model (for, e.g., hypothesis testing). This is because of two reasons: (i)
many MPIs contain values that are not completely independent of one another
(due to higher-order dependencies), and (ii) the set of all MPIs are
themselves not completely independent (having occasionally overlapping
dynamics). As discussed in Sec. S2, the fact that the data used in computing
the index are not independent may result in fewer effective samples than
expected, which could bias the estimate of an absolute Spearman correlation
and make many SPI pairs appear to have a higher correlation than they actually
do. Nevertheless, the effective sample size of the dataset does not affect
empirical clustering.

###  S3.2 Modular versus literature categorization of SPIs

We aimed to explore whether the literature categories or data-driven modules
provided a more succinct representation of algorithm performance on the three
MTS classification datasets studied in this paper. Grouping the SPIs by the
six literature categories (of Fig. LABEL:MT-fig:schematic) revealed a wide
variability of classification performance within each category. By contrast,
as shown in Fig. LABEL:MT-fig:classificationF, grouping SPIs by the fourteen
data-driven modules (from Fig. LABEL:MT-fig:dendrogram) results in far more
characteristic performance levels within modules, suggesting our data-driven
modular representation as a useful way to capture differential performance of
SPIs on a given task. For example, on the smartwatch activity dataset, SPIs
from each literature category display a wide range of performance, including
the information-theoretic SPIs (green in Fig. S3A), which vary in accuracy
from 45% (for the undirected SPI stochastic interaction with a
Kozechenko–Leonenko estimator, si_kozachenko_k-1) to 92% (for the directed SPI
causally conditioned entropy with a Kozachenko–Leonenko estimator,
cce_kozachenko). By contrast, we identify two modules containing consistently
high-performing SPIs: M4 (transfer entropy) and M5 (parametric Granger
causality and integrated information), in Fig. S3B. These methods relate to
measuring the communication and integration of information in a
multidimensional process, adopting the Wiener–Granger theory of ‘causality’
and ‘feedback’ discussed earlier, suggesting that this way of conceptualizing
and quantifying dependencies (predicated on self-predictability) is a useful
one for analyzing these human behavioral recordings. We find qualitatively
similar results in the EEG state dataset (Figs S3C,D) and fMRI film dataset
(Figs S3E,F).

Figure S2:  A high-resolution dendrogram illustrates the empirical
organization of SPIs at a fine-grained scale. The color scheme illustrates
clusters of SPIs below a cut-off of 0.76, and leaf nodes are labeled according
to the SPI identifiers given in Sec S1.

Figure S3:  Grouping SPIs by literature category can result in greater
performance variability on a given MTS classification problem than grouping
SPIs by modules. Classification accuracy distributions over SPIs are shown as
a violin plot when grouping SPIs as literature-derived categories, see Fig.
LABEL:MT-fig:schematic (left panels A, C, E) and by modular categories, see
Fig. LABEL:MT-fig:dendrogram (right panels, B, D, F). Results are shown for
each datasets analyzed: smartwatch movement (A, B), EEG state (C, D), and fMRI
film (E, F). Classification accuracy using the combination of all SPIs is
shown as a horizontal red line.

## References

  * Kendall [1938] M. G. Kendall.  A new measure of rank correlation.  _Biometrika_ 30, 81 (1938).
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, et al.  Scikit-learn: Machine learning in python.  _J Mach Learn Res_ 12, 2825 (2011).
  * Rousseeuw and Driessen [1999] P. J. Rousseeuw and K. V. Driessen.  A fast algorithm for the minimum covariance determinant estimator.  _Technometrics_ 41, 212 (1999).
  * Ledoit and Wolf [2004] O. Ledoit and M. Wolf.  A well-conditioned estimator for large-dimensional covariance matrices.  _J Multivar Anal_ 88, 365 (2004).
  * Chen et al. [2010] Y. Chen, A. Wiesel, Y. C. Eldar, and A. O. Hero.  Shrinkage algorithms for MMSE covariance estimation.  _IEEE Trans Signal Process_ 58, 5016 (2010).
  * Virtanen et al. [2020] P. Virtanen, R. Gommers, T. E. Oliphant, et al.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nat Methods_ 17, 261 (2020).
  * Gilpin [1993] A. R. Gilpin.  Table for conversion of kendall’s tau to spearman’s rho within the context of measures of magnitude of effect for meta-analysis.  _Educ Psychol Meas_ 53, 87 (1993).
  * Brockwell and Davis [2009] P. J. Brockwell and R. A. Davis.  _Time series: theory and methods_.  Springer Science & Business Media (2009).
  * Afyouni et al. [2019] S. Afyouni, S. M. Smith, and T. E. Nichols.  Effective degrees of freedom of the Pearson’s correlation coefficient under autocorrelation.  _NeuroImage_ 199, 609 (2019).
  * Cliff et al. [2021] O. M. Cliff, L. Novelli, B. D. Fulcher, J. M. Shine, and J. T. Lizier.  Assessing the significance of directed and multivariate measures of linear dependence between time series.  _Phys Rev Res_ 3, 013145 (2021).
  * Panda et al. [2019] S. Panda, S. Palaniappan, J. Xiong, et al.  hyppo: A multivariate hypothesis testing python package.  _arXiv preprint arXiv:1907.02088_ (2019).
  * Tavenard et al. [2020] R. Tavenard, J. Faouzi, G. Vandewiele, et al.  Tslearn, a machine learning toolkit for time series data.  _J Mach Learn Res_ 21, 1 (2020).
  * Székely et al. [2007] G. J. Székely, M. L. Rizzo, N. K. Bakirov, et al.  Measuring and testing dependence by correlation of distances.  _Ann Stat_ 35, 2769 (2007).
  * Shen et al. [2020] C. Shen, C. E. Priebe, and J. T. Vogelstein.  From distance correlation to multiscale graph correlation.  _J Am Stat Assoc_ 115, 280 (2020).
  * Gretton et al. [2008] A. Gretton, K. Fukumizu, C. Teo, et al.  A kernel statistical test of independence.  In J. Platt, D. Koller, Y. Singer, and S. Roweis, eds., _Advances in Neural Information Processing Systems_ , volume 20. Curran Associates, Inc. (2008).
  * Heller et al. [2013] R. Heller, Y. Heller, and M. Gorfine.  A consistent multivariate test of association based on ranks of distances.  _Biometrika_ 100, 503 (2013).
  * Sakoe and Chiba [1978] H. Sakoe and S. Chiba.  Dynamic programming algorithm optimization for spoken word recognition.  _IEEE Trans Signal Process._ 26, 43 (1978).
  * Itakura [1975] F. Itakura.  Minimum prediction residual principle applied to speech recognition.  _IEEE Trans Acoust_ 23, 67 (1975).
  * Vlachos et al. [2002] M. Vlachos, G. Kollios, and D. Gunopulos.  Discovering similar multidimensional trajectories.  In _Proceedings 18th international conference on data engineering_ , pp. 673–684. IEEE (2002).
  * Cuturi and Blondel [2017] M. Cuturi and M. Blondel.  Soft-dtw: a differentiable loss function for time-series.  In _International Conference on Machine Learning_ , pp. 894–903. PMLR (2017).
  * Petitjean et al. [2011] F. Petitjean, A. Ketterlin, and P. Gançarski.  A global averaging method for dynamic time warping, with applications to clustering.  _Pattern Recognit_ 44, 678 (2011).
  * Schultz and Jain [2018] D. Schultz and B. Jain.  Nonsmooth analysis and subgradient methods for averaging in dynamic time warping spaces.  _Pattern Recognit_ 74, 340 (2018).
  * Park and Smith [2021] J. Park and C. Smith.  pyEDM.  https://github.com/SugiharaLab/pyEDM (2021).
  * Kalainathan and Goudet [2019] D. Kalainathan and O. Goudet.  Causal discovery toolbox: Uncover causal relationships in python.  _arXiv preprint arXiv:1903.02278_ (2019).
  * Hoyer et al. [2008] P. Hoyer, D. Janzing, J. M. Mooij, J. Peters, and B. Schölkopf.  Nonlinear causal discovery with additive noise models.  _Adv Neural Inf Process Syst_ 21, 689 (2008).
  * Daniusis et al. [2012] P. Daniusis, D. Janzing, J. Mooij, et al.  Inferring deterministic causal relations.  _arXiv preprint arXiv:1203.3475_ (2012).
  * Fonollosa [2019] J. A. Fonollosa.  Conditional distribution variability measures for causality detection.  In _Cause Effect Pairs in Machine Learning_ , pp. 339–347. Springer (2019).
  * Blöbaum et al. [2018] P. Blöbaum, D. Janzing, T. Washio, S. Shimizu, and B. Schölkopf.  Cause-effect inference by comparing regression errors.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 900–909. PMLR (2018).
  * Sugihara et al. [2012] G. Sugihara, R. May, H. Ye, et al.  Detecting causality in complex ecosystems.  _Science_ 338, 496 (2012).
  * Takens [1981] F. Takens.  Detecting strange attractors in turbulence.  In D. Rand and L.-S. Young, eds., _Dynamical Systems and Turbulence_ , volume 898 of _Lecture Notes in Mathematics_ , pp. 366–381. Springer, Berlin / Heidelberg (1981).
  * Pao et al. [2021] G. M. Pao, C. Smith, J. Park, et al.  Experimentally testable whole brain manifolds that recapitulate behavior.  _arXiv preprint arXiv:2106.10627_ (2021).
  * Lizier [2014] J. T. Lizier.  JIDT: An information-theoretic toolkit for studying the dynamics of complex systems.  _Front Robot AI_ 1, 11 (2014).
  * MacKay [2003] D. J. MacKay.  _Information Theory, Inference and Learning Algorithms_.  Cambridge University Press (2003).
  * Bossomaier et al. [2016] T. Bossomaier, L. Barnett, M. Harré, and J. T. Lizier.  Transfer entropy.  In _An introduction to transfer entropy_ , pp. 65–95. Springer (2016).
  * Schreiber [2000] T. Schreiber.  Measuring information transfer.  _Phys Rev Lett_ 85, 461 (2000).
  * Theiler [1986] J. Theiler.  Spurious dimension from correlation algorithms applied to limited time-series data.  _Phys Rev A_ 34, 2427 (1986).
  * Bartlett [1946] M. S. Bartlett.  On the theoretical specification and sampling properties of autocorrelated time-series.  _Supplement to the Journal of the Royal Statistical Society_ 8, 27 (1946).
  * Barnett et al. [2009] L. Barnett, A. B. Barrett, and A. K. Seth.  Granger causality and transfer entropy are equivalent for Gaussian variables.  _Phys Rev Lett_ 103, 238701 (2009).
  * Staniek and Lehnertz [2008] M. Staniek and K. Lehnertz.  Symbolic transfer entropy.  _Phys Rev Lett_ 100, 158101 (2008).
  * Granger [1969] C. W. Granger.  Investigating causal relations by econometric models and cross-spectral methods.  _Econometrica_ 37, 424 (1969).
  * Geweke [1982] J. Geweke.  Measurement of linear dependence and feedback between multiple time series.  _J Am Stat Assoc_ 77, 304 (1982).
  * Kramer [1998] G. Kramer.  _PhD Dissertation: Directed information for channels with feedback_.  Citeseer (1998).
  * Amblard and Michel [2011] P.-O. Amblard and O. J. J. Michel.  On directed information theory and Granger causality graphs.  _J Comput Neurosci_ 30, 7 (2011).
  * Ziebart et al. [2013] B. D. Ziebart, J. A. Bagnell, and A. K. Dey.  The principle of maximum causal entropy for estimating interacting processes.  _IEEE Trans Inf Theory_ 59, 1966 (2013).
  * Massey [1990] J. Massey.  Causality, feedback and directed information.  In _Proc. Int. Symp. Inf. Theory Applic.(ISITA-90)_ , pp. 303–305 (1990).
  * Lizier [2013] J. T. Lizier.  _The Local Information Dynamics of Distributed Computation in Complex Systems_.  Springer, Berlin / Heidelberg (2013).
  * Ay [2015] N. Ay.  Information geometry on complexity and stochastic interaction.  _Entropy_ 17, 2432 (2015).
  * Ay and Wennekers [2003] N. Ay and T. Wennekers.  Temporal infomax leads to almost deterministic dynamical systems.  _Neurocomputing_ 52, 461 (2003).
  * Tononi [2004] G. Tononi.  An information integration theory of consciousness.  _BMC Neurosci_ 5, 1 (2004).
  * Tononi et al. [2016] G. Tononi, M. Boly, M. Massimini, and C. Koch.  Integrated information theory: From consciousness to its physical substrate.  _Nat Rev Neurosci_ 17, 450 (2016).
  * Balduzzi and Tononi [2008] D. Balduzzi and G. Tononi.  Integrated information in discrete dynamical systems: motivation and theoretical framework.  _PLoS Comput Biol_ 4, e1000091 (2008).
  * Kitazono and Oizumi [2020] J. Kitazono and M. Oizumi.  Practical $\Phi$ Toolbox for Integrated Information Analysis.  https://github.com/oizumi-lab/PhiToolbox/ (2020).
  * Oizumi et al. [2016a] M. Oizumi, S.-i. Amari, T. Yanagawa, N. Fujii, and N. Tsuchiya.  Measuring integrated information from the decoding perspective.  _PLoS Comput Biol_ 12, e1004654 (2016a).
  * Merhav et al. [1994] N. Merhav, G. Kaplan, A. Lapidoth, and S. S. Shitz.  On information rates for mismatched decoders.  _IEEE Trans Inf Theory_ 40, 1953 (1994).
  * Oizumi et al. [2016b] M. Oizumi, N. Tsuchiya, and S.-i. Amari.  Unified framework for information integration based on information geometry.  _Proc Natl Acad Sci U S A_ 113, 14817 (2016b).
  * Denovellis et al. [2020] E. Denovellis, M. Myroshnychenko, and D. Ulianych.  Eden-kramer-lab/spectral_connectivity: v0.2.5.dev0 (2020).
  * Dhamala et al. [2008a] M. Dhamala, G. Rangarajan, and M. Ding.  Analyzing information flow in brain networks with nonparametric Granger causality.  _NeuroImage_ 41, 354 (2008a).
  * Dhamala et al. [2008b] M. Dhamala, G. Rangarajan, and M. Ding.  Estimating Granger causality from Fourier and wavelet transforms of time series data.  _Phys Rev Lett_ 100, 018701 (2008b).
  * Bastos and Schoffelen [2016] A. M. Bastos and J.-M. Schoffelen.  A tutorial review of functional connectivity analysis methods and their interpretational pitfalls.  _Front Syst Neurosci_ 9 (2016).
  * Schlögl and Supp [2006] A. Schlögl and G. Supp.  Analyzing event-related EEG data with multivariate autoregressive parameters.  _Prog Brain Res_ 159, 135 (2006).  Publisher: Elsevier.
  * Gotman [1983] J. Gotman.  Measurement of small time differences between EEG channels: method and application to epileptic seizure propagation.  _Electroencephalogr Clin Neurophysiol_ 56, 501 (1983).
  * Nolte et al. [2008] G. Nolte, A. Ziehe, V. V. Nikulin, et al.  Robustly estimating the flow direction of information in complex physical systems.  _Phys Rev Lett_ 100, 234101 (2008).
  * Larson et al. [2021] E. Larson, A. Gramfort, D. A. Engemann, et al.  mne-tools/mne-python: v0.23.0 (2021).
  * Nolte et al. [2004] G. Nolte, O. Bai, L. Wheaton, et al.  Identifying true brain interaction from EEG data using the imaginary part of coherency.  _Clin Neurophysiol_ 115, 2292 (2004).
  * Lachaux et al. [1999] J.-P. Lachaux, E. Rodriguez, J. Martinerie, and F. J. Varela.  Measuring phase synchrony in brain signals.  _Hum Brain Mapp_ 8, 194 (1999).
  * Vinck et al. [2010] M. Vinck, M. van Wingerden, T. Womelsdorf, P. Fries, and C. M. Pennartz.  The pairwise phase consistency: A bias-free measure of rhythmic neuronal synchronization.  _NeuroImage_ 51, 112 (2010).
  * Stam et al. [2007] C. J. Stam, G. Nolte, and A. Daffertshofer.  Phase lag index: Assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources.  _Hum Brain Mapp_ 28, 1178 (2007).
  * Vinck et al. [2011] M. Vinck, R. Oostenveld, M. Van Wingerden, F. Battaglia, and C. M. Pennartz.  An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias.  _NeuroImage_ 55, 1548 (2011).
  * Kaminski and Blinowska [1991] M. J. Kaminski and K. J. Blinowska.  A new method of the description of the information flow in the brain structures.  _Biol Cybern_ 65, 203 (1991).
  * Korzeniewska et al. [2003] A. Korzeniewska, M. Mańczak, M. Kamiński, K. J. Blinowska, and S. Kasicki.  Determination of information flow direction among brain structures by a modified directed transfer function (dDTF) method.  _J Neurosci Methods_ 125, 195 (2003).
  * Baccala et al. [1998] L. Baccala, K. Sameshima, G. Ballester, A. Do Valle, and C. Timo-Iaria.  Studying the interaction between brain structures via directed coherence and Granger causality.  _Applied Sig Process_ 5 (1998).
  * Baccalá and Sameshima [2001] L. A. Baccalá and K. Sameshima.  Partial directed coherence: a new concept in neural structure determination.  _Biol Cybern_ 84, 463 (2001).
  * Baccala et al. [2007] L. A. Baccala, K. Sameshima, and D. Y. Takahashi.  Generalized partial directed coherence.  In _2007 15th International conference on digital signal processing_ , pp. 163–166. IEEE (2007).
  * Rokem et al. [2009] A. Rokem, M. Trumpis, and F. Perez.  Nitime: Time-series analysis for neuroimaging data.  In _Proceedings of the 8th Python in Science Conference_ , pp. 68–75 (2009).
  * Seabold and Perktold [2010] S. Seabold and J. Perktold.  statsmodels: Econometric and statistical modeling with python.  In _9th Python in Science Conference_ (2010).
  * Friedman et al. [2010] J. Friedman, T. Hastie, and R. Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _J Stat Softw_ 33, 1 (2010).
  * Andersen et al. [2011] M. Andersen, J. Dahl, Z. Liu, et al.  Interior-point methods for large-scale cone programming.  In S. Sra, S. Nowozin, and S. J. Wright, eds., _Optimization for machine learning_. The MIT Press (2011).
  * Williams and Rasmussen [1995] C. Williams and C. Rasmussen.  Gaussian processes for regression.  _Advances in neural information processing systems_ 8 (1995).
  * Reinsel [2003] G. C. Reinsel.  _Elements of Multivariate Time Series Analysis_.  Springer Science & Business Media (2003).
  * Engle and Granger [1987] R. F. Engle and C. W. Granger.  Co-integration and error correction: Representation, estimation, and testing.  _Econometrica_ 55, 251 (1987).
  * Pesaran and Shin [1999] M. H. Pesaran and Y. Shin.  An autoregressive distributed lag modelling approach to cointegration analysis.  In S. S, ed., _Econometrics and Economic Theory in the 20th Century: The Ragnar Frisch Centennial Symposium_. Cambridge University Press: Cambridge (1999).
  * Hipp et al. [2012] J. F. Hipp, D. J. Hawellek, M. Corbetta, M. Siegel, and A. K. Engel.  Large-scale cortical correlation structure of spontaneous oscillatory activity.  _Nat Neurosci_ 15, 884 (2012).
  * Khan et al. [2018] S. Khan, J. A. Hashmi, F. Mamashli, et al.  Maturation trajectories of cortical resting-state networks depend on the mediating frequency band.  _NeuroImage_ 174, 57 (2018).
  * Heitmann and Breakspear [2017] S. Heitmann and M. Breakspear.  Handbook for the brain dynamics toolbox.  _Brisbane, QLD: QIMR Berghofer Medical Research Institute_ (2017).
  * Kaneko and Tsuda [2011] K. Kaneko and I. Tsuda.  _Complex systems: Chaos and Beyond_.  Springer Science & Business Media (2011).
  * Kuramoto [2003] Y. Kuramoto.  _Chemical oscillations, waves, and turbulence_.  Courier Corporation (2003).
  * Novikov [2019] A. Novikov.  Pyclustering: Data mining library.  _J Open Source Softw_ 4, 1230 (2019).
  * Chik et al. [2009] D. Chik, R. Borisyuk, and Y. Kazanovich.  Selective attention model with spiking elements.  _Neural Netw_ 22, 890 (2009).
  * Heitmann et al. [2018] S. Heitmann, M. J. Aburn, and M. Breakspear.  The brain dynamics toolbox for Matlab.  _Neurocomputing_ 315, 82 (2018).
  * Wang and Terman [1995] D. Wang and D. Terman.  Locally excitatory globally inhibitory oscillator networks.  _IEEE Trans Neural Netw_ 6, 283 (1995).
  * Wang and Terman [1997] D. Wang and D. Terman.  Image segmentation based on oscillatory correlation.  _Neural Comput_ 9, 805 (1997).
  * Zaghloul et al. [2012] M. E. Zaghloul, J. L. Meador, and R. W. Newcomb.  _Silicon implementation of pulse coded neural networks_ , volume 266\.  Springer Science & Business Media (2012).
  * Ruiz et al. [2021] A. P. Ruiz, M. Flynn, J. Large, M. Middlehurst, and A. Bagnall.  The great multivariate time series classification bake off: A review and experimental evaluation of recent algorithmic advances.  _Data Min Knowl Discov_ 35, 401 (2021).
  * Van Essen et al. [2012] D. C. Van Essen, K. Ugurbil, E. Auerbach, et al.  The Human Connectome Project: a data acquisition perspective.  _NeuroImage_ 62, 2222 (2012).
  * Shine et al. [2019] J. M. Shine, M. Breakspear, P. T. Bell, et al.  Human cognition involves the dynamic integration of neural activity and neuromodulatory systems.  _Nat Neurosci_ 22, 289 (2019).
  * Poldrack et al. [2016] R. Poldrack, E. Congdon, W. Triplett, et al.  A phenome-wide examination of neural and cognitive function.  _Sci Data_ 3, 160110 (2016).
  * Zerbi et al. [2021] V. Zerbi, M. Pagani, M. Markicevic, et al.  Brain mapping across 16 autism mouse models reveals a spectrum of functional connectivity subtypes.  _Mol Psychiatry_ pp. 1–11 (2021).
  * Ko et al. [2005] M. H. Ko, G. West, S. Venkatesh, and M. Kumar.  Online context recognition in multisensor systems using dynamic time warping.  In _2005 International Conference on Intelligent Sensors, Sensor Networks and Information Processing_ , pp. 283–288. IEEE (2005).
  * pandas development team [2020] T. pandas development team.  pandas-dev/pandas: Pandas (2020).
  * Runge et al. [2020] J. Runge, X.-A. Tibau, M. Bruhns, J. Muñoz-Marí, and G. Camps-Valls.  The causality for climate competition.  In _NeurIPS 2019 Competition and Demonstration Track_ , pp. 110–120. PMLR (2020).
  * for Seismology [2021] I. R. I. for Seismology.  IRIS Wilber 3 system.  https://ds.iris.edu/wilber3/ (2021).
  * Times [2021] T. N. Y. Times.  Coronavirus (Covid-19) data in the United States.  https://github.com/nytimes/covid-19-data (2021).