# Automated Discovery of Interactions and Dynamics for Large Networked
Dynamical Systems

Yan Zhang  Beijing Normal University, Beijing, China  Yu Guo  State Key
Laboratory for Novel Software Technology at Nanjing University, Nanjing, China
Software Institute, Nanjing University, Nanjing, China  Zhang Zhang  Beijing
Normal University, Beijing, China  Mengyuan Chen  Beijing Normal University,
Beijing, China  Shuo Wang  Beijing Normal University, Beijing, China  Jiang
Zhang  Beijing Normal University, Beijing, China

###### Abstract

Understanding the mechanisms of complex systems is very important. According
to specific dynamic rules, a networked dynamical system, understanding a
system as a group of nodes interacting on a given network, is a powerful tool
for modelling complex systems. However, finding such models according to time
series of behaviours is difficult. Conventional methods can work well only on
small networks and for some types of dynamics. This paper proposes a unified
framework for an automated interaction network and dynamics discovery (AIDD)
on various network structures and dynamics, based on a stochastic gradient
descent algorithm. The experiments show that AIDD can be applied to large
systems with thousands of nodes and is robust against noise and missing
information. We further propose a new method to test data-driven models based
on control experiments. The results show that AIDD is able to learn the real
network dynamics correctly.

###### keywords:

network inference, dynamics learning, control, networked dynamical systems,
causal discovery

## Introduction

Living cells, brains, human society, stock markets, global climate systems,
and so forth are complex systems composed of many nonlinear interactive units
[1, 2, 3, 4, 5]. By decomposing a complex system into a static network with
dynamics on nodes, networked dynamical system models are powerful tools to
describe complex systems, playing a paramount role in understanding their
collective behaviours and controlling their functions [2, 3, 6]. However,
building such models requires professional knowledge and modelling experience,
which hinders the wide application of these methods. The reconstruction of
such networked dynamical systems in a data-driven way remains a fundamental
problem, i.e., to retrieve the interaction network structure and the node
dynamics from time-series data of complex system behaviours without any
subjective biases [7, 8].

Although many classical approaches to time series forecasting have been
proposed [9, 10], prediction of the behaviours of complex systems with highly
nonlinear and long-range correlations, especially those based on a complex
network structure, had not been resolved until the recent introduction of
graph (neural) network (GNN) models [11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22]. GNNs are designed particularly for networked dynamical systems. By
learning complex functions of information aggregation and propagation on a
given network, GNNs can simulate any complex dynamics on such networks.
However, a complete graph is always required for most GNN models, which
hinders their wider applications [19, 23, 24, 25, 8, 26].

How to reveal network structure from time series data is of great importance
because the revealed interaction network can not only help us to understand
the behaviours of complex systems but also can improve algorithms’
explainability and transparency in terms of causality [27, 28, 29, 30, 5]. The
interdependence relations or causal structure can be obtained by directly
calculating some statistical measures [31, 32, 33], perturbing the system [34,
35], optimising a score function [36, 37], or expanding the complex
interaction dynamics on a set of basal functions [7, 38, 39], and other
methods [40, 41, 42].

Among these algorithms, the algorithm for revealing network interactions
(ARNI) [43] is one of the most prominent methods. It can not only infer a
network with high accuracy but can also be adopted for various nonlinear
dynamics. However, one disadvantage is that the performance of the model
strongly depends on the choice of the basal functions. If the prior biases on
basal functions are missing, this approach becomes very time-consuming,
limiting its application to larger systems.

Very few studies have been proposed to perform both network inference and time
series forecasting tasks together, although some network inference algorithms
are capable of forecasting. The implicit and time-variant network structures
can also be obtained from deep learning models for forecasting based on an
attention mechanism [25, 44, 45, 8, 46, 47]. The first framework to derive an
explicit network is NRI (neural relation inference) [19], in which an
encoder–decoder framework is used. However, the complicated encoding process
to infer the connections from time series data has limited scalability and
accuracy on larger networks [24]. Further studies promote NRI on several
aspects, such as the consideration of constraints [48, 23], and inferring
hidden nodes [49], whereas the problem of scalability and accuracy remains.

However, despite the large number of experiments on evaluating and comparing
various data-driven models that have been conducted based on prediction tasks,
much room for improvement remains. According to the three-layer causal
hierarchy, intervention and counterfactual experiments rather than predictions
are the gold standards for testing a trained data-driven model as a complete
substitution of the original system [27, 50]. That is, downstream tasks such
as the control experiments of both the learned model and the original systems
should be tested and compared [51]. However, we found only a few studies have
performed this kind of task [51, 6].

Figure 1: The workflow and evaluation of Automated Interaction and Dynamics
Discovery(AIDD). Workflows showing how our proposed AIDD framework models a
complex system, and how can be evaluated on tasks of time series forecasting,
interaction network inference, and control experiments. The framework of the
AIDD is also shown in the inset box. A column of the adjacency matrix for the
candidate network is sampled by the network generator. It can be regarded as a
mask vector to filter out the supposed unrelated nodes. Then, the time series
information for related nodes is input into the dynamics learning module,
which then outputs a prediction for the new states of all nodes. After that,
the prediction is compared against the data. The loss function can be
calculated, and the gradient information can be back-propagated directly.
After optimisation, a learned networked dynamical system represented by neural
networks can be obtained.

As shown in Fig. 1, this paper proposes a unified framework for automated
interactions and dynamics discovery (AIDD). This is a universal framework for
learning both the interaction structure and dynamics of a complex system from
time series data. The design of a lightweight network generator and a
universal dynamics learning component based on Markov dynamics makes AIDD not
only applicable to various networks and dynamics, but also enables it to
reconstruct very large networks with high accuracy and robustness. The entire
framework is differentiable so that it can be optimised directly by automatic
differentiation and machine learning techniques [52]. Beyond tasks of network
inference and time series forecasting, we propose a new method to test a
learned data-driven model based on control experiments. Finally, we test the
validity of our framework on real gene regulatory networks under noisy,
incomplete, and interrupted data, which is close to realistic situations. The
results demonstrate that a high performance can be obtained.

## Results

### Problem Formulation

Suppose the complex system to be considered evolves under discrete time steps.
Thus, the dynamics to be reconstructed can be described by a mapping.

$\mathbf{X}^{t+1}=f(\mathbf{X}^{t},A)+\mathbf{\zeta}^{t},$ (1)

where
$\mathbf{X}^{t}=[X_{1}^{t},X_{2}^{t},\cdot\cdot\cdot,X_{N}^{t}]\in\mathbb{R}^{N\times
D}$ is the state of the system at time $t$, $N$ is the number of nodes, $D$ is
the dimension of the state space of each single node, $A$ is the adjacency
matrix of the interaction network to be inferred, and
$\mathbf{\zeta^{t}}\in\mathbb{R}^{N\times D}$ is the noise imposed on nodes.
However, Equation 1 can only describe the dynamical processes with explicit
mathematical forms and cannot be applied to those defined by rule tables or
transitional probabilities, such as cellular automata, Boolean dynamics, or
Markov chains. Therefore, instead of Equation 1, we use a more general form, a
Markov chain $\\{\mathbf{X}^{t}\\}$, to describe the dynamics.

