Inferring connectivity in networked dynamical systems: Challenges using Granger causality

Determining the interactions and causal relationships between nodes in an unknown networked dynamical system from measurement data alone is a challenging, contemporary task across the physical, biological, and engineering sciences. Statistical methods, such as the increasingly popular Granger causality, are being broadly applied for data-driven discovery of connectivity in fields from economics to neuroscience. A common version of the algorithm is called pairwise-conditional Granger causality, which we systematically test on data generated from a nonlinear model with known causal network structure. Specifically, we simulate networked systems of Kuramoto oscillators and use the Multivariate Granger Causality Toolbox to discover the underlying coupling structure of the system. We compare the inferred results to the original connectivity for a wide range of parameters such as initial conditions, connection strengths, community structures, and natural frequencies. Our results show a significant systematic disparity between the original and inferred network, unless the true structure is extremely sparse or dense. Specifically, the inferred networks have significant discrepancies in the number of edges and the eigenvalues of the connectivity matrix, demonstrating that they typically generate dynamics which are inconsistent with the ground truth. We provide a detailed account of the dynamics for the Erdős-Rényi network model due to its importance in random graph theory and network science. We conclude that Granger causal methods for inferring network structure are highly suspect and should always be checked against a ground truth model. The results also advocate the need to perform such comparisons with any network inference method since the inferred connectivity results appear to have very little to do with the ground truth system.

Figure 1

Inferring network connectivity via Granger causality analysis. (a) Schematics of a coupled dynamical system where the directed network architecture plays an important role. The time series generated by each node is influenced by its connectivity to other nodes. (b) In several applications, the connectivity structure is unknown, but noisy measurements from each node are available. (c) We use Granger causality to infer the original network structure from the noisy data.

Figure 2

Increase in synchrony as connection strength increases. We generate random 12-node Erdős-Rényi networks with a range of connection probabilities. We then solve the Kuramoto model on each network for varying connection strengths. For each network and connection strength pair, we calculate the average order parameter $r\left(t\right)$ [Eq. (6)]. We see that as the connection strength increases, the synchrony also increases. However, for sparse networks, the network remains unsynchronized $\left(r\left(t\right)\approx \frac{1}{\sqrt{n}}\right)$ and for dense networks, the network synchronizes for moderate connection strength. This data was generated in Experiment C1 (see Sec. 5).

Figure 3

Pair of coupled Kuramoto oscillators with distinct natural frequencies. We show the four possible network architectures in panels (a)–(d). We first plot ${\theta }_1$ and ${\theta }_2$, the solution of the differential equations in Eq. (5). We then plot $cos\left({\theta }_1\right)$ and $cos\left({\theta }_2\right)$, the more natural way to view oscillators. In panel (a), the oscillators are uncoupled, so they merely oscillate with their natural frequency. However, in panels (b)–(d), we see cases leading to synchronization. The overall synchronization of the network can be summarized by the parameter $r\left(t\right)$ with full synchronization achieved when $r\left(t\right)=1$ [see Eq. (6)].

Figure 4

Example of synchronicity in structured Kuramoto networks. (a) We have two disjoint subnetworks. The blue oscillators have frequencies with average $-0.2$ while the green oscillators have average frequency 0.5. As we see in the right panel (c), the individual trajectories collapse. The blue nodes synchronize to frequency ${\omega }_A$ and the green nodes synchronize to frequency ${\omega }_B$. In the lower left panel (b), we see that the measure of synchronicity for each community, $r\left(t\right)$ [Eq. (6)], approaches one but at different synchronization times. However, when $r\left(t\right)$ is evaluated on the entire network, we do not achieve total synchronicity because the two communities are not connected.

Figure 5

Overview of steps in our methodology. We will experiment with varying the decisions made in each step—see Sec. 4.

Figure 6

Percentage of inferred edges as the connection strength varies. We try all four possible two-node networks (0, 1, or 2 edges), and we vary the connection strength across the horizontal axis. We also try two amounts of noise. The teal circles are for Experiment A1 (low noise) and the orange closed circles are for Experiment A2 (higher noise). The specific parameter choices for these experiments are in Table 6. The error is higher for low noise. As the connection strength increases, so does the number of inferred edges, a pattern that will continue for larger networks.

