Inferring network properties from time series using transfer entropy
                    and mutual information: Validation of multivariate versus bivariate
                    approaches

Novelli, Leonardo; Lizier, Joseph T.

doi:10.1162/netn_a_00178

Abstract

Functional and effective networks inferred from time series are at the core of network neuroscience. Interpreting properties of these networks requires inferred network models to reflect key underlying structural features. However, even a few spurious links can severely distort network measures, posing a challenge for functional connectomes. We study the extent to which micro- and macroscopic properties of underlying networks can be inferred by algorithms based on mutual information and bivariate/multivariate transfer entropy. The validation is performed on two macaque connectomes and on synthetic networks with various topologies (regular lattice, small-world, random, scale-free, modular). Simulations are based on a neural mass model and on autoregressive dynamics (employing Gaussian estimators for direct comparison to functional connectivity and Granger causality). We find that multivariate transfer entropy captures key properties of all network structures for longer time series. Bivariate methods can achieve higher recall (sensitivity) for shorter time series but are unable to control false positives (lower specificity) as available data increases. This leads to overestimated clustering, small-world, and rich-club coefficients, underestimated shortest path lengths and hub centrality, and fattened degree distribution tails. Caution should therefore be used when interpreting network properties of functional connectomes obtained via correlation or pairwise statistical dependence measures, rather than more holistic (yet data-hungry) multivariate models.

Author Summary

We compare bivariate and multivariate methods for inferring networks from time series, which are generated using a neural mass model and autoregressive dynamics. We assess their ability to reproduce key properties of the underlying structural network. Validation is performed on two macaque connectomes and on synthetic networks with various topologies (regular lattice, small-world, random, scale-free, modular). Even a few spurious links can severely bias key network properties. Multivariate transfer entropy performs best on all topologies for longer time series.

Directed connectivity, Functional connectivity, Network inference, Multivariate transfer entropy, Information theory, Complex networks

INTRODUCTION

Functional and effective network inference in neuroscience typically involves preprocessing the data, defining the parcellation, extracting the time series, inferring the links in the model network, and measuring network properties, for example, to compare patients and controls or to predict phenotype (Bassett & Sporns, 2017; Fornito, Zalesky, & Bullmore, 2016). Each step in the pipeline requires making modeling and analysis choices, whose influence on the final results is the subject of ongoing research (Aquino, Fulcher, Parkes, Sabaroedin, & Fornito, 2020; Cliff, Novelli, Fulcher, Shine, & Lizier, 2020; Zalesky et al., 2016). As part of this effort, we study how the choice of different inference algorithms affects the properties of the resulting network model in comparison to the underlying structural network, and whether the ability to accurately reflect these properties changes across different underlying structural networks.

The structure (or topology) of a network can be described at multiple scales (Bassett & Sporns, 2017): from the microscopic (individual links), to the mesoscopic (modules and motifs) and the macroscopic (summary statistics, such as average shortest path length and measures of small-worldness; Rubinov & Sporns, 2010). At each scale, the structure is associated with development, aging, cognition, and neuropsychiatric diseases (Xia et al., 2020). Previous studies have assessed the performance of different network inference algorithms in identifying the structural links at the microscale (Kim, Rogers, Sun, & Bollt, 2016; Novelli, Wollstadt, Mediano, Wibral, & Lizier, 2019; Razi, Kahan, Rees, & Friston, 2015; Runge, Nowack, Kretschmer, Flaxman, & Sejdinovic, 2018; Sun, Taylor, & Bollt, 2015). The goal of this work is to extend the assessment to all scales and to a variety of topologies, across a range of related network inference algorithms. We link the performance at the microscale to the resulting network properties at the macroscale, and describe how this changes as a function of the overall topology.

We compare bivariate and multivariate approaches for inferring network models, employing statistical dependence measures based on information theory (Shannon, 1948). These approaches include functional network inference, which produces models of networks of pairwise or bivariate statistical relationships between nodes, and can either quantify undirected statistical dependence, in the case of mutual information (MI; Cover & Thomas, 2005), or directed dependence, in the case of transfer entropy (TE; Bossomaier, Barnett, Harré, & Lizier, 2016; Schreiber, 2000). These approaches also include effective network inference, which is intended to produce the simplest possible circuit models that explain the observed responses (Aertsen, Gerstein, Habib, & Palm, 1989). In this class, we evaluate the use of multivariate TE, which, in contrast to the bivariate approaches, aims to minimize spurious links and infer minimal models of the parent sets for each target node in the network.