$f(\mathbf{X}^{t+1}|\mathbf{X}^{t},A)\equiv
P(\mathbf{X}^{t+1}|\mathbf{X}^{t},A),$ (2)

where $f$ is the dynamics to be discovered, $\mathbf{X}^{t}$ is the
abbreviation for the event that the random variable $\mathbf{x}^{t}$ takes
value $\mathbf{x}^{t}\in S^{N}$, where S is the state space of each single
node and can be either a finite set of discrete values or an infinite set with
continuous values. $P$ is the conditional probability. Equation 2 is
compatible with Equation 1 but more general [53]. It can even be extended to
non-Markov random processes with finite histories by adding more hidden
auxiliary variable variables[54].

However, it is difficult to infer the probabilities in Equation 2,
particularly when $N$ is large. Fortunately, the interactions of complex
systems are always localised, which means that
$P(\mathbf{X^{t+1}}|\mathbf{X^{t}},A)$ can be factorised into local
transitional probabilities [55]:

$P(\mathbf{X^{t+1}}|\mathbf{X^{t}},A)=\prod\limits_{i=1}^{N}P(X^{t+1}_{i}|\mathbf{X^{t}}\odot
A_{\cdot i}),$ (3)

where $\odot$ is the element-wise product, and $A_{\cdot i}$ represents the
$i^{th}$ column of matrix $A$, and $\mathbf{X^{t}}\odot A_{\cdot i}$ is a
vector representing the state combination of all neighbour nodes of $i$. Then

$f_{i}(X^{t+1}_{i}|\mathbf{X^{t}}\odot A_{\cdot i})\equiv
P(X^{t+1}_{i}|\mathbf{X^{t}}\odot A_{\cdot i})$ (4)

represents the local dynamics of node $i$, which is also called a causal
mechanism in the literature [55].

Therefore, our task becomes the reconstruction of the network $A$ and learning
the local dynamics $f_{i}$ according to the observed time series
$\mathbf{x}=(\mathbf{x}^{1},\mathbf{x}^{2},\cdot\cdot\cdot,\mathbf{x}^{T})$
with $T$ time steps on multiple samples.

### Model

We build a neural network framework consisting of two modules to solve the
reconstruction problem, as shown in the inset panel of Fig. 1. The first
module is a network generator that can generate a candidate network adjacency
matrix $\hat{A}(\Theta)$ with the parameters $\Theta$. The second module then
attempts to simulate the system dynamics $f$ by using a set of neural networks
$\hat{f}_{i}(\hat{X}_{i}^{t+1}|\mathbf{x}^{t}\odot\hat{A}_{\cdot
i}(\theta_{\cdot i}),\phi_{i})$ for any node $i$, which are parameterized by
$\Phi=(\phi_{1},\cdot\cdot\cdot\phi_{N})$ to predict the future state
$\hat{\mathbf{X}}^{t+1},\hat{\mathbf{X}}^{t+2},\cdot\cdot\cdot$ according to
the candidate matrix $\hat{A}(\Theta)$ and the observed state of the previous
time step $\mathbf{x}^{t}$.

Instead of using the complicated graph network architecture to generate the
candidate network as described in [19], we directly sample each element in the
adjacency matrix $\hat{A}_{ij}\sim Bernoulli(\theta_{ij})$, where
$\theta_{ij}\in[0,1]$ represents the probability that the entry of the
$i^{th}$ row and the $j^{th}$ column in $\hat{A}$ takes the value $1$. To make
the sampling process differentiable, we use the Gumbel-softmax technique [56,
19] to generate the adjacency matrix [24].

$\hat{A}_{ij}=\sigma((\log(\theta_{ij})+\xi_{ij})/\tau),$ (5)

where $\xi_{ij}\sim Gumbel(0,1)$ is a random number following the standard
Gumbel distribution, $\sigma$ is the softmax function, and $\tau$ is the
parameter of temperature to adjust the softness of the sampling process. The
random numbers generated by Equation 5 have a similar distribution as
$Bernoulli(\theta_{ij})$, especially when $\tau$ is large. The simulated
sampling process is differentiable such that the gradients can be passed by.
When $\tau\rightarrow\infty$, $\hat{A}_{ij}$ exactly equals 1 with probability
$\theta_{ij}$, or 0 with probability $1-\theta_{ij}$.

Compared to other network generation mechanisms based on the Hadmard product
of two $V$ dimensional node feature vectors [19, 8, 44], where $V\ll N$, our
method has higher accuracy in inferring links because more parameters
($N\times N$ v.s. $N\times V$) are used.

However, optimising $N\times N$ parameters does not reduce the implementation
performance because the matrix elements are independent of each other, which
means that $\hat{A}$ can be generated column-by-column separately, where each
column represents the possible neighbourhood of a single node. In this way,
our framework has a large improvement in flexibility and computational
efficiency compared to the encoder-decoder frameworks such as [19] and can be
applied to very large networks. However, the limitations are that the networks
should be static, and the correlations between elements of $\hat{A}$ are
ignored. Further, the introduction of noise $\xi_{ij}$ can push the network
generator to jump out of local minimums during optimisation.

According to Equation 3, the dynamics learning module can also be decomposed
into local modules node by node. Each local dynamics learner can be modelled
by a feedforward or a recurrent neural network $\hat{f}_{i}(\cdot|\phi_{i})$.
Both the network structure and the parameters $\phi_{i}$ can be shared by
different nodes or not. If they share the dynamics learning module, they form
a graph network [11, 19], as shown in Fig. 1.

### Objective function

Finally, the network inference and dynamics learning problem can be solved by
minimising the following objective function.

$\begin{split}L(\Theta,\Phi)&=\mathbb{E}_{\hat{A}\sim
B(\Theta)}\left(-\sum\limits_{t=1}^{T}\log
P(\hat{\mathbf{X}}^{t+1}=\mathbf{x}^{t+1}|\mathbf{x}^{t},\hat{A},\Phi)\right)+\lambda\sum_{ij}\hat{A}_{ij}\\\
&\approx-\frac{1}{K}\sum_{k=1}^{K}\sum_{t=1}^{T}\sum_{i=1}^{N}L_{i}(\hat{X}_{i}^{t+1}=x_{i}^{t+1}|\mathbf{x}^{t},\hat{A}_{\cdot
i}^{k};\phi_{i})+\lambda\sum_{ij}\hat{A}_{ij},\end{split}$ (6)

where

$L_{i}(\hat{X}_{i}^{t+1}=x_{i}^{t+1}|\mathbf{x}^{t},\hat{A}_{\cdot
i}^{k};\phi_{i})=\log\hat{f}_{i}(\hat{X}_{i}^{t+1}=x_{i}^{t+1}|\mathbf{x}^{t}\odot\hat{A}^{k}_{\cdot
i};\phi_{i})$ (7)

