Inferring causal networks of dynamical systems through transient dynamics and perturbation

Inferring causal relations from time series measurements is an ill-posed mathematical problem, where typically an infinite number of potential solutions can reproduce the given data. We explore in depth a strategy to disambiguate between possible underlying causal networks by perturbing the network, where the forcings are either targeted or applied at random. The resulting transient dynamics provide the critical information necessary to infer causality. Two methods are shown to provide accurate causal reconstructions: Granger causality (GC) with perturbations, and our proposed perturbation cascade inference (PCI). Perturbed GC is capable of inferring smaller networks under low coupling strength regimes. Our proposed PCI method demonstrated consistently strong performance in inferring causal relations for small (2–5 node) and large (10–20 node) networks, with both linear and nonlinear dynamics. Thus, the ability to apply a large and diverse set of perturbations to the network is critical for successfully and accurately determining causal relations and disambiguating between various viable networks.

Figure 1

One simple example to demonstrate the fundamentally ill-posed nature of network inference. Specifically, two Kuramoto networks with coupling strength $K=3$ can easily produce identical time series data, and it is impossible to disambiguate between them from this data alone. We show the best-case predictions of GC where the angles of all oscillators are observed relative to the angle of node 1. When nodes 2 and 3 are perturbed (i.e., their initial angles and frequencies are allowed to be different), GC is able to perfectly resolve both graphs.

Figure 2

An illustration of two approaches to network inference. Given the true network shown on the left, one approach is to passively collect time series observations of each network node (top center) and use time series analysis to infer the network structure shown on the right. If perturbations of the network nodes are allowed, then another approach is to systematically force the network nodes and study their dynamics after each perturbation to infer the network structure.

Figure 3

(a) An example of a 5-node directed mass-spring network with the endpoints attached to two fixed walls. All masses in this simulation are set to 1. (b) The performance of GC on such mass-spring networks varies with network sizes and coupling strengths. When the ratio of coupling strength to mass is close to 1, GC is capable of inferring the network structure of the mass-spring system with 80–85% accuracy for network sizes between 2 and 10 nodes. GC is only able to achieve accuracy above 85% when the network size is 2. Furthermore, if we increase or decrease the ratio of coupling strength to mass away from 1, then the inference accuracy of GC drops below 80% for mass-spring networks of all sizes.

Figure 4

Comparison of Granger inferred networks to the true underlying Kuramoto oscillator network with varied simulation endtimes and coupling strengths. Seemingly large edge prediction accuracies do not resemble the true underlying network structure. Moreover, the eigenvalues of the network connectivity matrix are not well characterized. Therefore, it would not be possible to use an inferred network that is 75–80% accurate to simulate the time dynamics of a system.

Figure 5

GC network prediction of a 5-node Kuramoto oscillator network with a coupling strength of 5 when length of time series data analyzed is stopped at $T=5,10,15,20,25$ with sampling time $\mathrm{\Delta }t=0.1$. Each of the colored lines is a trajectory of one of the five oscillators. Note that the network inferences are best when the length of the simulation passed into GC is close to the transient time length of 50 timepoints.

Figure 6

All experiments are performed on ER networks with 0.5 connection probability and accuracy is averaged over 100 trial simulations where initial conditions are resampled. Time points are sampled at every $\mathrm{\Delta }t=0.1$. (a) Plot of transient time length of Kuramoto network vs. coupling strength and network size. The transient time length of the system depends inversely on the coupling strength, and this dependence stays approximately the same for all networks larger than 5 nodes. Right graph plots the transient time length vs. the inverse coupling strength of Kuramoto networks ranging from 1 to 20 nodes in size. A cluster of linear regions can be seen in the graph with a slope of approximately 20 and this cluster corresponds to all networks 5 nodes or larger. Therefore, for moderately large networks the transient time length is approximately 20 over the coupling strength of the system. (b) GC accuracy over Kuramoto networks with sizes $n=5,10,15$ and varying simulation endtimes and coupling strengths. Black line in heatmap plots corresponds to $T(n,K)=\frac{4.5n}{K}$ (optimal simulation time length is equal to 4.5 times the network size over the coupling strength). Note that for all network sizes, GC achieves a maximum accuracy around $T(n,K)$.

Figure 7