All of these inference techniques seek to infer a network model of the relationships between the nodes in a system. Different methods capture different aspects of these relationships and don’t necessarily seek to replicate the underlying structural topology, nor do we expect them to in general (particularly in neuroimaging experiments, where aspects of the structure may be expressed more or less or not at all, depending on the cognitive task). In spite of that, in this paper we do seek to evaluate and indeed validate these methods in inferring microscopic, mesoscopic, and macroscopic features of the underlying network structure. Crucially, we perform this validation under idealized conditions—including full observability, stationarity no subsampling, and so fourth—that allow us to establish a hypothesis that effective networks should be not just complementary to the structural but converge to it under these conditions, as our available data increase. Indeed, under these idealized conditions (specifically in the absence of hidden nodes, and other simplifying assumptions, including stationarity), effective networks inferred via multivariate TE are proven to converge to the underlying structure for sufficiently long time series (Runge, 2018; Sun et al., 2015). In gaining an understanding of these multivariate effective connectivity inference algorithms, it is important to validate that they perform to that expectation where it is applicable, and investigate how that performance varies with respect to factors such as sample size, and so on. In doing so, we also address the recent call for more extensive model diversity in testing multivariate algorithms: “To avoid biased conclusions, a large number of different randomly selected connectivity structures should be tested [including link density as well as properties such as small-worldness]” (Runge, 2018, p. 13).

Outside of these idealized circumstances though, we can no longer make a clear general hypothesis on how the effective network models are expected to reflect the underlying structure, yet a successful validation gives confidence that the directed statistical relationships they represent remain accurate as an effective network model at the microscale. Furthermore, it is at least desirable for not only effective networks but also functional networks to recognize important features in the underlying network structure: to track overall regime changes in the macroscopic structure reliably, and to reflect the mesoscopic properties of distinctive nodes (or groups of nodes) in the structure. The desire for recognition of important features in the network is applicable whether the inference is made under idealized conditions or not.

This motivates our validation study, which is primarily conducted under idealized conditions as above and based on synthetic datasets involving ground truth networks of 100–200 nodes with different topologies, from regular lattice to small-world, random, scale-free, and modular. Many of these structural properties are incorporated in the macaque connectomes analyzed in the last section. At the macroscale, we measure several fundamental and widely used properties, including shortest path length, clustering coefficient, small-world coefficient, betweenness centrality, and features of the degree distributions (Rubinov & Sporns, 2010). These properties of the inferred network models are compared with those of the real underlying structural networks in order to validate and benchmark different inference algorithms in terms of their ability to capture the key properties of the underlying topologies. At the microscale, the performance is assessed in terms of precision, recall, and specificity of the inferred model in classifying the links of the underlying structural network. As above, while we do not expect all approaches to strictly capture the microscale features, these results help to explain their performance at the macroscale.

For most of our experiments, the time series of node activity on these networks are generated by vector autoregressive (VAR) dynamics, with linearly coupled nodes and Gaussian noise. Both the VAR process and the inference algorithms are described in detail in the Methods section, where we also discuss how MI and the magnitude of Pearson correlation are equivalent for stationary VAR processes. This implies that the undirected networks obtained via the bivariate MI algorithm are equivalent to the widely employed undirected functional networks obtained via correlation, extending the implications of our results beyond information-theoretic methods. Further, our results based on TE extend to Granger causality, which is equivalent to TE for stationary VAR processes (Barnett, Barrett, & Seth, 2009). Networks inferred using bivariate TE are typically referred to as directed functional networks, to emphasize the directed nature of their links. The extension to multivariate TE for effective network inference can also be viewed as an extension to multivariate Granger causality for the stationary VAR processes here.

We find that multivariate TE performs better on all network topologies at all scales, for longer time series. Bivariate methods can achieve better recall with limited amount of data (shorter time series) in some circumstances, but the precision and the ability to control false positives are neither consistent nor predictable a priori. On the other hand, thanks to recent statistical improvements, multivariate TE guarantees high specificity and precision regardless of the amount of data available, and the recall steadily increases with more data. We discuss how the mesoscopic properties of the underlying structural network—particularly the network motifs—can influence the precision and recall of the model at the microscale. In turn, we show how the performance at the microscale affects the inferred network properties at the macroscale. We observe that bivariate methods are often unable to capture the most distinctive topological features of the networks under study (including path length, clustering, degree distribution, and modularity), largely because of their inability to control false positives at the microscale.

Our final section moves beyond the validation under idealized conditions to extend the experiments to time series from a neural mass model employed on the CoCoMac connectome. Although the incorporation of complexities such as nonlinear interactions and subsampled time series do somewhat reduce the performance of the methods, the superior performance of multivariate TE aligns with the results above from the validation in idealized conditions. While further and more wide-ranging experiments are required, our experiments provide substantial evidence for the validity of multivariate TE in providing effective network models, which still retain meaningful network insights in more realistic conditions.

METHODS

Generating Dynamics on Networks

Two models are used to generate time series dynamics on networks of coupled variables. Vector autoregressive processes are employed for validation studies under idealized conditions, and a neural mass model is used on the weighted CoCoMac connectome as a final investigation going beyond these conditions. The simulation and analysis pipeline is illustrated in Figure 1.