is the local log-likelihood, and $K$ is the number of samples for matrix
$\hat{A}$ under a given $\Theta$, and $x_{i}^{t}$ is the observational vector
of states of node $i$ at time $t$. Thus, the objective function contains two
terms, the former being the log-likelihood, which can be decomposed into local
terms. The latter is the structural loss to conform the network to be sparse
while avoiding over-fitting [48]. The parameter $\lambda$ can adjust the
relative importance of the structural loss. When we perform multiple time
predictions, we can use $\hat{\mathbf{X}}^{t}$ predicted by the previous time
series instead of $\mathbf{x}^{t}$ in Equation 7 to calculate the loss.

If the state space $S$ of each node is real, then the local log-likelihood
Equation 7 can be taken as a mean-absolute error (MAE) form,

$L_{i}(\hat{X}_{i}^{t+1}=x_{i}^{t+1}|\mathbf{x}^{t},\hat{A}_{\cdot
i}^{k})=\left|x_{i}^{t+1}-\hat{f}_{i}(\mathbf{x}^{t}\odot\hat{A}^{k}_{\cdot
i};\phi_{i})\right|,$ (8)

by assuming that $\hat{f}_{i}(x^{t+1}_{i}|\mathbf{x}^{t}\odot A_{\cdot i})$ is
an independent Laplacian distribution for any data point
$\frac{1}{(2\nu)}\exp{\left[-\frac{1}{\nu}\cdot\left|x_{i}^{t+1}-\mu_{i}\left(\mathbf{x}^{t}\odot\hat{A}_{\cdot
i}|\phi_{i}\right)\right|\right]}$, where, the mean value of the Laplacian
distribution $\mu_{i}$ is modelled by a neural network
$\hat{f}_{i}(\mathbf{x}^{t}\odot\hat{A}_{\cdot i};\phi_{i})$, which can be
regarded as a new form of the dynamics learner [57]. We let $\nu=1$ in this
paper.

Then, the network dynamics to fit the observational data can be obtained by
optimising the objective functions Equation 8 node by node [58]. We use the
stochastic gradient descent algorithm to optimise. More details about training
and testing can be found in the Methods section and the Supplementary Section
1.

### Performances and Comparisons

To test our model and compare it with others, we generate time series data by
a set of network dynamics on a variety of graphs as the ground truths. All
continuous, discrete, and binary dynamics are included. Spring (spring
dynamics) [19], SIR (an inter-city meta-population SIR epidemic model) [59],
and Michaelis–Menten kinetics [60, 61] all are examples of continuous
dynamics. The coupled map network (CMN) [62, 63] and the voter model [64] are
representatives of discrete and binary dynamics, respectively. The details of
the dynamical models can be found in the Supplementary Section 6. The graphs
to be inferred are either generated by models (ER [65] for Erdos Renyi, WS [4]
for Watts-Strogatz, and BA [66] for Barabasi-Albert) or from empirical data
including a gene network (for S. cerevisiae, Gene) [67], an inter-city traffic
network (City) of China, three social networks (email, Dorm, and Blog), and a
road network(Road) [68]. The details of the empirical networks can be found in
the Supplementary Section 7. For each model and network, we ran the simulation
for various lengths of time. All the reported results are on the testing
dataset.

We compare our model to a series of baseline methods for both network
inference and single step forecasting tasks. ARNI [43] and NRI [19] are both
state-of-the-art models for network inference and time series forecasting. The
former is based on the block-orthogonal regression method, and the latter is
based on deep learning and graph networks. Two other frequently used
statistical metrics, partial correlation and mutual information, are also
compared on the network inference task. In addition, the classic time series
forecasting model, long short-term memory (LSTM) [54], is compared on the
prediction task. The details of these comparison models can be found in the
Supplementary Section 5. In all the experiments, we used a graph network with
128 hidden units as the dynamics learning module (see the method section for
details). The neural network is shared by all nodes. Other parameters are
shown in the footnote of Table 1 and Supplementary Section 9.

Table 1: Comparisons of performance on network inference and dynamics
prediction tasks between AIDD and other selected methods (columns) on
different dynamics and networks (rows)

Type | Model | Network | ARNI | MI | PC | NRI | LSTM | OURS
---|---|---|---|---|---|---|---|---
AUC | MSE | AUC | AUC | MSE/ACC | AUC | MSE/ACC | MSE/ACC | AUC
Con. | Spring | ER-10 | 0.5853 | 1.33E-03 | 0.7500 | 0.8250 | 2.60E-08 | 0.9998* | 2.98E-04 | 2.70E-04* | 1.0
WS-10 | 0.5125 | 1.58E-03 | 0.6875 | 0.7875 | 8.40E-08 | 0.9997* | 3.35E-04 | 3.31E-04* | 1.0
BA-10 | 0.5169 | 1.10E-03 | 0.6422 | 0.6571 | 7.00E-10 | 0.9999* | 2.14E-04 | 2.90E-05* | 1.0
ER-2000 | - | - | 0.4997 | - | - | - | 2.25E-03 | 1.18E-05 | 0.9886
WS-2000 | - | - | 0.5002 | - | - | - | 5.89E-03 | 8.51E-06 | 0.9933
BA-2000 | - | - | 0.5010 | - | - | - | 4.54E-03 | 2.09E-03 | 0.9523
SIR | SIR-371(D) | 0.5424* | 8.25E-03 | 0.5027 | 0.5119 | - | - | 2.28E-03* | 2.98E-05 | 0.9156
Menten | Gene-100(D) | 1.0 | 9.71E-03 | 0.5416 | 0.6574 | - | - | 2.29E-03* | 4.37E-05 | 0.9960*
Dis. | CMN | ER-10 | 1.0 | 2.33E-09 | 0.5745 | 0.7804 | 1.40E-05 | 0.8850 | 2.60E-04 | 5.60E-06* | 1.0
WS-10 | 1.0 | 2.35E-09 | 0.6875 | 0.8375 | 9.40E-06 | 0.9331 | 2.40E-04 | 2.80E-06* | 1.0
BA-10 | 1.0 | 2.40E-09 | 0.4390 | 0.7439 | 1.30E-05 | 0.6753 | 9.21E-05 | 6.90E-06* | 1.0
ER-200 | 0.8441* | 4.17E-02 | 0.5774 | 0.7648 | - | - | 5.91E-05 | 2.04E-06 | 0.9987
WS-200 | 1.0 | 2.36E-09 | 0.6969 | 0.7506 | - | - | 1.63E-04 | 1.95E-06 | 0.9987*
BA-200 | 0.8840* | 2.45E-02 | 0.5533 | 0.7493 | - | - | 1.46E-04 | 2.57E-06 | 0.9874
WS-1000 | - | - | 0.5670 | - | - | - | 3.54E-05 | 2.92E-06 | 0.9795
BA-1000 | - | - | 0.5290 | - | - | - | 3.46E-05 | 5.48E-05 | 0.9105
Bin. | Voter | ER-10 | - | - | - | 0.4552 | 0.8447* | 0.5000 | 0.5242 | 0.9647 | 1.0
WS-10 | - | - | - | 0.5250 | 0.9062* | 0.5037 | 0.6007 | 0.9463 | 1.0
BA-10 | - | - | - | 0.4607 | 0.9588* | 0.4999 | 0.6917 | 0.9866 | 1.0
WS-1000 | - | - | 0.5470* | - | - | - | 0.5317 | 0.6650 | 0.9996
BA-1000 | - | - | 0.5030 | - | - | - | 0.5208* | 0.6758 | 0.9942
EMAIL-1133 | - | - | 0.4999 | - | - | - | 0.5333* | 0.7212 | 0.9576
ROAD-1174 | - | - | 0.5004 | - | - | - | 0.5455* | 0.8942 | 0.9996
DORM-217(D) | - | - | 0.5219 | - | - | - | 0.5735* | 0.6951 | 0.9901
BLOG-1224(D) | - | - | 0.4995 | - | - | - | 0.5295* | 0.6793 | 0.8603

  * •