A simple example illustrates the strategy to infer the structure of a network of four coupled nodes using systematic perturbations. We make the important assumption that when a node is activated, it perturbs only its nearest neighbors at the next point in time. In the first perturbation, N1 is activated which causes N2 and N3 to become perturbed which in turn cause N4 to become perturbed. Originally, all edges in the network are given a probability of 0.5. This cascade of perturbations provides new evidence that N1 causes N2 and N3, either N2 or N3 causes N4, and N1 does not cause N4. Therefore, the edges from N1 to N2 and N3 are updated to 1, the edges from N2 and N3 to N4 are updated to 0.66, and the edge from N1 to N4 is updated to 0.05. This process is repeated for the second perturbation of N2 and the third perturbation of N3 with the adjacency matrix of edge probabilities updated each time. The final probability matrix is thresholded at 0.5 which results in the fully accurate network prediction shown in the top right.

Figure 8

Accuracy of PCI, GC, and CCM on perturbed mass-spring and Kuramoto oscillator ER networks (0.5 connection probability) with varying sizes, perturbation forcing magnitudes, and coupling strengths. All results are averaged over 100 simulation trials. For the mass-spring system, PCI achieves $\approx 80\%$ accuracy for all networks of at most 10 nodes and this accuracy increases slightly if the force magnitude by which the nodes are perturbed is increased from 10 to 50. For larger network sizes, PCI requires lower coupling strengths to reach the same predictive accuracy. As an example, for the 20-node Erdős-Rényi networks, PCI achieves this accuracy only when all springs in the mass-spring network have spring constants equal to 0.1. For the Kuramoto oscillator system PCI requires large coupling and forcing strengths $(>20)$ to correctly infer networks of 2 nodes or larger. Therefore, if the coupling and forcing of the system are large, then PCI successfully infers Kuramoto oscillator networks with an accuracy of 80–90%. GC is capable of inferring network structure of a 10-node mass-spring system with 85% accuracy when the coupling strengths are 0.1 or smaller. Similar to PCI, its performance drops below 80% for larger coupling strengths and network sizes and is not significantly affected by the magnitude with which the nodes are perturbed. GC consistently gets below 75% accuracy on Kuramoto oscillator networks with 3 nodes or larger. Note that even with these perturbed systems, CCM is incapable of accurately predicting any structure and for networks of larger than 3 nodes, achieves an accuracy of close to 50%.

Figure 9

The accuracy of PCI on inferring 20-node mass-spring networks, when nodes are perturbed randomly or in order of their outdegree or centrality. We perform this experiment on three different network types: Erdős-Rényi, Erdős-Rényi with cluster, and scale-free. We observe one case (namely, on the clustered ER graphs), where perturbing less than half of the network nodes (8/20) results in an optimistic 85% inference accuracy.

Figure 10

Analysis of accuracy, true positive rate, and false positive rate of PCI over varying number of perturbations and observations of 10 node mass-spring and Kuramoto oscillator networks. No data above the diagonal for each plot is shown because nodes that are perturbed must be observed.

Figure 11

Examples of possible network inference errors when not all nodes in the network are perturbed. (a) When node 1 is perturbed, we do not observe the response of node 2 because it is not observed. Therefore, we mistakingly assume that node 3 becomes activated directly after node 1 so a directed edge from 1 to 3 is inferred. (b) After node 1 is perturbed, we observe that node 3 is perturbed immediately after so we correctly infer a directed edge from 1 to 3. However, node 1 also perturbs node $2^{\prime }$, which activates node $3^{\prime }$, which in turn activates node 4. Therefore, we observe that node 4 is activated directly after node 3 and mistakingly infer a directed edge from 3 to 4.

Figure 12

Analysis of accuracy, true positive rate, and false positive rate of PCI on 5 node and 10 node mass-spring ER networks with connection probability $p=0.5$ and heterogeneous coupling strengths. For every existing edge in each ER network, we choose the coupling strengths uniformly at random from the interval $[0.5-\delta ,0.5+\delta ]$ where $\delta$ is the coupling strength deviation. The left figure shows that PCI maintains an accuracy of 85% regardless of how widely we deviate the coupling strengths of every edge in a 5 node ER network. For 10 node networks, the right figure shows that at $\delta =0$, the accuracy of PCI is approximately 82% and slowly drops to 75% as the coupling strengths become more heterogeneous.

Figure 13

Analysis of PCI for increasing levels of noise and forcing on the mass-spring and Kuramoto systems. Noise is added to the time series trajectories as white Gaussian noise in the frequency domain and y-axes correspond to the variance of the white Gaussian noise.