Figure 1.

Pipeline of the network inference comparison, using the CoCoMac
connectome as an example (see Methods and Numerical Simulations in the Macaque
Connectome section for full details). With the real
adjacency matrix defining the links, synthetic time series are generated
using either a linear autoregressive (VAR) system with uniform link
weights, or a nonlinear neural mass model with realistic weights. Three
network inference algorithms (bivariate MI, bivariate TE, multivariate
TE) are then employed to analyze the time series. Bivariate methods
consider pairs of nodes independently of each other, either taking the
past states into account (bivariate TE) or not (bivariate MI). On the
other hand, multivariate TE analyzes pairs of nodes in the context of
the whole network. At the microscale, the links in the inferred networks
are classified as true/false positives/negatives and these scores are
used to compute standard performance measured (precision, recall, and
specificity). At the macroscale, the performance is instead measures
according to the ability of each algorithm to faithfully reflect network
properties (summary statistics) of the underlying structural network
serving as ground truth.

View large Download slide

Pipeline of the network inference comparison, using the CoCoMac connectome as an example (see ec1Methods and ec3Numerical Simulations in the ec2Macaque Connectome section for full details). With the real adjacency matrix defining the links, synthetic time series are generated using either a linear autoregressive (VAR) system with uniform link weights, or a nonlinear neural mass model with realistic weights. Three network inference algorithms (bivariate MI, bivariate TE, multivariate TE) are then employed to analyze the time series. Bivariate methods consider pairs of nodes independently of each other, either taking the past states into account (bivariate TE) or not (bivariate MI). On the other hand, multivariate TE analyzes pairs of nodes in the context of the whole network. At the microscale, the links in the inferred networks are classified as true/false positives/negatives and these scores are used to compute standard performance measured (precision, recall, and specificity). At the macroscale, the performance is instead measures according to the ability of each algorithm to faithfully reflect network properties (summary statistics) of the underlying structural network serving as ground truth.

Networks of linearly coupled Gaussian variables.

We simulate discrete-time, stationary, first-order vector autoregressive processes (VAR) on underlying structural networks of N nodes. A VAR process is described by the recurrence relation

Z (t + 1) = Z (t) \cdot C + ε (t),

(1)

where Z (t) is a row vector and Z_i(t) is the activity of node i at time t. The Gaussian noise ε(t) is spatially and serially uncorrelated, with standard deviation θ = 0.1. The N × N weighted adjacency matrix C = [C_ij] describes the network structure, where C_ij is the weight of the directed connection from node i to node j. These dynamics are generated on various network topologies, as detailed in the following sections. The choice of the weights C_ij (detailed in the following sections) guarantees the stability of the system, which is a sufficient condition for stationarity (Atay & Karabacak, 2006). Since stationary VAR processes have multivariate Gaussian distributions, the information-theoretic measures we use can be directly related to Pearson correlation and Granger causality (Granger, 1969). The simple VAR dynamics are chosen as the primary model for our validation studies instead of nonlinear alternatives because the main goal is not to prove the superiority of nonlinear dependence measures on nonlinear systems, which has been shown elsewhere (Novelli et al., 2019). We rather aim to show that, even on linearly coupled Gaussian variables—perfectly suitable to be studied via cross-correlations—multivariate approaches are better able to infer the macroscopic network properties. In addition, the VAR dynamics are amenable to be investigated using the faster Gaussian estimator for the information-theoretic measures, allowing us to carry out more extensive simulations over a wider range of parameters.

Neural mass model on the CoCoMac connectome.

To provide an extension beyond the linear VAR dynamics, neural activity in various brain regions is modeled (following Shine, Aburn, Breakspear, & Poldrack, 2018; Li et al., 2019) as an oscillating two-dimensional neural mass model derived by mode decomposition from the Fitzhugh-Nagumo single neuron model (FitzHugh, 1961). As previously presented (Li et al., 2019; Shine et al., 2018), the CoCoMac connectome (Kötter, 2004) is used to provide directed coupling between 76 regions, with axonal time delays between these regions based on the length of fiber tracts as estimated by diffusion spectrum imaging (Sanz Leon et al., 2013). Data provided from Li et al. (2019) and Shine et al. (2018) was simulated using the open-source framework The Virtual Brain (Sanz Leon et al., 2013), using code implementing the model freely available at https://github.com/macshine/gain_topology (Shine, 2018).

Langevin equations () specify the neural mass model, via the dynamics of local mean membrane potential (V) and the slow recovery variable (W) at each regional node i:

\begin{matrix} {\dot{V}}_{i} (t) = 20 (W_{i} (t) + 3 V_{i} {(t)}^{2} - V_{i} {(t)}^{3} + γ I_{i}) + ξ_{i} (t), \\ {\dot{W}}_{i} (t) = 20 (- W_{i} (t) - 10 V_{i} (t)) + η_{i} (t) . \end{matrix}