In the network column, we use network - size format. The networks marked with
“D” means that they are directed graphs. All networks generated by models
share the same edge density value, which is 1% for large networks (size > 10),
and it is 20% and 3% for small networks with sizes smaller than 10, and ER
networks with size = 200, respectively, to avoid isolated nodes. All the
results are the averages of five repeated experiments. The same data volume is
shared for different methods in one row. The items marked by “-” indicate that
valid results of the model cannot be obtained due to the limitations of the
specific method on dynamics, memory, or time consumption. The best results
among all the compared algorithms in the same row are boldfaced, and the
second-best results are marked “*”. More parameter settings are shown in the
Supplementary Section 9.

Figure 2: Multi-step prediction results of AIDD on CMN dynamics data in a
10-node ER network. In (a), we show the time series data of multi-step
predictions and the ground truths for two selected nodes. In (b), we show how
the mean square error (MSE) increases with time for CMN dynamics. The
parameters are the same as in Table 1.

As shown in Table 1, our model outperforms all the methods on large networks
for both tasks. Compared with the ARNI model, AIDD does not rely on the choice
of basis functions; as a result, it can be applied to very diverse dynamics.
By using neural networks as a universal estimator for dynamics, AIDD avoids
the problem of combinatorial explosion of basis functions. This enables AIDD
to have competitive performances on space and time complexity ($O(N^{2})$ v.s.
$O(N^{3})$ for ARNI. In ARNI model, the time complexity of finding the
interactions of a single node is $O(N^{2})$, so the time complexity of
revealing the entire network interactions is $O(N^{3})$). In order to compare
performance on time series forecasting, we slightly modified the original
algorithm of ARNI such that it can also output one-step prediction [43].
Compared to the NRI framework, our model has a much lighter network generation
architecture. NRI cannot output any result under the same limitations of time
and space on networks with sizes larger than 30 owing to the computational
complexity [19].

The model can also output multi-step prediction results by feeding the result
of the one-step prediction output back to the model. Figure 2 shows the
results for the selected dynamics.

Figure 3: The performance of AIDD under different factors in network inference
and dynamics learning. (a) shows how the number of nodes and the volume of
data (the number of samples$\times$ the number of time steps, which was fixed
to 100) jointly influence the network inference accuracy on WS networks under
CMN dynamics. With the exception of 100 nodes with a 4% edge density, all
nodes shared the same edge density value, which was 1%. (b) shows how
performance decreases with edge density. For the experiments in (b), we set
the number of nodes to 100, and the sparse matrix parameter $\lambda$ was set
to 0.

In general, AIDD works very well on large sparse networks, and the performance
on both tasks decreases as the edge density increases, as shown in Fig. 3(b).

We can improve the accuracy by increasing the amount of data. Figure 3(a)
shows how the area under curve (AUC) depends on both network size and data
volume to be fed into the model systematically, and similar contours can be
obtained on MSE errors (see Supplementary Section 2). There is a trade-off
between network size and data volume under a given accuracy, as shown in Fig.
3(a). It is interesting to observe that data volume is sensitive to network
size only when the number of nodes is between 300 and 500, and beyond that, a
minimum amount of data volume is sufficient to obtain an acceptable accuracy
(e.g., $AUC=0.7$), and this almost does not depend on how large the network
is. We suppose that the given total number of epochs is sufficient for
training the model only for networks smaller than 300, which is the reason why
a sudden increase in data volume is observed. When the size is larger than
500, the model can converge quickly when sufficient data volume is provided;
therefore, the curves become insensitive again.

### Robustness against noise and hidden nodes

A good data-driven model must be robust against noise and unobservable nodes
such that it can be applied to the real world. To show the robustness against
noise of AIDD, we plot changes in AUC with the magnitude of noise on
Michaelis–Menten kinetics [60, 61], which can describe the dynamics of Gene
regulatory networks, as shown in Fig. 4. Our model can recover the network
structure with 0.85 AUC when the mean magnitude of noise is 0.3.

In real applications, we can only obtain partial information of the entire
system owing to the limitations of observation. Thus, a certain proportion of
nodes are unobservable or hidden. This requires the inference algorithm to be
robust to hidden nodes. Thus, we test the AIDD on an incomplete network. To
generate the incomplete network data as the ground truth, we randomly select a
certain percentage of nodes as the hidden nodes (Fig. 4(a)), and the time
series data of these nodes are removed. AUC decreases and MSE increases as the
fraction of the number of unobserved nodes increases on both spring and voice
dynamics, as shown in Fig. 4(c); however, the sensitivity depends on various
types of dynamics. It is found that when the proportion of missing nodes
reaches 50$\%$, the inference accuracy is still above 95$\%$, which proves
that our model can achieve superior results in the absence of normally
sufficient amounts of data.

Furthermore, we test the ability of AIDD to reveal unknown network structures
of unobservable nodes on CMN and Voter dynamics, with only the number of
hidden nodes available. We completed this task by performing the same
interaction inference task, setting the states for unknown nodes to random
values. Figure 4(d) shows the AUCs of the link structures of unknown networks
on Voter and CMN dynamics. The results reveal that the network inference
accuracy is robust for missing nodes. The algorithm can recover the
interactions even for unobservable nodes with over 80% accuracy. The details
of the algorithm can be found in [69] and Supplementary Section 3.