Figure 7

Results of Granger causality inference on the two-community network. Panel (a) depicts the true network. The resulting network from Experiment B1 in panel (b) has many extra connections and even connects the two separate communities, but the MVGC Toolbox [29] provides warnings. In Experiment B2, we increase the noise and try again, producing the network in panel (c) without warnings. This network is missing many edges but also connects the two communities. In Experiment B3, we keep the higher level of noise but halve the time step, resulting in the network in panel (d) without warnings. We again have vast overestimation of edges and the community structure is lost. The specific parameter choices for these experiments are in Table 6.

Figure 8

Autocovariance decay for Experiments B2 and B3. In these two experiments, the toolbox does not provide warnings, but there are significant errors (Fig. 7). Here we plot the autocovariance sequence for each experiment to demonstrate that it decays exponentially, as required.

Figure 9

Results from Experiment C1. Here $n=12$ and twenty different Erdős-Rényi networks are generated while varying the percentage of connections. We also vary the connection strength $K$. In panel (a), we plot the true percentage of connections versus the estimated percentage of connections for five values of $K$. If the percentage of connections was correct, our points would be on the dashed diagonal line. However, they may still have the wrong edges even if the correct number are inferred. For the sparse and dense cases, a varying connection strength $K$ is considered in panels (b) and (c). The correct percentage of connections is plotted as a horizontal dashed line for reference. In the inset plots, we see some examples of the inferred networks. These are colored visualizations of the adjacency matrices. White squares denote zeros (no edge) and colored squares denote ones (an edge).

Figure 10

Optimal bands of synchrony. We consider the results of Experiment C1 in terms of average $r$, the synchrony measure in Eq. (6). For each connection probability $p$ across the horizontal axis, we plot a gray line showing the range of average synchrony $r$ attained as we varied the connection strength $K$. We then plot green circles in panel (a) for the values of $r$ for which the error was less than 10%. We also plot orange circles in panel (b) for the values of $r$ for which the error was less than 20%.

Figure 11

Varying number of oscillators. Here we plot the percent connectivity vs estimated percent connectivity and five values of connection strength $K$ for three numbers of oscillators. Left to right, we compare six, twelve, and twenty-four oscillators. The general pattern is consistent, but the average error seems to grow with the number of oscillators.

Figure 12

Varying system parameters. We plot the percent connectivity vs estimated percent connectivity and five values of connection strength $K$ while changing the distribution of initial conditions, the distribution of natural frequencies, and the number of trials. In the first column, we vary the distributions of random initial conditions, and in the second column, we vary the distributions of random natural frequencies. In the third column, we change the number of trials. The general pattern is consistent except when the number of trials is varied; the number of edges inferred grows as the number of trials grows. Results accompanied by a warning are marked with an “$\times$” instead of a circle.

Figure 13

Varying data sampling in time. We plot the percent connectivity vs estimated percent connectivity and five values of connection strength $K$ while changing the data sampling in time. We use data from time 0 to $T$, where $T$ varies down the rows. We use a time step of $\mathrm{\Delta }t$ where $\mathrm{\Delta }t$ varies across the columns. The general pattern is consistent except when the end time $T$ is small. Results accompanied by a warning are marked with an “$\times$” instead of a circle.

Figure 14

Comparing implementations. Here we plot the percent connectivity vs estimated percent connectivity and five values of connection strength $K$. In Experiment D1, we repeat Experiment C1 with the GCCA implementation and observe very little difference.

Figure 15

Comparing implementations. Here we plot the percent connectivity vs estimated percent connectivity and five values of connection strength $K$. In Experiment D2, we compare the MVGC Toolbox to other implementations. We observe that each method infers very few edges when considering only one trial at a time and voting over 50 sets.

Figure 16

Eigenvalue comparison. One measure of accuracy is how well the estimated network $\tilde{A}$ would recreate the same dynamics as the true network $A$. We therefore compare the eigenvalues of $\tilde{A}$ (open colored circles) and $A$ (closed black circles) for the six estimated networks in Fig. 9. We see many cases of eigenvalues being significantly wrong. For example, in the first plot, the inferred network has only eigenvalues of about 0, missing eigenvalues for significant growth.