(2)

In the above, ξ_i and η_i are independent standard Wiener noises and I_i is the synaptic current

I_{i} = \sum_{j} C_{j i} S_{j} (t - τ_{j i}),

(3)

with C_ji indicating the connection weight from j to i and incorporating time delays τ_ji from j to i (estimated as described above). The CoCoMac connectome network contains 1,560 directed connections (including 66 self-links), with τ_ji on non-self links having an average of 19.8 ms (standard deviation 8.32 ms). The membrane potentials V_i are converted to normalized firing rates S_i via a sigmoid activation function

S_{i} (t) = \frac{1}{1 + e^{- σ (V_{i} (t) - m)}},

(4)

with parameter m = 1.5 chosen to align the sigmoid with its typical input. The parameters for gain σ = 0.5 (in ) and excitability γ = 0.3 (in ) are selected to simulate activity in the integrated regime of dynamics identified by Shine et al. (2018).

Finally, the time series of membrane voltage V_i(t) (originally obtained with a 0.5-ms temporal resolution via stochastic Heun integration) are subsampled at 15 ms (selected as half the median time for the autocorrelation functions to decay to 1/e).

Network Inference Algorithms

As illustrated in Figure 1, three algorithms are employed to infer network models from the time series using the IDTxl Python package (Wollstadt et al., 2019):

Bivariate mutual information for functional connectivity.

Mutual information (MI) is computed between all pairs of nodes independently, in a bivariate fashion, and only the measurements that pass a strict statistical significance test (described below) are interpreted as undirected links.

MI is a measure of statistical dependence between random variables (Cover & Thomas, 2005), introduced by Shannon in laying the foundations of information theory (Shannon, 1948). Formally, the MI between two continuous random variables X and Y with joint probability density function μ(x, y) and marginal densities μ_X(x) and μ_Y(y) is defined as

I (X; Y) : = \int \int μ (x, y) log \frac{μ (x, y)}{μ_{X} (x) μ_{Y} (y)} d x d y,

(5)

where the integral is taken over the set of pairs (x, y) such that μ(x, y) > 0. The strength of MI lies in its model-free nature, meaning that it doesn’t require any assumptions on the distribution or the variables (e.g., Gaussian). Being able to capture nonlinear relationships, MI is typically presented as a generalized version of the Pearson correlation coefficient. However, for the VAR processes considered here () with stationary multivariate Gaussian distributions, the MI between two variables X and Y under these dynamics is completely determined by the magnitude of their Pearson correlation coefficient ρ (Cover & Thomas, 2005):

I (X; Y) = - ln (1 - ρ^{2}) .

(6)

Crucially, this one-to-one relationship between MI and the absolute value of ρ (for VAR processes) implies that the networks inferred via the bivariate MI algorithm are equivalent to the functional networks obtained via cross-correlation—widely employed in neuroscience. This equivalence persists whenever a Gaussian estimator for MI is used (which models the processes as VAR), even for nonlinear dynamics, as is used in our experiments. Differences may lie in how the raw MI values are transformed into a network structure. Early approaches often used a fixed threshold aimed at obtaining a prescribed link density, while the bivariate MI algorithm used here adopts an adaptive threshold (different for each link) to meet a desired statistical significance level. The statistical significance is computed via null hypothesis testing to reflect the probability of observing a larger MI from the same samples if their temporal relationship were destroyed (the p value is obtained from a chi-square test, as summarized in Lizier, 2014). The critical level for statistical significance is set to α = 0.01/N, where N is the network size. This produces a Bonferroni correction for the inference of parent nodes for each target (i.e., for each target, there is a 0.01 chance under the null hypothesis that at least one spurious parent node is selected, assuming independent sources).

Bivariate transfer entropy for directed functional connectivity.

Transfer entropy (TE) is computed between all pairs of nodes independently, in a bivariate fashion, and only the measurements that pass a strict statistical significance test (described below) are interpreted as links.

TE is a model-free measure of statistical dependence between random variables (Schreiber, 2000); however, different from MI and cross-correlation, it is a directed and not a symmetric measure (i.e., the TE from a source node X to a target node Y is not necessarily the same as the TE from Y to X), and specifically considers information about the dynamic state updates of the target Y. Thus, employing TE has the advantage of generating directed networks and providing a more detailed model of the dynamics of the system under investigation. Formally, the TE from a source stochastic process X to a target process Y is defined as (Schreiber, 2000)

T_{X \to Y} (t) : = I (X_{t - 1}; Y_{t} | Y_{< t}),

(7)

where I(X_t−1; Y_t|Y_<t) is the conditional mutual information (Cover & Thomas, 2005) between the previous sample X_t−1 of the source and the next sample Y_t of the target, in the context of (the vector of) the target’s past values Y_<t. The directed and dynamic nature of TE derives specifically from taking the past of the target into account when measuring the lagged statistical dependence between X and Y.