Figure 4: The robustness evaluation of AIDD against noise and missing nodes.
(a) shows a schematic ground truth network with missing information on the
unobserved nodes (grey nodes). (b) shows the influence of proportion of
unobserved nodes on the accuracy of interaction inference on the partial
network with observed nodes measured by AUC, and the accuracy of dynamic
predictions (inset) measured by the MSE of the observable nodes on Spring,
CMN, and the AUC of the observable nodes on Voter dynamics. All the
experiments were conducted on an ER network with 100 nodes, and all networks
generated by models share the same edge density value, which is 4%. (c) shows
the dependence of AUC and MSE on the mean of noise added on each node for the
Michaelis–Menten kinetics (Gene dynamics) on the yeast S. cerevisiae gene
network with 100 nodes. (d) shows the ability of AIDD to infer interactions on
the entire network (the light colour bars) and the unobserved partial networks
(the dark colour bars). All the experiments are conducted on CMN and Voter
dynamics with ER, WS, and BA networks, and all networks contain 100 nodes with
10% unobservable nodes selected randomly, and all networks generated by models
share the same edge density value, which is 4%. The parameters in these
experiments can be referred to in the Supplementary Section 9.

### Control

To further verify that AIDD has learned the ground truth dynamics and that the
trained model can replace the original system, we design control experiments.
The reasons why we choose the control problem as the test bed for our model
include (1) the control problem in a complex network is very important, and it
has relevance to many engineering fields [6]; (2) the control problem is more
difficult than predictions. This is true because to control a system means to
intervene in it. As a result, we have stood at least on the second level of
the causal hierarchy [50].

Here, our control problem is to find a control law based on the learned
network dynamics such that if the rule is applied to the ground truth model,
we can synchronise the whole system by regulating a few nodes. To do this, we
divide the nodes in the network into two groups: the driver nodes, who can be
manipulated by the controller directly, and the target nodes, whose states
cannot be directly controlled, but can be influenced indirectly. The
controller is designed by optimising a set of parameters such that the control
objective, that is, the synchronisation of the whole system, is achieved.

The control experiments consisted of two stages as shown in the bottom part of
Fig. 1. In the first stage, we find the optimised controller’s parameters on
the learned network dynamics to achieve the designed objective. In the second
stage, we apply the optimised controller to the ground truth model. In order
to verify the learned model, we also optimize another controller on the ground
truth model. After that, we compare the deviation curves of the two
controllers from the control target on the same ground truth model. If the
curves are similar, we conclude that the learned model can substitute for the
real model. The detailed methods can be referred to in the Supplementary
Section 4.

Two control experiments were designed. The first problem is to synchronise the
movement directions of all masses with spring dynamics on a small-size BA
network with 10 nodes (see Fig. 5(a)). Three nodes with the largest degrees
are selected as the driver nodes. The controller implemented by a neural
network adjusts the forces imposed on the drivers according to the current
state of the system at each time step. The experimental results are shown in
Fig. 5. The two MSE curves in Fig. 5(b) depict the degree of control achieved
by the trained model and the ground truth model, respectively. They overlap to
show that the learned dynamics can be a good substitute for the real system.
Both curves approached 0 within 13 timesteps, indicating that the controllers
achieved the goal within the given time.

The second control experiment involves querying all the oscillators in a CMN
model on a WS model network with 10 nodes (see Fig. 5(d)) to take the same
value of $0.6$, which is the mean value an oscillator can take with a given
range. The controller is a parameterized neural network, which maps the
current state of the system into the control signals (see Supplementary
Section 4). The signals are the forces imposed on the two drivers (with the
largest degrees) directly. From Fig. 6 (e) and (f), the controls are not well
achieved for all nodes because the MSE curves do not converge to zeros.
However, the two MSE curves overlap very well, indicating that the surrogate
behaves identically to the ground truth model.

Figure 5: The control experiments on learned models. (a) shows the spring
network that we studied for the first experiment. Three large nodes are driver
nodes, and the others are target nodes. The control objective is to request
all masses to have the same movement direction. (b) shows the final movement
states of all the target nodes under the controls. (c) shows two MSE curves
for evaluating goal achievement versus time steps of the controls. One
represents the results of learned model, and the other is the ground truth.
(d) is a coupled mapping network that we studied in the second experiment. Two
large nodes were selected as driver nodes. The control objective is to ask all
oscillators to have the same value of 0.6, which is the mean of the value
range for all nodes. (e) shows the oscillations of all target nodes during
control. (f) shows two MSE curves for evaluating goal achievement versus time
steps of the controls. One is for the trained model, and the other for the
ground truth. The parameters in these experiments are given in the
Supplementary Section 9.

### Gene regulatory network inference

To verify that our algorithm can be applied to actual scenarios and not only
on toy models, we attempt to infer the real subnetwork structure from the
known transcriptional network of yeast S. cerevisiae according to the time
series data of mRNA concentrations generated by GeneNetWeaver (GNW) [67], a
famous simulator for gene dynamics.

The networks used by GNWs are extracted from known biological interaction
networks (Escherichia coli, Saccharomyces cerevisiae, etc.). On these
networks, GNW uses a set of dynamical equations to simulate the transcription
and translation processes, and it has considered many factors close to real
situations (see Supplementary Section 6 for details). Therefore, GNW is a
famous platform for benchmarking and performance assessment of network
inference methods.