In practice, in order to estimate the TE from the time series, Y_<t is usually constructed as an embedding vector (Takens, 1981), and a maximum lag must be specified to build a finite embedding of the target’s past (using either uniformly or nonuniformly spaced variables, (Faes, Nollo, & Porta, 2011; Kugiumtzis, 2013; Vlachos & Kugiumtzis, 2010). Here, using the IDTxl Python package (Wollstadt et al., 2019), a nonuniform embedding of the target’s past Y_<t is built via iterative greedy selection of statistically significant elements of Y_<t (as detailed in the next section). While similar multivariate embeddings and lags larger than one time step can be used for the source X in principle, here only a single sample at lag 1 for the source is used for the VAR model in alignment with its dynamics.

Analogously to the MI, the bivariate TE values are transformed into a directed network structure by testing their statistical significance. This is computed (using a theoretical null distribution, as summarized in Lizier, 2014) to reflect the probability of observing a larger TE from the same samples if the source samples were temporally decoupled from the target and its past. A critical level of α = 0.01/N is used as per MI.

Bivariate TE has found wide application in studies of directed functional connectivity (e.g., Honey, Kotter, Breakspear, & Sporns, 2007; Lizier, Heinzle, Horstmann, Haynes, & Prokopenko, 2011; Marinazzo, Wu, Pellicoro, Angelini, & Stramaglia, 2012; Orlandi, Stetter, Soriano, Geisel, & Battaglia, 2014; Stetter, Battaglia, Soriano, & Geisel, 2012; Wibral et al., 2011). Importantly, TE and Granger causality are equivalent for Gaussian variables (Barnett et al., 2009), which applies to the VAR processes considered here (Equation 1), suggesting that a Gaussian estimator for TE be employed (Bossomaier et al., 2016). The networks inferred via the bivariate TE algorithm with this estimator are equivalent (for any dynamics) to those obtained via Granger causality, which is also widely employed in neuroscience.

Multivariate transfer entropy for effective connectivity.

Different from the bivariate approaches above, the multivariate TE approach does not consider pairs of nodes in isolation, but focusses on modeling dynamic updates in each target process in the system by selecting a minimal set of sources that collectively contribute to the computation of the target’s next state. More formally, for each target Y, this method aims at identifying the minimal set of sources X_<t that maximize the collective TE to Y, defined as

T_{X \to Y} (t) : = I (X_{< t}; Y_{t} | Y_{< t}) .

(8)

The multivariate TE network inference algorithm is described in full in Novelli et al. (2019), synthesizing together components of Faes et al. (2011), Lizier and Rubinov (2012), Montalto, Faes, and Marinazzo (2014), Sun et al. (2015), and Vlachos and Kugiumtzis (2010). A greedy approach is used to iteratively select the candidate variables to add to X_<t from a candidate set {X_{i, t−l}} of lagged variables from the past of each source X_i ∈ Z ∖{Y}, up to some maximum lag l ≤ L (L = 1 for the VAR model and L = 4 for the neural mass model used here). Testing the statistical significance of the conditional mutual information for a new candidate source to be included at each step—conditioning on the previously selected sources—provides an adaptive stopping condition for the greedy algorithm. The family-wise error rate for each target is set to α = 0.01 using the max statistic (Novelli et al., 2019), meaning that for each target there is a 0.01 chance that, under the null hypothesis, at least one spurious parent node is selected. The conditioning on previously selected sources serves to prevent the selection of spurious sources that are correlated with true sources because of common driver as well as pathway or chain effects (referred to as holding redundant information only; Lizier & Rubinov, 2012; Stramaglia, Cortes, & Marinazzo, 2014). Such conditioning also enables capturing multivariate or synergistic effects on the target that cannot be detected by examining individual sources in isolation (Lizier & Rubinov, 2012; Stramaglia et al., 2014). Furthermore, in contrast to always conditioning on all potential sources, by conditioning only on previously selected sources the iterative approach defers overly high-dimensional analysis until it is genuinely required, buying statistical sensitivity (“recall”; Lizier & Rubinov, 2012; Novelli et al., 2019). Every node is studied as a target (in parallel on a computing cluster, using IDTxl; Wollstadt et al., 2019) and the results are then combined into a directed network describing the information flows in the system. Similar to the bivariate TE discussed above, a nonuniform embedding of the target’s past Y_<t is built first (Vlachos & Kugiumtzis, 2010), before the second step of selecting sources via the same iterative greedy algorithm (Novelli et al., 2019). While multiple past samples of any given source can be considered (e.g., as has be done in Novelli et al., 2019), only one past value is examined here (L = 1) for the VAR experiments in line with their known structure in Equation 1 and in order to focus on network structure effects only.

Given that TE and Granger causality are equivalent for Gaussian variables (Barnett et al., 2009), using the Gaussian estimator with the multivariate TE algorithm can be viewed as extending Granger causality in the same multivariate/greedy fashion.

Evaluation Metrics

At the microscale (individual links), the network inference performance is evaluated against the known underlying network structure as a binary classification task, using standard statistics based on the number of true positives (TP, i.e., correctly classified existing links), false positives (FP, i.e., absent links falsely classified as existing), true negatives (TN, i.e., correctly classified absent links), and false negatives (FN, i.e., existing links falsely classified as absent). The following standard statistics are employed in the evaluation:

Precision = TP/(TP + FP);
Recall (true positive rate) = TP/(TP + FN);
Specificity (true negative rate) = TN/(TN + FP).
False positive rateFP/(FP + TN).

Intuitively, the precision measures how often an inferred link is actually present in the underlying structure, the recall measures the proportion of true links that are detected, and the specificity measures the proportion of absent links that are correctly not inferred. For a properly controlled family-wise error rate (α), the expected specificity is 1 − α.

At the macroscale, the performance is evaluated in terms of the accuracy in measuring network properties of interest on the inferred network, as compared with their real values when measured on the underlying structural network. These properties include the following;

Characteristic path length. The average shortest distance between all pairs of nodes (Rubinov & Sporns, 2010). A shorter average value is typically interpreted as an indication of the efficiency of the network in propagating information. The characteristic path length is only well defined for connected networks—a problem that is avoided by construction in this study by only generating connected networks. This limitation could be alternatively overcome by replacing the characteristic path length with the analogous global efficiency (Latora & Marchiori, 2001), reported in the Supporting Information for completeness.
Clustering coefficient. The clustering coefficient of a node is the fraction of the node’s neighbours that are also neighbours of each other (Fagiolo, 2007; Watts & Strogatz, 1998). The mean clustering coefficient hence reflects the average prevalence of clustered connectivity around individual nodes.
Small-worldness coefficient. Small-world networks are formally defined as networks that are significantly more clustered than random networks, yet have approximately the same characteristic path length as random networks (Watts & Strogatz, 1998). The small-worldness coefficient was proposed by Humphries and Gurney (2008) to capture this effect in the following single statistic (although improvements on the original measure have been recently suggested; Neal, 2017; Zanin, 2015):
$σ = \frac{\frac{C}{C_{rand}}}{\frac{L}{L_{rand}}} .$
(9)
In Equation 9, C and L respectively denote the average clustering coefficient and the characteristic path length; analogously, C_rand and L_rand denote the average clustering coefficient and the characteristic path length of Erdős-Rényi networks having the same size and number of links as the network under study.
Degree distribution. The probability distribution of the in- and out-degree over the whole network (where the in- and out-degree of a node are defined as the number of its incoming and outgoing connections, respectively).
Betweenness centrality. The fraction of all shortest paths in the network that pass through a given node (Rubinov & Sporns, 2010), excluding paths that start and end on the given node.
Modularity. A measure of the separation of a network into specific groups (or modules), defined as the fraction of the edges that fall within the given groups minus the expected fraction if the edges were distributed at random (Newman & Girvan, 2004).
Rich-club coefficient. The “rich-club” phenomenon refers to the tendency of hubs (nodes with high degree) to form tightly interconnected communities (Colizza, Flammini, Serrano, & Vespignani, 2006). The rich-club coefficient was first proposed by Zhou and Mondragon (2004) to quantify this effect:
$ϕ (k) = \frac{2 E_{> k}}{N_{> k} (N_{> k} - 1)},$
(10)
whereby E_>k denoted the number of edges among the N_>k nodes having degree higher than a given value k.

SMALL-WORLD NETWORKS

Numerical Simulations

The first experiment is aimed at testing the robustness of the three inference algorithms with respect to vast changes in network structure. The Watts-Strogatz model is used to generate a spectrum of topologies, ranging from regular lattices to random networks (similar to Erdős-Rényi networks, although not equivalent; Maier, 2019) through a small-world transition (Watts & Strogatz, 1998). Each simulation starts with a directed ring network of 100 nodes with uniform link weights C_ij = C_ii = 0.15 and fixed in-degree d_in = 4 (i.e., each node is linked to two neighbors on each side, as well as to itself via a self-loop). The source of each link is then rewired with a given probability p ∈ [0, 1], so as to change the overall network topology while keeping the in-degree of each node fixed. Only rewiring attempts that keep the network connected are accepted, in order to allow the measurement of the average shortest path length. The simulations for each p are repeated 10 times on different network realizations and with random initial conditions.

Results

At the microscale, the performance is evaluated in terms of precision, recall, and specificity in the classification of the links (present or absent) in the inferred network compared with the underlying structural network. In the case of bivariate MI, each undirected link in the inferred network is represented as two directed links in opposite directions. For longer time series of 10,000 samples, multivariate TE is the most accurate method, achieving optimal performance according to all metrics on all the network topologies generated by the Watts-Strogatz rewiring model (Figure 2, right column). Bivariate TE also achieves nearly optimal recall and high specificity on all topologies; however, despite the strict statistical significance level, the precision is significantly lower on lattice-like topologies (low rewiring probability) than on random ones (high rewiring probability). The opposite trend is shown by the bivariate MI algorithm, whose precision and recall drastically decrease with increasing rewiring probability. As expected, the recall of all methods decreases when shorter time series of 1,000 samples are provided (Figure 2, left column). However, the recall for multivariate TE is consistent across topologies, while it decreases with higher rewiring probability when bivariate methods are used. This results in the bivariate TE having larger recall for lattice-like topologies, while multivariate TE has larger recall than bivariate for more random topologies (i.e., for a rewiring probability larger than p = 0.2). A further interesting effect is that bivariate TE attains better precision on shorter time series than on longer ones.

Figure 2.

View large Download slide

Performance as a function of the rewiring probability in Watts-Strogatz ring networks (N = 100 nodes). Multivariate TE is consistent across network structure and guarantees high precision regardless of the amount of data available (third row). Bivariate TE (second row) can have better recall than multivariate TE for shorter time series (T = 1,000, left column) but its precision is not consistent (it drops when T = 10,000, right column) and the optimal time series length cannot be determined a priori. Bivariate MI (first row) has lower precision and recall than TE-based methods, for both T = 1,000 and T = 10,000. For each value of the rewiring probability, the results for 10 simulations on different networks are presented (low-opacity markers) in addition to the mean values (solid markers).

At the macroscale, the three algorithms are tested on their ability to accurately measure three fundamental network properties relevant through the small-world transition, using the longer time series of 10,000 samples. Multivariate TE is able to closely approximate the real shortest path length on all the network topologies generated by the Watts-Strogatz rewiring model, while the bivariate MI and TE algorithms produce significant underestimates, particularly on lattice-like topologies (Figure 3). Similarly, multivariate TE is able to closely match the real mean clustering on all the network topologies generated by the Watts-Strogatz rewiring model, while the bivariate MI and TE algorithms consistently overestimate it (Figure 4). The related measure of local efficiency (Latora & Marchiori, 2001) is reported in the Supporting Information.

Figure 3.

View large Download slide

Characteristic path length as a function of the rewiring probability in Watts-Strogatz ring networks (N = 100 nodes and T = 10,000 time samples). Multivariate TE is able to closely approximate the characteristic path length of the real topologies (ground truth). On the other hand, bivariate MI and TE produce significant underestimates due to spurious links creating shortcuts across the network, particularly on lattice-like topologies (low rewiring probability). The results for 10 simulations on different network realizations are presented (low-opacity markers) in addition to the mean values (solid markers).

Figure 4.

View large Download slide

Average clustering coefficient as a function of the rewiring probability in Watts-Strogatz ring networks (N = 100 nodes and T = 10,000 time samples). The multivariate TE algorithm closely matches the average clustering coefficient of the real networks (ground truth), which is instead overestimated by bivariate MI and TE. The results for 10 simulations on different network realizations are presented (low-opacity markers) in addition to the mean values (solid markers).

Given the above results on the characteristic path length and the mean clustering coefficient, it is not surprising that the bivariate MI and TE algorithms significantly overestimate the real small-worldness coefficient, while the multivariate TE method produces accurate estimates on all the network topologies generated by the Watts-Strogatz rewiring model (Figure 5). Equivalent results are found if the alternative measures of “small-world index” (Neal, 2017) or “double-graph normalized index” (Telesford, Joyce, Hayasaka, Burdette, & Laurienti, 2011) are computed instead (not shown).

Figure 5.

View large Download slide

Small-world coefficient as a function of the rewiring probability in Watts-Strogatz ring networks (N = 100 nodes and T = 10,000 time samples). The multivariate TE algorithm produces accurate estimates of the small-world coefficient of the real topologies (ground truth), which is instead strongly overestimated by bivariate MI and TE. The results for 10 simulations on different network realizations are presented (low-opacity markers) in addition to the mean values (solid markers).

Discussion

At the microscale, the results concerning the bivariate TE can be explained in the light of the recent theoretical derivation of TE from network motifs for VAR dynamics (Novelli, Atay, Jost, & Lizier, 2020). For a fixed in-degree, the TE decreases with the rewiring probability, making it harder for candidate links to pass the statistical significance tests when only short time series are available. This explains why the recall for the bivariate TE slightly drops with higher rewiring probability for T = 1000 (Figure 2). We can speculate that a similar mechanism could be responsible for the more drastic drop in the recall for the bivariate MI (via evidence from derivations for covariances from network structure for similar processes; Pernice, Staude, Cardanobile, & Rotter, 2011; Schwarze & Porter, 2020). The fact that bivariate TE is larger for regular lattice structures also explains why its recall is slightly higher than for multivariate TE here: The redundancy between close sources that elevates their bivariate TE is explicitly conditioned out of the multivariate TE for secondary sources. On the other hand, as the rewiring increases, the higher recall for multivariate TE must be due to this method capturing synergistic effects that (more disparate) multiple sources have on the target, which the bivariate method does not.

Comparing the results between shorter and longer time series raises another question: Why is the precision of the bivariate TE worse for longer time series than for shorter ones, especially for lattice-like topologies? More complex motifs involving common parents and multiple walks, which are more prevalent in regular lattice topologies, can result in nonzero TE on spurious links. These indirect effects are typically weak; however, for long enough time series, the low TE values can be distinguished from noise and thus pass the statistical significance tests. The resulting spurious links (false positives) decrease the precision and the specificity as the time series length is increased, with the effect being stronger in regular lattice topologies. In other words, the Bonferroni correction of the statistical significance level (i.e., dividing α by the network size N) does not result in a well-calibrated test for bivariate inference methods—the sources are correlated, and the tests on them are not independent. The differences in the specificity on the plots are subtle because the networks are sparse; however, they manifest in large differences in the precision. Crucially, this effect is not seen for the multivariate TE, which maintains specificity consistent with the requested α for all topologies and time series lengths. Thus, lower recall achieved by multivariate TE on regular lattice networks for short time series (compared with bivariate TE) can be viewed as a compromise to control the specificity in a consistent fashion. A compelling argument in favor of controlling the specificity is provided by Zalesky et al. (2016 p. 407), who conclude that “specificity is at least twice as important as sensitivity [i.e., recall] when estimating key properties of brain networks, including topological measures of network clustering, network efficiency and network modularity.” Unfortunately, there is currently no consistent a priori way (and no reasonable candidate) to determine the optimal time series length for bivariate TE to attain high precision.

Moving to the macroscale results, it is clear that the ability to control the false positives while building connectomes is a crucial prerequisite for the application of complex network measures. Adding only a few spurious links leads to significant underestimation of the average shortest path length—an effect that has previously been reported for lattice-like networks using MI (Bialonski, Horstmann, & Lehnertz, 2010) and extended here to TE and across a range of topologies (Figure 3). Together with the clustering coefficient, the shortest path length is a defining feature of small-world networks. Although evidence of small-world properties of functional networks obtained from fMRI recordings have been provided in several studies (e.g., van den Heuvel, Stam, Boersma, & Hulshoff Pol, 2008), whether or not the brain is a small-world network is still being debated (Hilgetag & Goulas, 2015; Papo, Zanin, Martínez, & Buldú, 2016). Following Papo et al. (2016), the question addressed here is of a pragmatic rather than an ontological nature: Independently of whether the brain is a small-world network, to what extent can neuroscientists using standard system-level neuroimaging techniques interpret the small-world construct in the context of functional brain networks? An indication that the interpretation is problematic was provided by Hlinka, Hartman, and Paluš (2012), who showed that functional connectivity matrices of randomly coupled autoregressive processes show small-world properties. The effect is due to intrinsic properties of correlation rather than just to the finite sample size problem or spatial oversampling. Specifically, correlation has a transitivity property: For any node X with neighbors Y and Z (and respective correlations ρ_XY and ρ_XZ), a lower bound can be derived for the correlation between the neighbors (Langford, Schwertman, & Owens, 2001):

ρ_{Y Z} \geq ρ_{X Y} ρ_{X Z} - \sqrt{1 - ρ_{X Y}^{2}} \sqrt{1 - ρ_{X Z}^{2}} .

(11)

In particular, a strong positive correlation between two pairs of them implies a positive correlation within the third pair: $ρ_{X Y}^{2}$ + $ρ_{X Z}^{2}$ > 1 implies $ρ_{Y Z}^{2}$ > 0 (Langford et al., 2001). The problem was further investigated by Zalesky, Fornito, and Bullmore (2012), who showed that functional connectivity matrices of independent processes also exhibit small-world properties and that—in practice—the correlation between neighbors is much higher than the theoretical lower bound in Equation 11. These considerations on correlation extend to bivariate MI, given the one-to-one relationship between MI and the absolute value of Pearson’s correlation coefficient for the Gaussian variables considered in this study (see the ec5Bivariate mutual information for functional connectivity section in the Methods Section). This transitivity property results in more triangular cliques in functional networks, that is, an inflated clustering coefficient across the whole spectrum of networks in Figure 4. Together with the underestimate of the shortest path length discussed above, the outcome is an overestimate of the small-worldness coefficient (Figure 5). As shown, the limitations of bivariate methods can be overcome by multivariate TE, to a large degree for shorter time series and certainly when sufficiently long time series are available.