In the experiment, we used yeast S. cerevisiae gene network with 100 nodes as
the benchmark gene network, and we used the default parameters of DREAM4_In-
Silico in GeneNetWeaver software to generate data. For the dynamics learner,
we use different neural networks for each node because of the heterogeneity of
node dynamics and the existence of latent variables, noise, and perturbations
[67]. We compare our method with partial correlation, Bayesian network
inference, and mutual information algorithms. Our method outperforms others on
network inference (Fig. 6(a)) on the AUC(0.82). It can also predict the
dynamics with a relatively high accuracy (the average absolute error (MAE) is
0.038, see Fig. 6. This indicates that our method can perform well realistic
gene regulatory dynamics.

Figure 6: Performances of AIDD and other compared methods on network inference
for the gene regulatory network of yeast S. cerevisiae with 100 nodes. (a) ROC
curves of different network inference methods. The comparison methods include
Bayesian network (BN), partial correlation (PC), mutual information (MI), and
AIDD. The AUC for different methods is marked in the legend. (b) shows the
comparison between the observed time series of the expression data (real) and
the predicted data on selected genes. In this plot, the solid lines represent
the predictions and the dotted lines represent the observed data.

## Discussion

In this paper, we propose a unified framework of automatic interaction and
dynamics discovery, called AIDD, for large network dynamical systems. We also
propose a new standard based on control tasks to evaluate whether the true
network dynamics can be learned.

The main highlights of AIDD include scalability, universality, and robustness.
The high scalability is reflected by the fact that the model can be applied to
large networks with more than thousands of nodes with more than 90% accuracy
because the training procedure can be taken node by node. AIDD is a universal
framework because it can be applied to various types of dynamics, including
continuous, discrete, and binary. AIDD is robust not only on noisy input
signals, but also on unobservable nodes. AIDD can recover an entire network
even when time series data is missing with more than 90% accuracy. AIDD was
also shown to work well on datasets generated by GeneNetWeaver, which emulates
the real environment of gene regulatory network dynamics.

Furthermore, we propose a new method based on controls to test the validity of
the learned model. We have optimised two controllers, one of which is
optimised on the real system and the other is optimised on the model. Then, we
apply them to a real system. If they behave similarly, we conclude that the
learned model can be a substitute for the real system. Control experiments on
spring and CMN dynamics of small networks have proved that well-trained AIDD
models can replace the real systems.

This framework has many potential applications. For example, AIDD can be used
to infer missing links according to the dynamics information. AIDD can also be
used in time series forecasting. In contrast to other forecasting models, a
clear binary network can be output by AIDD, which can provide deeper insights
into element interactions and potential causal links, increasing the
explanability of the model.

However, some drawbacks are present in AIDD. First, a large amount of training
data, especially the time series in diverse initial conditions, is required to
obtain a good model. Nevertheless, it is difficult to obtain different time
series under a variety of initial conditions. Although we can split a long
time series into segments, and the first values on each segment can be treated
as a new initial condition, the diversity of the initial conditions is always
not high enough to train an AIDD model with high quality. Hence we may develop
new models that are suitable for small data.

Second, all the dynamics considered in this paper are Markovian, but this
property is hardly satisfied in real cases. New extensions and experiments on
non-Markovian dynamics should be conducted. For example, we can use a
recurrent neural network instead of a feed-forward network as the dynamics
learning component.

Finally, our network generator samples networks according to the naive mean
field assumption. Although good results have been obtained on network
inference, correlations between nodes are ignored. Thus, we can use generative
graph models to replace the Bernoulli network generator such that correlations
and inductive bias on structures can be considered.

## Methods

### AIDD

The framework consists of two parts: a network generator and a dynamics
learner. The input of the model is the state information of all nodes at time
$t$, and the output of the model is the predicted state information of all
nodes at time $t+1$. The inferred adjacency matrix $\hat{A}$ can also be
retrieved from the network generator.

The network generator is simply a differential matrix sampler parameterized by
$N^{2}$ parameters $\Theta$, as described in the model sub-section of the main
text. However, the dynamics learners are different for different tasks. We
will illustrate the details as follows.

#### Architecture of the Dynamics Learning module

In all experiments except the final one (the gene regulatory network inference
from GeneNetWeaver generated data), we used the same neural network structure
to construct the dynamics learning module, as shown in Fig. 7. The network
structure and weights are shared by all nodes. This design has been verified
to be suitable for learning various complex dynamics [19, 24].

In all experiments, we set the hidden layer size of the dynamics learner to
128. The parameter $\tau$ in the network generator was set to 1. In the
experiments, we set $K$ to be the number of epochs because $\hat{A}_{\cdot
i}(\Theta)$ are sampled at each epoch.

The $\lambda$ in the objective function is the coefficient of the sparse
matrix. In the experiments where the number of nodes is less than $1000$,
$\lambda$ is set to $0.0001$, while $\lambda$ is set to $0.001$ for larger
networks. More details about the parameters can be found in the Supplementary
Section 9.

Figure 7: The feedforward process of the dynamics learning module. It can be
divided into four steps. (1) Node to Edge: aggregating the original
information of nodes to form representations of edges; (2) Edge to Edge:
update the edge representations; (3) Edge to Node: aggregate all information
on neighbouring edges of each node to form a new feature vector of the current
node; (4) node to node: update the node representations; (5) Output: finally,
concatenate the node representations and the input state vectors of node $i$
to feed into a feedforward network, and output the prediction of the next
state of $i$.

#### Architecture of AIDD for gene regulatory network inference

We adopt a heterogeneous network structure in the task of gene regulatory
network inference because the regulation dynamics and parameters are
heterogeneous for different genes. In detail, the dynamics learning module
consists of several multi-layer perceptrons (MLPs), with each MLP
corresponding to a single gene. The input of the MLP is a vector produced by
the element-wise product of the gene expression vector $\mathbf{X}^{t}$ and
the column vector $\hat{A}_{i}$ of the adjacency matrix, which represents the
TF(transcription factors) regulation acting on the corresponding gene $i$. The
output of the MLP is $\hat{X}_{i}^{t+1}$ , which is the estimated
concentration of gene $i$ at time $t+1$. The concatenation of all MLP outputs
is the gene expression value $\hat{\mathbf{X}}^{t+1}$ at time $t+1$. Then, we
can compare the output estimation $\hat{\mathbf{X}}^{t+1}$ and the real
expression value $\mathbf{X}^{t+1}$ to compute the loss function. Then, the
stochastic gradient descent algorithm can be applied to minimise the loss
function. The structure diagram of the model is shown in Fig. 8.

Figure 8: The architecture of AIDD in the task of gene regulatory network
inference. The same network generator is used, but different dynamics learners
are set for each node. In the forward process, the adjacency matrix is
generated by the network generator through Gumbel softmax sampling. Then, the
element-wise products between gene expression vector $X^{t}$ at time $t$ and
column $i$ of the adjacency matrix are calculated as the input for the
corresponding dynamics learner $i$. Subsequently, the dynamics learner $i$,
which is an MLP, computes the output $\hat{X}_{i}^{t+1}$, which is the
estimation of gene i’s expression at time $t+1$. The back-propagation process
updates the network generator and dynamics learners simultaneously.

### Training

We separate the training data by batches and use the stochastic gradient
descent algorithm to update the parameters $\Theta,\Phi$ step by step with
different learning rates $lr_{\theta}$ and $lr_{\phi}$, respectively, until
the epochs of training exceed a threshold. The gradients can be computed
directly by the automatic differentiation technique on the Pytorch platform
because all the steps in our framework are differentiable. To improve the
efficiency, we sampled one adjacency matrix at each epoch and updated the
parameters immediately. For the complete algorithm, readers are referred to
the Supplementary Section 1. We implemented the gradient descent algorithm
with the Adam optimizer, and all the algorithms were run on an RTX
2080Ti(11G).

### Network Inference and Time Series Forecasting

After sufficient training, the optimised parameters of $\Theta^{*}$ and
$\Phi^{*}$ can be obtained. Then, we can use a network generator to sample the
complete adjacency matrices as the inferred networks by setting
$\Theta=\Theta^{*}$ and $\tau\rightarrow\infty$ to obtain absolute 0 or 1.

A single time-step prediction of states for all nodes can be sampled by using
the dynamics learner $\hat{\mathbf{x}}^{t^{\prime}+1}\sim
P(\hat{\mathbf{X}}^{t^{\prime}+1}|\mathbf{x}^{t^{\prime}},\hat{A}(\Theta^{*}),\Phi^{*})$,
where $\hat{A}(\Theta^{*})$ is a sampled adjacency matrix by the optimal
parameters. Multiple time-step predictions can also be obtained in an
independent rollout manner [70], that is, sampling the state recursively,
$\hat{\mathbf{x}}^{t^{\prime}+1}\sim
P(\hat{\mathbf{X}}^{t^{\prime}+1}|\hat{\mathbf{x}}^{t^{\prime}},\hat{A}(\Theta^{*}),\Phi^{*})$
for all $t^{\prime}>0$. Note that $\hat{\mathbf{x}}^{t^{\prime}}$ represents
the sample of $\hat{\mathbf{X}}^{t^{\prime}}$ according to the dynamics
learner.

Areas under the curve (AUC) and mean square errors (MSE for continuous values)
or AUC (for binary values) are used to evaluate the results of network
inference and time series forecasting, respectively. More details on training
and evaluation can be found in the Supplementary Section 8 and 9.

## Data availability

We used the default parameters of DREAM4_In-Silico in GeneNetWeaver software
to generate gene network, GeneNetWeaver software can be download on
https://github.com/tschaffter/gnw [67]. Three social networks (email, Dorm,
and Blog), and a road network(Road) [68] can be found at
http://networkrepository.com/email-enron-only.php#collapse_ack.

## Code availability

AIDD code repository can be found at https://github.com/kby24/AIDD. The
repository includes example codes for generating data and Algorithms of AIDD.

## References

  * [1] Siegenfeld, A. F. & Bar-Yam, Y. An introduction to complex systems science and its applications. _Complexity_ 2020 (2020).
  * [2] Boccaletti, S., Latora, V., Moreno, Y., Chavez, M. & Hwang, D.-U. Complex networks: Structure and dynamics. _Phys.reports_ 424, 175–308 (2006).
  * [3] Bullmore, E. & Sporns, O. Complex brain networks: graph theoretical analysis of structural and functional systems. _Nat.reviews neuroscience_ 10, 186–198 (2009).
  * [4] Watts, D. J. & Strogatz, S. H. Collective dynamics of ’small world’ networks. _Nature_ 393, 440–442 (1998).
  * [5] Runge, J., Nowack, P., Kretschmer, M., Flaxman, S. & Sejdinovic, D. Detecting and quantifying causal associations in large nonlinear time series datasets. _Sci. Adv._ 5, eaau4996 (2019).
  * [6] Liu, Y.-Y. & Barabási, A.-L. Control principles of complex systems. _Rev. Mod. Phys._ 88, 035006 (2016).
  * [7] Wang, W.-X., Lai, Y.-C. & Grebogi, C. Data based identification and prediction of nonlinear and complex dynamical systems. _Phys. Reports_ 644, 1–76 (2016).
  * [8] Ha, S. & Jeong, H. Deep learning reveals hidden interactions in complex systems. Preprint at https://arxiv.org/abs/2001.02539 (2020).
  * [9] Brockwell, P. J., Davis, R. A. & Fienberg, S. E. _Time Series: Theory and Methods_ (Springer Science & Business Media, 1991).
  * [10] Benidis, K. et al. Neural forecasting: Introduction and literature overview. Preprint at https://arxiv.org/abs/2004.10240 (2020).
  * [11] Scarselli, F., Gori, M., Tsoi, A. C., Hagenbuchner, M. & Monfardini, G. The graph neural network model. _IEEE Transactions on Neural Networks_ 20, 61–80 (2008).
  * [12] Battaglia, P. W. et al. Relational inductive biases, deep learning, and graph networks. Preprint at https://arxiv.org/abs/1806.01261 (2018).
  * [13] Wu, Z. et al. A comprehensive survey on graph neural networks. _IEEE Transactions on Neural Networks Learn. Syst._ 32, 4–24 (2021).
  * [14] Zhou, J. et al. Graph neural networks: A review of methods and applications. Preprint at https://arxiv.org/abs/1812.08434 (2018).
  * [15] Zhang, Z., Cui, P. & Zhu, W. Deep learning on graphs: A survey. _IEEE Transactions on Knowl. Data Eng._ 1–1 (2020).
  * [16] Yu, B., Yin, H. & Zhu, Z. Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting. In _Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence IJCAI_ , 3634–3640 (ijcai.org,2018).
  * [17] Sanchez-Gonzalez, A. et al. Graph networks as learnable physics engines for inference and control. In _International Conference on Machine Learning_ , 4470–4479 (PMLR, 2018).
  * [18] Bapst, V. et al. Unveiling the predictive power of static structure in glassy systems. _Nat. Phys._ 16, 448–454 (2020).
  * [19] Kipf, T., Fetaya, E., Wang, K.-C., Welling, M. & Zemel, R. Neural relational inference for interacting systems. In _International Conference on Machine Learning_ , 2688–2697 (PMLR, 2018).
  * [20] Zhang, Z. et al. Neural gene network constructor: A neural based model for reconstructing gene regulatory network. Preprint at https://www.biorxiv.org/content/10.1101/842369v1.abstract (2019).
  * [21] Zheng, C., Fan, X., Wang, C. & Qi, J. Gman: A graph multi-attention network for traffic prediction. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , 34, 1234–1241 (2020).
  * [22] Guo, S., Lin, Y., Feng, N., Song, C. & Wan, H. Attention based spatial-temporal graph convolutional networks for traffic flow forecasting. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , 33, 922–929 (2019).
  * [23] Alet, F., Weng, E., Lozano-Pérez, T. & Kaelbling, L. P. Neural relational inference with fast modular meta-learning. _Adv.Neural Inf. Process. Syst._ 32, 11827–11838 (2019).
  * [24] Zhang, Z. et al. A general deep learning framework for network reconstruction and dynamics learning. _Appl. Netw. Sci._ 4, 1–17 (2019).
  * [25] Pareja, A. et al. Evolvegcn: Evolving graph convolutional networks for dynamic graphs. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , 34, 5363–5370 (2020).
  * [26] Franceschi, L., Niepert, M., Pontil, M. & He, X. Learning discrete structures for graph neural networks. In _International conference on machine learning_ , 1972–1982 (PMLR, 2019).
  * [27] Pearl, J. _Causality_ (Cambridge university press, 2009).
  * [28] Tank, A., Covert, I., Foti, N., Shojaie, A. & Fox, E. Neural granger causality for nonlinear time series. Preprint at https://arxiv.org/abs/1802.05842 (2018).
  * [29] Löwe, S., Madras, D., Zemel, R. & Welling, M. Amortized causal discovery: Learning to infer causal graphs from time-series data. Preprint at https://arxiv.org/abs/2006.10833 (2020).
  * [30] Glymour, C., Zhang, K. & Spirtes, P. Review of causal discovery methods based on graphical models. _Front. genetics_ 10, 524 (2019).
  * [31] Peng, J., Wang, P., Zhou, N. & Zhu, J. Partial correlation estimation by joint sparse regression models. _J. Am. Stat. Assoc._ 104, 735–746 (2009).
  * [32] Stuart, J. M., Segal, E., Koller, D. & Kim, S. K. A gene-coexpression network for global discovery of conserved genetic modules. _Science_ 302, 249–255 (2003).
  * [33] Stetter, O., Battaglia, D., Soriano, J. & Geisel, T. Model-free reconstruction of excitatory neuronal connectivity from calcium imaging signals. _PLoS Comput. Biol_ 8, e1002653 (2012).
  * [34] Nitzan, M., Casadiego, J. & Timme, M. Revealing physical interaction networks from statistics of collective dynamics. _Sci.advances_ 3, e1600396 (2017).
  * [35] Timme, M. & Casadiego, J. Revealing networks from dynamics: an introduction. _J. Phys. A: Math. Theor._ 47, 343001(2014).
  * [36] Liu, H., Kim, J. & Shlizerman, E. Functional connectomics from neural dynamics: probabilistic graphical models for neuronal network of caenorhabditis elegans. _Philos. Transactions Royal Soc. B: Biol. Sci._ 373, 20170377 (2018).
  * [37] Runge, J. Causal network reconstruction from time series: From theoretical assumptions to practical estimation. _Chaos:An Interdiscip. J. Nonlinear Sci._ 28, 075310 (2018).
  * [38] Casadiego, J., Nitzan, M., Hallerberg, S. & Timme, M. Model-free inference of direct network interactions from nonlinear collective dynamics. _Nat. communications_ 8, 1–10 (2017).
  * [39] Li, L., Xu, D., Peng, H., Kurths, J. & Yang, Y. Reconstruction of complex network based on the noise via QR decomposition and compressed sensing. _Sci. reports_ 7, 1–13 (2017).
  * [40] Granger, C. W. Investigating causal relations by econometric models and cross-spectral methods. _Econom. journal Econom.Soc._ 424–438 (1969).
  * [41] Sugihara, G. et al. Detecting causality in complex ecosystems. _Science_ 338, 496–500 (2012).
  * [42] Ma, C., Chen, H.-S., Lai, Y.-C. & Zhang, H.-F. Statistical inference approach to structural reconstruction of complex networks from binary time series. _Phys. Rev. E_ 97, 022301 (2018).
  * [43] Pathak, J., Hunt, B., Girvan, M., Lu, Z. & Ott, E. Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach. _Phys. review letters_ 120, 024102 (2018).
  * [44] Velickovic,P. et al. Graph attention networks. In _International Conference on Learning Representations_ (OpenReview.net,2018)
  * [45] Wang, D. et al. Learning attribute-structure co-evolutions in dynamic graphs. Preprint at https://arxiv.org/abs/2007.13004 (2020).
  * [46] Velickovic, P. et al. Pointer graph networks. _Stat_ 1050, 11 (2020).
  * [47] Wu, Z., Pan, S., Long, G., Jiang, J. & Zhang, C. Graph wavenet for deep spatial-temporal graph modeling. In _Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI-19_ , 1907–1913 (IJCAI, 2019).
  * [48] Li, Y., Meng, C., Shahabi, C. & Liu, Y. Structure-informed graph auto-encoder for relational inference and simulation. In _ICML Workshop on Learning and Reasoning with Graph-Structured Data_ 8, (2019).
  * [49] Ayed, I., de Bézenac, E., Pajot, A., Brajard, J. & Gallinari, P. Learning dynamical systems from partial observations. Preprint at https://arxiv.org/abs/1902.11136 (2019).
  * [50] Pearl, J. & Mackenzie, D. _The Book of Why: The New Science of Cause and Effect_ (Basic Books, 2018).
  * [51] Baggio, G., Bassett, D. S. & Pasqualetti, F. Data-driven control of complex networks. Preprint at https://arxiv.org/abs/2003.12189 (2020).
  * [52] Baydin, A. G., Barak, P., Radul, A. A. & Siskind, J. Automatic differentiation in machine learning: a survey. _J. Mach.Learn. Res._ 18, 1–43 (2018).
  * [53] Gardiner, C. W. et al. _Handbook of Stochastic Methods_ 3 (springer Berlin, 1985).
  * [54] Hochreiter, S. & Schmidhuber, J. Long short-term memory. _Neural computation_ 9, 1735–1780 (1997).
  * [55] Schölkopf, B. Causality for machine learning. Preprint at https://arxiv.org/abs/1911.10500 (2019).
  * [56] Jang, E., Gu, S. & Poole, B. Categorical reparameterization with gumbel-softmax. In _Fifth International Conference on Learning Representations, ICLR_ (OpenReview.net, 2017).
  * [57] Doersch, C. Tutorial on variational autoencoders. Preprint at https://arxiv.org/abs/1606.05908 (2016).
  * [58] Lee, J. et al. Deep neural networks as gaussian processes. In _Sixth International Conference on Learning Representations, ICLR_ (OpenReview.net, 2018).
  * [59] Brockmann, D. & Helbing, D. The hidden geometry of complex, network-driven contagion phenomena. _Science_ 342, 1337–1342 (2013).
  * [60] Karlebach, G. & Shamir, R. Modelling and analysis of gene regulatory networks. _Nat. Rev. Mol. Cell Biol._ 9, 770–780(2008).
  * [61] Barzel, B.& Barabási, A.-L. Network link prediction by global silencing of indirect correlations. _Nat. biotechnology_ 31, 720–725 (2013).
  * [62] Garcia, P., Parravano, A., Cosenza, M., Jiménez, J. & Marcano, A. Coupled map networks as communication schemes. _Phys. Rev. E_ 65, 045201 (2002).
  * [63] Jalan, S., Amritkar, R. & Hu, C.-K. Synchronized clusters in coupled map networks. i. numerical studies. _Phys. Rev. E_ 72, 016211 (2005).
  * [64] Campbell, A., Gurin, G. & Miller, W. E. _The Voter Decides_. (Row, Peterson, and Co., 1954).
  * [65] Erdos, P.& Rényi, A. On the evolution of random graphs. _Publ. Math. Inst. Hung. Acad. Sci_ 5, 17–60 (1960).
  * [66] Barabási, A.-L. & Albert, R. Emergence of scaling in random networks. _Science_ 286, 509–512 (1999).
  * [67] Schaffter, T., Marbach, D. & Floreano, D. Genenetweaver: in silico benchmark generation and performance profiling of network inference methods. _Bioinformatics_ 27, 2263–2270 (2011).
  * [68] Rossi, R. & Ahmed, N. The network data repository with interactive graph analytics and visualization. In _Proceedings of the AAAI Conference on Artificial Intelligence_ 29, 4292–4293 (AAAI Press, 2015)
  * [69] Mengyuan Chen. et al. Inference for network structure and dynamics from time series data via graph neural network. Preprint at https://arxiv.org/abs/2001.06576 (2020).
  * [70] Tang, C. & Salakhutdinov, R. R. Multiple futures prediction. In _Advances in Neural Information Processing Systems_ , 15398–15408 (2019).

## Acknowledgements (not compulsory)

We acknowledge Prof. Qinghua Chen, Dr. Lifei Wang, and the workshops in Swarma
Club for the helpful discussion. We acknowledge the support of the National
Natural Science Foundation of China (NSFC) under grant numbers 61673070.

## Author contributions statement

J.Z. conceived the experiments, Y.Z. and Y.G. conducted most of the
experiments. Z.Z. and S.W. conducted the experiments on gene networks, M.Y.C.
conducted the experiments of robustness. J.Z., Y.Z., and Y.G. wrote, and all
authors reviewed the manuscript.

## Competing interests

The authors declare no competing interests.