# Inferring networks from time series:
a neural approach

Thomas Gaskin[1,2,⋆]
Grigorios A. Pavliotis[2]
Mark Girolami[3,4]

###### Abstract

Network structures underlie the dynamics of many complex phenomena, from gene
regulation and foodwebs to power grids and social media. Yet, as they often
cannot be observed directly, their connectivities must be inferred from
observations of their emergent dynamics. In this work we present a powerful
and fast computational method to infer large network adjacency matrices from
time series data using a neural network. Using a neural network provides
uncertainty quantification on the prediction in a manner that reflects both
the non-convexity of the inference problem as well as the noise on the data.
This is useful since network inference problems are typically underdetermined,
and a feature that has hitherto been lacking from network inference methods.
We demonstrate our method’s capabilities by inferring line failure locations
in the British power grid from observations of its response to a power cut.
Since the problem is underdetermined, many classical statistical tools (e.g.
regression) will not be straightforwardly applicable. Our method, in contrast,
provides probability densities on each edge, allowing the use of hypothesis
testing to make meaningful probabilistic statements about the location of the
power cut. We also demonstrate our method’s ability to learn an entire cost
matrix for a non-linear model from a dataset of economic activity in Greater
London. Our method outperforms OLS regression on noisy data in terms of both
speed and prediction accuracy, and scales as $N^{2}$ where OLS is cubic. Since
our technique is not specifically engineered for network inference, it
represents a general parameter estimation scheme that is applicable to any
parameter dimension.

  * 1Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge CB3 0WA, United Kingdom; 2Department of Mathematics, Imperial College London, London SW7 2AZ, United Kingdom; 3Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, United Kingdom; 4The Alan Turing Institute, London NW1 2DB, United Kingdom

⋆To whom correspondence should be addressed<EMAIL_ADDRESS>
###### Contents

  1. I Introduction
  2. II Results
  3. III Discussion
  4. IV Methods

## I. Introduction

Networks are important objects of study across the scientific disciplines.
They materialise as physical connections in the natural world, for instance as
the mycorrhizal connections between fungi and root networks that transport
nutrients and warning signals between plants [1, 2], human traffic networks
[3, 4], or electricity grids [5, 6]. However, they also appear as abstract,
non-physical entities, such as when describing biological interaction networks
and food webs [7, 8, 9], gene or protein networks [10, 11], economic cost
relations [12, 13], or social links between people along which information
(and misinformation) can flow [14, 15, 16]. In all examples, though the links
constituting the network may not be tangible, the mathematical description is
the same. In this work, we are concerned with inferring the structure of a
network from observations of its dynamics. The problem is of great scientific
bearing: for instance, one may wish to understand the topology of an online
social network from observing how information is passed through it, and some
work has been done on this question [17, 18, 19]. Another important
application is inferring the connectivity of neurons in the brain by observing
their responses to external stimuli [20, 21]. In an entirely different
setting, networks crop up in statistics in the form of conditional
independence graphs, describing dependencies between different variables,
which again are to be inferred from data [22, 23].

Our approach allows inferring network connectivities from time series data
with uncertainty quantification. Uncertainty quantification for network
inference is important for two reasons: first, the observations will often be
noisy, and one would like the uncertainty on the data to translate to an
uncertainty on the predicted network. Secondly however, completely inferring
large networks requires equally large amounts of data – typically at least
$N-1$ equations per node, $N$ being the number of nodes – and these
observations must furthermore be linearly independent. Data of such quality
and quantity will often not be available, leading to an underdetermined
inference problem. The uncertainty on the predicted network should thus also
reflect (at least to a certain degree) the non-convexity of the problem under
consideration: how many networks are compatible with the observed data? To the
best of our knowledge, no current network inference method is able to provide
this information.

Network inference can be performed using ordinary least squares (OLS)
regression [6, 24], but this is confined to the case where the dynamics are
linear in the adjacency matrix. In addition, without additional constraints
OLS usually only works when the network is uniquely identifiable, and its
computational cost typically scales as $N^{3}$, making it infeasible for
larger networks. Our method avoids these limitations. Computationally fast
network inference methods have been developed [17, 18, 19], but these tend to
be highly specialised to a particular type of observation data and give no
uncertainty quantification on the network prediction. Our method, by contrast,
is versatile, not having been specifically engineered to fit the network case,
and having previously been applied to learn low-dimensional vectors of
parameters for stochastic differential equations [25]. The presented approach
is thus not limited to networks, but constitutes a general parameter
estimation method.

In the following, we train a neural network to find a network adjacency matrix
$\mathbf{A}\in\mathbb{R}^{N\times N}$ that, when inserted into an underlying
physical equation, optimally reproduces a given dataset. Using a neural net
allows us to exploit the fact that, as the net trains, it traverses the space
of all graphs, calculating a loss function at each point. By tracking the its
path and gathering the loss values we are able to construct a probability
density on the adjacency matrix, allowing for uncertainty quantification. In
this work we learn networks with up to $1000$ nodes (or a million parameters
to be learned) while outperforming traditional regression methods on noisy
data, both in terms of predictive performance and computational efficiency. We
begin with two application studies: first, we infer locations of a power line
cut in the British power grid from observations of the network response to the
cut; and secondly, we infer economic cost relations between retail centres in
Greater London. Thereafter we conduct a comparative analysis of our method’s
performance, before finally analysing the relationship between uncertainty in
the data, the convexity of the problem, and the prediction uncertainty.

(a) Power grid topology

(b) Network response to line failure

Figure 1: (a) Approximate high-voltage electricity transmission grid of Great
Britain. Shown are 630 accurately placed nodes, representing power stations,
substations, and transmission line intersections, and their connectivity as of
January 2023 [26, 27, 28]. Colours indicate the operating voltage of the
lines. The size of the nodes indicate their power generation or consumption
capacity (absolute values shown). White ringed nodes indicate the 38 nodes
that are real power stations with capacities over 400 MW [29], with all other
nodes assigned a random capacity in $[-200,+200]$. The two dotted edges in the
northeast of England are the edges affected by a simulated power cut, labelled
by the indices of their start and end vertices. (b) The network response to
the simulated power line failure, measured at four different nodes in the
network (marked A–D). The equation parameters were tuned to ensure phase-
locking of the oscillators ($\alpha=1$, $\beta=0.2$, $\kappa=30$). Nodes
closer to the location of the line cut (A and B) show a stronger and more
immediate response than nodes further away (C and D). The shaded area
indicates the 4-second window we use to infer the line location.

(a) Densities on edges with highest relative error

(b) True (black) and predicted network response

Figure 2: Estimating the line failure location. (a) The densities on four
edges with the highest relative prediction error
$|\hat{a}_{ij}-a^{0}_{ij}|/a^{0}_{ij}$ and their respective $p$-values for
measuring the unperturbed value $a^{0}_{ij}$. Red dotted lines indicate the
values of the unperturbed network, green lines the expectation values of the
distributions. The marginals are smoothed using a Gaussian kernel. We use a
training set of length $L=400$ steps, and the batch size is $B=2$. CPU
runtime: 24 minutes. (b) True (black) and predicted network responses at six
different locations in the network. The responses are each normalised to the
value at $t=0$. The shaded area represents the 400 time steps used to train
the model. While the model is able to perfectly fit the response within the
training range, it is not able to learn the full network from insufficient
data, causing the time series to diverge for larger $t$.

## II. Results

### Inferring line failures in the British power grid

Power grids can be modelled as networks of coupled oscillators using the
Kuramoto model [30, 31, 32, 33, 34]. Each node $i$ in the network either
produces or consumes electrical power $P_{i}$ while oscillating at the grid
reference frequency $\Omega$. The nodes are connected through a weighted
undirected network $\mathbf{A}=(a_{ij})$, where the link weights $a_{ij}\sim
Y_{ij}U_{ij}^{2}$ are obtained from the electrical admittances $Y_{ij}$ and
the voltages $U_{ij}$ of the lines (see Methods). This coupling allows the
phases $\varphi_{i}(t)$ of the nodes to synchronise according to the
differential equation [33]

$\alpha\dfrac{\mathrm{d}^{2}\varphi_{i}}{\mathrm{d}t^{2}}+\beta\dfrac{\mathrm{d}\varphi_{i}}{\mathrm{d}t}=P_{i}+\kappa\sum_{j}a_{ij}\sin(\varphi_{j}-\varphi_{i}),$
(1)

where $\alpha$, $\beta$, and $\kappa$ are the inertia, friction, and coupling
coefficients respectively. A requirement for dynamical stability of the grid
is that $\sum_{i}P_{i}=0$, i.e. that as much power is put into the grid as is
taken out through consumption and energy dissipation [32].

A power line failure causes the network to redistribute the power loads,
causing an adjustment cascade to ripple through the network until equilibrium
is restored [5]. In this work we recover the location of a line failure in the
British power grid from observing these response dynamics. Figure 1a shows the
high-voltage transmission grid of Great Britain as of January 2023, totalling
630 nodes (representing power stations, substations, and line intersections)
and 763 edges with their operating voltages. Of the roughly 1300 power
stations dotted around the island, we include those 38 with installed
capacities of at least 400 MW that are directly connected to the national grid
[29]; following [5, 32] we give all other nodes a random value
$P_{i}\sim\mathcal{U}[-200,+200]$ such that $\sum_{i}P_{i}=0$. We simulate a
power cut in the northeast of England in the following way: first, we let the
unperturbed network with adjacency matrix $\mathbf{A}^{0}$ reach a steady
state of phase synchronisation, before removing two links and recording the
network response (fig. 1b). From the response we can infer the adjacency
matrix of the perturbed network $\mathbf{\tilde{A}}$ (with missing links) and,
by comparing with $\mathbf{A}^{0}$, the line failure locations.

We let a neural network output a (vectorised) adjacency matrix
$\mathbf{\hat{A}}$ and use this estimated adjacency matrix to run the
differential equation eq. [1], which will produce an estimate of the phase
vector $\bm{\hat{\varphi}}$. A loss function $J(\mathbf{\hat{A}})$ is used to
train the internal weights $\bm{\theta}$ of the neural network such that it
will output an adjacency matrix that reproduces the observed data. It includes
the prediction error $\|\bm{\hat{\varphi}}-\bm{\varphi}\|_{2}$, two terms to
enforce undirectedness and tracelessness, and a regularising expression to
ensure the prediction $\mathbf{\hat{A}}$ is ‘close to’ the unperturbed
network:

$\displaystyle J_{\theta}(\mathbf{\hat{A}})$
$\displaystyle=\|\bm{\hat{\varphi}}(\mathbf{\hat{A}})-\bm{\varphi}\|_{2}+\|\mathbf{\hat{A}}-\mathbf{\hat{A}}^{T}\|_{2}+\mathrm{tr}(\mathbf{\hat{A}})$
$\displaystyle+\nu\|\mathbf{\hat{A}}-\mathbf{A}^{0}\|_{2}.$ (2)

The $\theta$ subscript indicates that the loss function is conditioned on the
inner weights and biases of the neural net, which are optimised using the
gradient $\partial_{\theta}J$. The
$\nu\|\mathbf{\hat{A}}-\mathbf{A}^{0}\|_{2}$ term is a function that initially
pushes the prediction towards $\mathbf{A}^{0}$ (see Methods). The training
process produces many estimates $\mathbf{\hat{A}}$ of the network, and from
these a probability distribution $\rho(\mathbf{A})$ can be estimated; see
Methods for details.

The time series $\bm{\varphi}(t)$ contains $L$ observations for each node, and
we perform a gradient descent step after predicting $B$ steps, $B$ being the
batch size. In theory $L=N-1$ is needed to completely infer the network,
though symmetries in the data usually mean $L>N$ observations are required in
practice [35]. Traditional regression methods are not applicable in the
underdetermined case $L<N-1$, since they involve inversion of the Gram matrix
$\mathbf{G}_{i}\mathbf{G}_{i}^{T}$, $\mathbf{G}_{i}\in\mathbb{R}^{N\times L}$
containing all $L$ interactions of the $i$-th node with all other nodes, and
additional constraints (e.g. sparsity assumptions) must thus be placed on the
network to infer its topology [36]. We purposefully underdetermine the problem
by only using $L<N-1$ steps, leading to an average Gram matrix rank of 6
(about 100 times less than required for complete inference). Additionally, we
train the network on data recorded $1$ second after the power cut, where many
nodes will still be close to equilibrium. Our method does not involve matrix
inversion, and can thus deal naturally with non-convex problems. Though the
neural network may be unable to completely infer the network, it can
nevertheless produce a density $\rho(\mathbf{\hat{A}})$, recorded during the
training, that allows us to perform hypothesis testing on the line failure
location.

We show the results in fig. 2, where we plot the densities on the four network
edges with the highest relative prediction error
$|\hat{a}_{ij}-a^{0}_{ij}|/a^{0}_{ij}$, $\hat{a}_{ij}$ being the most likely
value. The advantage of obtaining uncertainty quantification on the network is
now immediately clear: even in the underdetermined case we are able to make
meaningful statistical statements about the line failure location. We see that
the missing edges consistently have the highest relative prediction errors,
and that the $p$-values for measuring the unperturbed value $a^{0}_{ij}$ are
around 5%, while being statistically insignificant for all other edges. It is
interesting to note that the other candidate locations are also within the
vicinity of the line failure, though their predicted values are much closer to
the unperturbed value. In fig. 2b, we see that the predicted network
reproduces the response dynamics for the range covered by the training data
when inserted into eq. [1], but, since the problem was purposefully
underdetermined, the errors in the prediction $\mathbf{\hat{A}}$ cause the
predicted and true time series to diverge for larger $t$. Densities were
obtained in about twenty minutes on a regular laptop CPU.

### Inferring economic cost networks from noisy data

(a) Model dynamics

(b) Initial data and travel times network

(c) Inferred weighted degree distribution

Figure 3: Inferring economic cost networks. (a) In the model, $N$ origin zones
(red) are connected to $M$ destination zones (blue) through a weighted
directed network. Economic demand flows from the origin zones to the
destination zones, which supply the demand. We model the origin zones
$O_{i}(t)$ as a Wiener process with variance $\sigma_{O}=0.1$. The resulting
cumulative demand at destination zone $j$ is given by $W_{j}$. Note that the
origin zone sizes fluctuate more rapidly than the destination zones, since
there is a delay in the destination zones’ response to changing consumer
patterns, controlled by the parameter $\epsilon$. We use the parameters as
estimated in [25], $\alpha=0.92$, $\beta=0.54$, $\kappa=8.3$, and set
$\epsilon=2$. (b) The initial origin and destination zone sizes, given by the
total household income of the $N=629$ wards in London (blue nodes) and the
retail floor space of $M=49$ major centres (red nodes) [37, 38]. As derived in
[39], the weight on each link is then $\exp(-d_{ij}/\tau)$, where we extract
the transport times $d_{ij}$ from the GoogleMaps API as in [25], and the time
scale $\tau=\max_{ij}d_{ij}$ (see Methods). Background map: [40]. (c)
Predicted degree distribution (sold line) of the inferred network, for a high
noise level of $\sigma=0.14$, with Hellinger uncertainty eq. [7] (shaded
area), and the true distribution (red dotted line). CPU runtime: 3 min 41 s.

In the previous example the underlying network was a physical entity, but in
many cases networks model abstract connections. We therefore now consider a
commonly used economic model of the coupling of supply and demand [41, 13, 12]
and a dataset of economic activity across Greater London. The goal is to learn
the entire coupling network, not just to infer the (non-)existence of
individual edges. In the model, $N$ origin zones, representing economic
demand, are coupled to $M$ destination zones, modelling the supply side,
through a network whose weights quantify the convenience with which demand
from zone $i$ can supplied from zone $j$: the higher the weight, the more
demand flows through that edge (see fig. 3a). Such a model is applicable e.g.
to an urban setting [12], the origin zones representing residential areas, the
destination zones e.g. commercial centres, and the weights quantifying the
connectivity between the two (transport times, distances, etc.). The resulting
cumulative demand at destination zone $j$ depends both on the current size
$W_{j}(t)$ of the destination zone and the network weights $c_{ij}$:

$D_{j}=\sum_{i=1}^{N}\dfrac{W_{j}(t)^{\alpha}c_{ij}^{\beta}}{\sum_{k=1}^{M}W_{k}(t)^{\alpha}c_{ik}^{\beta}}O_{i}(t).$
(3)

The sizes $W_{j}$ are governed by a system of $M$ coupled logistic
Stratonovich stochastic differential equations:

$\mathrm{d}W_{j}=\epsilon W_{j}(D_{j}-\kappa W_{j})\mathrm{d}t+\sigma
W_{j}\circ\mathrm{d}B_{j},$ (4)

with given initial conditions $W_{j}(0)$, see fig. 3a. $\alpha$, $\beta$,
$\kappa$, and $\epsilon$ are scalar parameters. Our goal is to infer the
connectivities $c_{ij}$ from observations of the time series $\mathbf{O}(t)$
and $\mathbf{W}(t)$. The model includes multiplicative noise $B_{j}$ with
variance $\sigma\geq 0$, with $\circ$ signifying Stratonovich integration.
Crucially, the model depends non-linearly on the network $\mathbf{C}$:
classical regression methods are thus not applicable here.

We apply this model to a previously studied dataset of economic activity in
Greater London [13, 25]. We use the ward-level household income from $N=625$
wards for 2015 [38] and the retail floor space of the $M=49$ largest
commercial centres in Greater London as the initial origin zone and
destination zone sizes respectively, i.e. $\mathbf{O}(0)$ and $\mathbf{W}(0)$,
and from this generate a synthetic time series using the parameters estimated
in [25] for a high noise level of $\sigma=0.14$. The network consists of
public transport and driving travel times between the zones, extracted using
the Google Maps Distance API, for June 2022 [25]; see fig. 3b and Methods. We
generate enough data to make the problem uniquely inferable, in order to
isolate the effect of the measurement noise on the prediction.

Figure 3c shows the inferred distribution $P(k)$ of the (weighted) origin zone
node degrees $k_{i}=\sum_{j}c_{ij}$. The solid line is the maximum likelihood
prediction, and the dotted red line the true distribution . Even with a high
level of noise, the model manages to perfectly predict the underlying
connectivity matrix, comprising over 30.000 weights, in under 5 minutes on a
regular laptop CPU. We quantify the uncertainty on the distribution using the
Hellinger metric eq. [7]. As we will discuss in the last section, this method
meaningfully captures the uncertainty due to the noise in the data.

(a) $L^{1}$ error, first-order dynamics

(b) $L^{1}$ error, second-order dynamics

(c) Compute times

(d) Degree distribution, $N=1000$

(e) Triangle distribution, $N=1000$

Figure 4: Computational performance analysis. (a) $L^{1}$ prediction error of
the neural scheme (black) and OLS as a function of the noise variance $\sigma$
on the training data. For very high noise levels, the training data is
essentially pure noise, and the prediction errors begin to plateau. First-
order Kuramoto dynamics are used. (b) The same plot as in (a), but using
second-order dynamics. (c) Compute times for a single epoch of the neural
method as a function of the network size $N$. Shown are the compute times on a
standard CPU (Apple M1, green line) and GPU (NVidia GeForce RTX 3090, pink
line), averaged over 100 runs, with the shaded areas showing one standard
deviation. Also shown is the compute time for ordinary least squares (black
dotted line), which scales as $N^{3}$ (note the logarithmic axis). On the
right axis, the average $L^{1}$ prediction error
$\frac{1}{N}\|\mathbf{\hat{A}}-\mathbf{A}\|_{1}$ of the neural scheme after 10
epochs is shown, which remains fairly constant as a function of $N$, showing
that the number of gradient descent steps required to achieve a given average
prediction error does not depend on $N$. (d) Predicted degree distribution and
(e) triangle distribution of an inferred network, trained on first-order noisy
Kuramoto data ($\sigma=0.001$). The blue shaded areas indicate the Hellinger
uncertainty eq. [7], and the red dotted lines are the true distributions. CPU
runtime: 1 hour 3 minutes.

### Performance analysis and comparison with OLS

We now analyse our method’s performance, both in terms of prediction quality
and computational speed, by comparing it to that of a classical regression
method, presented e.g. in [36, 6]. As mentioned in the introdution, more
efficient network learning methodologies have been developed for specific
problems; however, we compare our method with OLS since both methods are
general and do not rely on a specific data structure. Figures 4a–b show our
method’s prediction accuracy alongside that of an OLS regression scheme
(outlined in Methods) as a function of the noise $\sigma$ on the training
data; the accuracy here is defined as the $L^{1}$ error

$\|\mathbf{\hat{A}}-\mathbf{A}\|_{1}=\sum_{i,j}|\hat{a}_{ij}-a_{ij}|.$ (5)

For the practically noiseless case of $\sigma<10^{-5}$, the regression scheme
on average outperforms the neural approach; however, even for very low noise
levels $\sigma\geq 10^{-5}$ and above, the neural approach proves far more
robust, outperforming OLS by up to one order of magnitude and maintaining its
prediction performance up to low noise levels of $\sigma\leq 10^{-3}$. These
results hold both for first-order ($\alpha=0$) and second-order Kuramoto
equations [1] (see fig. 4b); in the second-order case, the neural method
begins outperforming OLS at even lower levels of $\sigma$ than in the first-
order case, though the improvement is not as significant.

(a) Convexity uncertainty on $P(k)$

(b) Noise uncertainty on $P(k)$

(c) Noise uncertainty on $P(t)$

Figure 5: Quantifying the two types of uncertainty: (a) Total Hellinger error
on the degree distribution $P(k)$ as a function of $\mathfrak{c}$ (eq. [6]) in
the noiseless case. The error is normalised to the value at
$\mathfrak{c}=0.21(N-1)$. As $\mathfrak{c}$ increases, the error on the
prediction decreases almost linearly. We run the model from 10 different
initialisations and average over each (shaded area: standard deviation). (b)
and (c): Prediction uncertainty due to noise in the data. Shown are the
average Hellinger error (eq. [7]) and average relative entropy (eq. [8]) to
the maximum likelihood estimate for the degree distribution $P(k)$ and
triangle distribution $P(t)$ as a function of the noise $\sigma$ in the data.
Each line is an average over 10 different initialisations. In all cases,
training was conducted on synthetic, first-order Kuramoto data (eq. [20] with
$\alpha=0$, $\beta=1$, $\kappa=1$).

In figure 4c we show the average training time per epoch as a function of $N$,
with training conducted both on a standard laptop CPU and a standard GPU. Each
epoch of the model equation requires $\mathcal{O}(LN^{2})$ operations for the
vector-matrix multiplication in eq. [22], and $\mathcal{O}(LN^{2}/B)$ for the
stochastic gradient descent update, where we are holding $L/B$ constant to
ensure comparability. We thus obtain an algorithmic complexity of
$\mathcal{O}(LN^{2})$. As is visible, the average $L^{1}$ error per edge
weight remains constant over $N$, indicating that the number of gradient
descent steps required to achieve a given node-averaged prediction accuracy is
independent of $N$. The total complexity of our method is thus
$\mathcal{O}(n_{E}\times LN^{2})$, with $n_{E}$ the number of training epochs.
By comparison, OLS in general has a cubic dependency on $N$, due to the
required matrix inversion. Being a neural approach, our method can make use of
GPU accelerated training, leading to an order of magnitude faster performance.

Lastly, figures 4d–e show the estimated weighted degree and triangle
distributions of a graph with 1000 nodes, or 1 million edge weights to be
estimated, for noisy training data. The number of weighted, undirected
triangles on each node $i$ is given by
$\frac{1}{2}\sum_{jk}a_{ij}a_{jk}a_{ki}$. The model robustly finds the true
adjacency matrix, and we again quantify uncertainty on the prediction using
the Hellinger distance (eq. [7]) between each distribution estimated during
training and the maximum likelihood estiamte. Estimating a network with 1000
nodes on a standard laptop CPU took about 1 hour. Most high-performance
network inference techniques demonstrate their viability on graphs with at
most this number of nodes, e.g. ConNIe [17] and NetINF [19]. In [17], the
authors state that graphs with 1000 nodes can typically be inferred from
cascade data in under 10 minutes on a standard laptop. Similarly, the authors
of NetINF [19] state that it can infer a network with 1000 nodes in a matter
of minutes, though this algorithm does not infer edge weights, only the
existence of edges, and neither technique provides uncertainty quantification.

### Quantifying uncertainty

There are two sources of uncertainty when inferring adjacency matrices: the
non-convexity of the loss function $J$, and the noise $\sigma$ on the data. In
general, it is hard to quantify the convexity of $J$, since we do not know how
many networks fit the equation at hand. However, when the dynamics are linear
in the adjacency matrix $\mathbf{A}$, we can do so using the Gram matrix of
the observations of each node $i$, $\mathbf{G}_{i}\mathbf{G}_{i}^{T}$. For
regression methods to be applicable, $\mathbf{G}_{i}\mathbf{G}_{i}^{T}$ must
be invertible for each $i$. The (non-)convexity of the problem can thus be
quantified for example by the minimum rank of all the Gram matrices,

$\mathfrak{c}:=\min_{i}\mathrm{rk}\left(\mathbf{G}_{i}\mathbf{G}_{i}^{T}\right).$
(6)

The problem is fully determined if $\mathfrak{c}=N-1$. In figure 5a we show
the Hellinger uncertainty on the predicted degree distribution as a function
of $\mathfrak{c}$. The Hellinger error of the degree distribution is given by

$\displaystyle\bar{d}_{\mathrm{H}}(k)\sim$
$\displaystyle\int\left(\sqrt{\hat{\rho}(k)}-\sqrt{\hat{\rho}\left(k\;\middle|\;\hat{\mathbf{T}}\right)}\right)^{2}\exp(-\|\hat{\mathbf{T}}-\mathbf{T}\|_{2})\mathrm{d}\hat{\mathbf{T}}.$
(7)

Here, we define the maximum likelihood estimate as the predicted density
$\hat{\rho}(k)$ for which uncertainty is to be calculated. As is visible, the
error on the distribution decreases almost linearly as $\mathfrak{c}$ tends to
its maximum value $N-1$. For $\mathfrak{c}=N-1$, some residual uncertainty
remains due to the uncertainty on the neural network parameters $\bm{\theta}$.

In figures 5b–c we again show the Hellinger error on the predicted degree and
triangle distributions, this time as a function of the noise $\sigma$. To
demonstrate that the uncertainty is not an artefact of the choice of metric,
we also show the behaviour of the relative entropy,

$\displaystyle\bar{d}_{\mathrm{KL}}(k)\sim$
$\displaystyle\int\hat{\rho}\left(k\;\middle|\;\hat{\mathbf{T}}\right)\log\left(\dfrac{\hat{\rho}\left(k\;\middle|\;\hat{\mathbf{T}}\right)}{\hat{\rho}(k)}\right)\exp(-\|\hat{\mathbf{T}}-\mathbf{T}\|_{2})\mathrm{d}\hat{\mathbf{T}}.$
(8)

Both metrics reflect the noise on the training data, providing similarly
behaved, meaningful uncertainty quantification. As the noise tends to 0, some
residual uncertainty again remains. Our method thus manages to capture the
uncertainty arising from both sources: the non-convexity of $J$ and the noise
$\sigma$ on the data.

## III. Discussion

In this work we have demonstrated a performative method to estimate network
adjacency matrices from time series data. We showed its effectiveness at
correctly and reliably inferring networks in a variety of scenarios: convex
and non-convex cases, low to high noise regimes, and equations that are both
linear and non-linear in $\mathbf{A}$. We were able to reliably infer power
line failures in the national power grid of Great Britain, and the
connectivity matrix of an economic system covering all of Greater London. We
showed that our method is well able to handle inference of hundreds of
thousands to a million edge weights, while simultaneously giving uncertainty
quantification that meaningfully reflects both the non-convexity of the loss
function as well as the noise on the training data. Our method outperforms
classical ordinary least squares regression, both in terms of prediction
accuracy and computational speed, while extending to the case of non-linear
dynamics. In conjunction with our previous work [25], we have now also
demonstrated the viability of using neural networks for parameter calibration
in both the low- and high-dimensional case. Our method is simple to implement
as well as highly versatile, giving excellent results across a variety of
problems. All experiments in this work were purposefully conducted on a
standard laptop CPU, typically taking on the order of minutes to run; however,
being a neural scheme, our method can make use of GPU acceleration, further
reducing the compute times.

Many lines for future research open up from this work. Firstly, a thorough
theoretical investigation of the method and a more quantitative comparison to
other, similar methods is warranted, e.g. approximate Bayesian computing [42,
43], physics-informed neural networks [44], and using dropout for uncertainty
quantification in neural networks [45]. Another direction is further reducing
the amount of data required to learn parameters, which in many applications
may not be abundantly available, and in future research the authors aim to
address the question of learning systems properties from observations of a
single particle trajectory at the mean-field limit [46, 47].

## IV. Methods

### Neural estimation scheme

We employ the same method as in our previous work [25], used to estimate
parameter densities for ODE or SDE systems. We train a neural net to find a
graph adjacency matrix $\mathbf{\hat{A}}$ that, when inserted into the model
equations, reproduce a given time series
$\mathbf{T}=(\mathbf{x}_{1},...,\mathbf{x}_{L})$. A _neural network_ is a
function $u_{\theta}:\mathbb{R}^{N\times q}\to\mathbb{R}^{p}$, where $q\geq 1$
represents the number of time series steps that are passed as input. Its
output is the (vectorised) estimated adjacency matrix $\mathbf{\hat{A}}$,
which is used to run a numerical solver for $B$ iterations ($B$ is the _batch
size_) to produce an _estimated time series_
$\mathbf{\hat{T}}(\mathbf{\hat{A}})=(\mathbf{\hat{x}}_{i},...,\mathbf{\hat{x}}_{i+B})$.
This in turn is used to train the internal parameters $\bm{\theta}$ of the
neural net (the _weights and biases_) via a _loss functional_
$J(\hat{\mathbf{T}},\mathbf{T})$. As
$\mathbf{\hat{A}}=\mathbf{\hat{A}}(\bm{\theta})$, we may calculate the
gradient $\nabla_{\bm{\theta}}J$ and use it to optimise the internal
parameters of the neural net using a backpropagation method of choice; popular
choices include stochastic gradient descent, Nesterov schemes, or the Adam
optimizer [48]. Calculating $\nabla_{\bm{\theta}}J$ thus requires
differentiating the predicted time series $\mathbf{\hat{T}}$, and thereby the
system equations, with respect to $\mathbf{\hat{A}}$. In other words: the loss
function contains knowledge of the dynamics of the model. Finally, the true
data is once again input to the neural net to produce a new parameter estimate
$\mathbf{\hat{A}}$, and the cycle starts afresh. A single pass over the entire
dataset is called an _epoch_ (see also fig. 6).

The method is in effect a Gibbs sampling scheme, where we wish to obtain the
joint distribution

$p\left(\mathbf{\hat{A}},\mathbf{\hat{T}}\;\middle|\;\mathbf{T}\right)=\int
p\left(\mathbf{\hat{A}},\mathbf{\hat{T}},\bm{\theta}\;\middle|\;\mathbf{T}\right)\mathrm{d}\bm{\theta}$
(9)

by marginalising over the training steps. We sample from this joint density by
sampling the estimated network $\mathbf{\hat{A}}$, the estimated data
$\hat{\mathbf{T}}$, and the neural network weights and biases $\bm{\theta}$
from conditional densities and use each sample for the next conditional,

$\displaystyle\mathbf{\hat{A}}^{i+1}\sim
p\left(\mathbf{\hat{A}}\;\middle|\;\mathbf{\hat{T}}^{i},\bm{\theta}^{i},\mathbf{T}\right)$
(10) $\displaystyle\mathbf{\hat{T}}^{i+1}\sim
p\left(\mathbf{\hat{T}}\;\middle|\;\mathbf{\hat{A}}^{i+1},\bm{\theta}^{i},\mathbf{T}\right)$
(11) $\displaystyle\bm{\theta}^{i+1}\sim
p\left(\bm{\theta}\;\middle|\;\mathbf{\hat{A}}^{i+1},\mathbf{\hat{T}}^{i+1},\mathbf{T}\right).$
(12)

Neural Net Estimated network $\mathbf{\hat{A}}$ Numerical solver Predicted
data $\hat{\mathbf{T}}$ Observed data $\mathbf{T}$ Loss functional $J$
outputsused to runproducesinput toused to train

Figure 6: The methodological pipeline proposed in this work. The neural net
$u_{\theta}$ takes $q$ time series elements as input and outputs a predicted
network adjacency matrix. These predictions are fed into a numerical solver,
which produces a predicted time series. The true and predicted time series are
used to generate a loss functional, which in turn can be used to train the
neural net’s internal parameters $\bm{\theta}$. The goal is to retrieve the
true network $\mathbf{A}$ from the data. A single pass over the entire dataset
is called an _epoch_. The dataset is processed in _batches_ , meaning the loss
is calculated over $B$ steps of the time series before the weights are
updated. If $B=L$, training is equivalent to batch gradient descent; if $B=1$,
training is equivalent to stochastic gradient descent.

Here, $i$ is the training iteration index. Each of the three densities
represent different components in the cycle: eq. [10] is the neural net’s
output, which is deterministic:

$\mathbf{\hat{A}}^{i+1}=u_{\theta^{i}}(\mathbf{\hat{T}}^{i}).$ (13)

Eq. [11] is the output of the numerical solver, which may be stochastic,
depending on whether the numerical solver is run with noise:

$p\left(\mathbf{\hat{T}}\;\middle|\;\mathbf{\hat{A}}^{i+1},\bm{\theta}^{i},\mathbf{T}\right)\sim\rho_{NS}(\mathbf{\hat{A}}^{i+1}).$
(14)

Finally, the updated neural network parameters are drawn using the gradient
descent scheme.

By Bayes’ rule,

$p\left(\mathbf{\hat{A}},\hat{\mathbf{T}},\bm{\theta}\;\middle|\;\mathbf{T}\right)\sim\exp(-\|\hat{\mathbf{T}}-\mathbf{T}\|_{2})\pi^{0}(\mathbf{\hat{A}})\pi^{0}(\bm{\theta})\pi^{0}(\hat{\mathbf{T}}),$

with $\pi^{0}$ the priors [49].

### Simulating power grid dynamics

We simulate the power grid by initialising all phases $\varphi_{i}(t=0)$ at 1
and iterating the Kuramoto dynamics until $|\dot{\varphi}_{i}|/\varphi_{i}\leq
0.01\ \forall i$. We then cut two edges in the network and record the
response. The network weights $a_{ij}$ are given by $U_{ij}^{2}Y_{ij}$, were
$U_{ij}$ is the voltage on the line, and $Y_{ij}$ its electrical admittance in
mho. In general, calculating the admittance of a power line is difficult,
since it depends strongly on the material and geometry of the cable (cross-
section, number of cores, core radii, etc.), whether the line is overground or
underground, and the length of line. Since the neural network outputs are in
$[0,1]$, we scale the network weights $a_{ij}\to\lambda a_{ij}$ such that
$a_{ij}\in[0,1]$, and absorb $\lambda$ into the coupling constant $\kappa$;
see Supplementary Information for details on the calculations.

### Inferring power line failure locations

A hyperparameter sweep on synthetic data showed that using a deep neural
network with 5 layers, 20 nodes per layer, and no bias provides optimal
results. We use the hyperbolic tangent as an activation function on each layer
except the last, where we use the ‘hard sigmoid’

$\sigma(x)=\begin{cases}0,\ x\leq-3,\\\ 1,\ x\geq+3,\\\ x/6+1/2,\
\mathrm{else},\end{cases}\ $ (15)

which allows neural net output components to actually become zero, and not
just asymptotically close, thereby ensuring sparsity of the adjacency matrix.
We use the Adam optimizer [48] with a learning rate of $0.002$ for the
gradient descent step. We initialise the neural network’s weights with a prior
$\pi^{0}(\bm{\theta})$ in such a way that the prior
$\pi^{0}(\mathbf{\hat{A}})$ is a delta distribution on the complete graph,
$\pi(\hat{a}_{ij})\sim\delta(1)\ \forall i,j$. This maximises the sampling
domain on each edge.

Figure 7: The total loss $J$ and its derivatives $\partial_{s}J$ and
$\partial_{ss}J$, averaged over a window of 20 iterations. The red dotted line
indicates the value at which $\nu$ is set to 0.

We rewrite equation [1] as a vector-matrix equation,

$\alpha\dfrac{\mathrm{d}^{2}\bm{\varphi}}{\mathrm{d}t^{2}}+\beta\dfrac{\mathrm{d}\bm{\varphi}}{\mathrm{d}t}=\mathbf{P}+\mathrm{diag}(\mathbf{A}\bm{\Gamma}(\bm{\varphi})),$
(16)

which is more amenable to machine learning purposes, since it can make use of
fast matrix-multiplication operations. $\bm{\Gamma}$ is the interaction
kernel, $\Gamma_{ij}=\sin(\varphi_{j}-\varphi_{i})$. We use the following loss
function:

$\displaystyle J_{\theta}(\mathbf{\hat{A}})=$
$\displaystyle\|\bm{\hat{\varphi}}(\mathbf{\hat{A}})-\bm{\varphi}\|_{2}+\|\mathbf{\hat{A}}-\mathbf{\hat{A}}^{T}\|_{2}+\mathrm{tr}(\mathbf{\hat{A}})$
$\displaystyle+\nu\|\mathbf{\hat{A}}-\mathbf{A}^{0}\|_{2}.$ (17)

The first term is the error on the data, the second penalises asymmetry, and
the third sets the diagonal to zero (which cannot be inferred from the data).
$\nu$ is a function designed to push the neural network in the general
direction of the unperturbed adjacency matrix $\mathbf{A}^{0}$, but become
zero as soon as the neural network reaches $\mathbf{A}^{0}$. We let $\nu=10$
while $|\langle\partial_{s}J\rangle|>10^{-10}$ and
$|\langle\partial_{ss}J\rangle|>10^{-10}$, and $\nu=0$ thereafter, where $s$
is the iteration count (number of training steps), and $\langle J\rangle$ the
total loss averaged over a window of $20$ iterations, see fig. 7. In other
words, we push the neural network towards a stable minimum in the vicinity of
$\mathbf{A}^{0}$ and, once the loss stabilises, permanently set $\nu=0$.

### Inferring economic cost networks

A hyperparameter sweep on synthetic data showed using a neural network with
$2$ layers, 20 nodes per layer, and no bias to yield optimal results. We use
the hyperbolic tangent as the activation function on all layers except the
last, where we use the standard sigmoid function (since the network is
complete, there is no need to use the hard sigmoid, since all edges are
nonzero).

We generate a synthetic time series of length $L=10000$ from which we
subsample $2500$ 2-step windows, giving a total training set size of 5000.
This is to ensure we sample a sufficiently broad spectrum of the system’s
dynamics, thereby fully determining the inference problem. Equation parameters
are as in the high-noise case in [25]: $\alpha=0.92$, $\beta=0.54$,
$\kappa=8.3$, $\sigma=0.14$, and we set $\epsilon=2$. The origin zone and
destination zone sizes are as in [25]. For the network $\mathbf{C}$ we use the
Google Distance Matrix
API111developers.google.com/maps/documentation/distance-matrix to extract
travel data from Google Maps. The API can be used to extract both travel
distances and travel times for large batches of search queries for all the
transport modes available on Google Maps (driving, public transport, walking,
cycling). We define the distance $d_{ij}$ between two nodes as the _shorter_
of the two transport modes considered: for example, travelling from Kentish
Town to the West End takes about 19.7 minutes by public transport, and 22.35
minutes by car, so we set the distance to be 19.7 minutes (of course, this
assumes that everyone can choose between both modes of transport, and does not
factor in the added cost of driving in London.) Note that the Google Maps
public transport mode defaults to walking for very short distances. The
network weights are derived in [39] as

$c_{ij}=e^{-d_{ij}/\tau}$ (18)

where $d_{ij}$ is the distance in the metric in question and
$\tau=\sup_{i,j}d_{ij}$ the time/length scale, ensuring a unitless exponent.

To train the neural network we use the simple loss function

$J_{\theta}(\mathbf{\hat{A}})=\|\bm{\hat{\varphi}}(\mathbf{\hat{A}})-\bm{\varphi}\|_{2}.$
(19)

Since the dynamics are invariant under scaling of
$\mathbf{C}\to\lambda\mathbf{C}$, we normalise the rowsums of the predicted
and true networks, $\sum_{j}c_{ij}=1$.

### Comparison with OLS regression

Consider noisy Kuramoto dynamics,

$\alpha\dfrac{\mathrm{d}^{2}\varphi_{i}}{\mathrm{d}t^{2}}+\beta\dfrac{\mathrm{d}\varphi_{i}}{\mathrm{d}t}=P_{i}+\sum_{j}a_{ij}\sin(\varphi_{j}-\varphi_{i})+\xi_{i},$
(20)

with $\Gamma_{ij}=\sin(\varphi_{j}-\varphi_{i})$,
$\xi_{i}\overset{iid}{\sim}\mathcal{N}(0,\sigma)$, and $\omega_{i}$ the
eigenfrequencies of the nodes. Gathering the $L$ terms on the left side into a
single vector $\mathbf{X}_{i}$ for each node, we obtain $N$ equations

$\mathbf{X}_{i}=\mathbf{A}_{i}\cdot\mathbf{G}_{i}+\epsilon_{i},$ (21)

with $\mathbf{X}_{i}\in\mathbb{R}^{1\times L}$,
$\mathbf{A}_{i}\in\mathbb{R}^{1\times N}$ the $i$-th row of $\mathbf{A}$, and
$\mathbf{G}_{i}\in\mathbb{R}^{N\times L}$ the $L$ observations of the $i$-th
column of $\bm{\Gamma}$. From this we can then naturally estimate the $i$-th
row of $\mathbf{A}$ using ordinary least squares:

$\mathbf{\hat{A}}_{i}=\underset{\bm{\beta}\in\mathbb{R}^{1\times
N}}{\mathrm{argmin}}\
\|\mathbf{X}_{i}-\bm{\beta}\cdot\mathbf{G}_{i}\|_{2}^{2}=\mathbf{X}_{i}\mathbf{G}_{i}^{T}\left(\mathbf{G}_{i}\mathbf{G}_{i}^{T}\right)^{-1}.$
(22)

Given sufficiently many linearly independent observations, the Gram matrix
$\mathbf{G}_{i}\mathbf{G}_{i}^{T}$ will be invertible and the network can be
inferred (with the diagonal manually set to 0). We use synthetic first- and
second order Kuramoto data to compare our method’s performance to that of OLS.
In order to ensure invertibility of $\mathbf{G}_{i}\mathbf{G}_{i}^{T}$, we
generate $L\gtrsim N/2$ datasets of 2 time steps each.

### Code and data availability

Code and synthetic data can be found under
https://github.com/ThGaskin/NeuralABM. It is easily adaptable to new models
and ideas. The code uses the utopya package222utopia-project.org,
utopya.readthedocs.io/en/latest [50, 51] to handle simulation configuration
and efficiently read, write, analyse, and evaluate data. This means that the
model can be run by modifying simple and intuitive configuration files,
without touching code. Multiple training runs and parameter sweeps are
automatically parallelised. The neural core is implemented using
pytorch333pytorch.org. All synthetic datasets as well as the London dataset
have been made available, together with the configuration files needed to
reproduce the plots. Detailed instructions are provided in the supplementary
material and the repository. The British power grid data [26, 27, 28] is
property of the respective organisations and cannot be made available without
permission; however, as of early 2023 it is freely available from those
organisations upon request. The code used to run the experiments is available
in the repository.

Acknowledgements The authors are grateful to Dr Andrew Duncan for the fruitful
discussions on power grid dynamics. TG was funded by the University of
Cambridge School of Physical Sciences VC Award via DAMTP and the Department of
Engineering, and supported by EPSRC grants EP/P020720/2 and EP/R018413/2. The
work of GP was partially funded by EPSRC grant EP/P031587/1, and by J.P.
Morgan Chase & Co through a Faculty Research Award 2019 and 2021. MG was
supported by EPSRC grants EP/T000414/1, EP/R018413/2, EP/P020720/2,
EP/R034710/1, EP/R004889/1, and a Royal Academy of Engineering Research Chair.

## References

  * [1] SW Simard, et al., Mycorrhizal networks: Mechanisms, ecology and modelling.  Fungal Biology Reviews 26, 39–60 (2012).
  * [2] C Hettenhausen, et al., Stem parasitic plant Cuscuta australis (dodder) transfers herbivory-induced signals among plants.  Proceedings of the National Academy of Sciences 114, E6703–E6709 (2017).
  * [3] D Brockmann, D Helbing, The Hidden Geometry of Complex, Network-Driven Contagion Phenomena.  Science 342, 1337–1342 (2013).
  * [4] N Molkenthin, M Schröder, M Timme, Scaling Laws of Collective Ride-Sharing Dynamics.  Physical Review Letters 125 (2020).
  * [5] I Simonsen, L Buzna, K Peters, S Bornholdt, D Helbing, Transient Dynamics Increasing Network Vulnerability to Cascading Failures.  Physical Review Letters 100 (2008).
  * [6] SG Shandilya, M Timme, Inferring network topology from complex dynamics.  New Journal of Physics 13, 013004 (2011).
  * [7] U Stelzl, et al., A Human Protein-Protein Interaction Network: A Resource for Annotating the Proteome.  Cell 122, 957–968 (2005).
  * [8] SR Proulx, DE Promislow, PC Phillips, Network thinking in ecology and evolution.  Trends in Ecology & Evolution 20, 345–353 (2005) Special issue: Bumper book review.
  * [9] S Allesina, D Alonso, M Pascual, A General Model for Food Web Structure.  Science 320, 658–661 (2008).
  * [10] J Tegnér, MKS Yeung, J Hasty, JJ Collins, Reverse engineering gene networks: Integrating genetic perturbations with dynamical modeling.  Proceedings of the National Academy of Sciences 100, 5944–5949 (2003).
  * [11] BO Palsson, Systems Biology: Properties of Reconstructed Networks.  (Cambridge University Press, New York, NY, USA), (2006).
  * [12] M Batty, R Milton, A new framework for very large-scale urban modelling.  Urban Studies 58, 3071–3094 (2021).
  * [13] L Ellam, M Girolami, GA Pavliotis, A Wilson, Stochastic modelling of urban structure.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 474, 20170700 (2018).
  * [14] MD Vicario, et al., The spreading of misinformation online.  Proceedings of the National Academy of Sciences 113, 554–559 (2016).
  * [15] S Aral, L Muchnik, A Sundararajan, Engineering social contagions: Optimal network seeding in the presence of homophily.  Network Science 1, 125–153 (2013).
  * [16] S Vosoughi, D Roy, S Aral, The spread of true and false news online.  Science 359, 1146–1151 (2018).
  * [17] SA Myers, J Leskovec, On the Convexity of Latent Social Network Inference in Proceedings of the 23rd International Conference on Neural Information Processing Systems - Volume 2, NIPS’10.  (Curran Associates Inc., Red Hook, NY, USA), pp. 1741–1749 (2010).
  * [18] M Gomez-Rodriguez, D Balduzzi, B Schölkopf, Uncovering the Temporal Dynamics of Diffusion Networks in Proceedings of the 28th International Conference on International Conference on Machine Learning, ICML’11.  (Omnipress, Madison, WI, USA), pp. 561–568 (2011).
  * [19] M Gomez-Rodriguez, J Leskovec, A Krause, Inferring Networks of Diffusion and Influence.  ACM Transactions on Knowledge Discovery from Data 5, 1–37 (2012).
  * [20] VA Makarov, F Panetsos, O de Feo, A method for determining neural connectivity and inferring the underlying network dynamics using extracellular spike recordings.  Journal of Neuroscience Methods 144, 265–279 (2005).
  * [21] FV Bussel, Inferring synaptic connectivity from spatio-temporal spike patterns.  Frontiers in Computational Neuroscience 5 (2011).
  * [22] N Meinshausen, P Bühlmann, High-dimensional graphs and variable selection with the Lasso.  The Annals of Statistics 34, 1436 – 1462 (2006).
  * [23] M Yuan, Y Lin, Model Selection and Estimation in the Gaussian Graphical Model.  Biometrika 94, 19–35 (2007).
  * [24] M Timme, J Casadiego, Revealing networks from dynamics: an introduction.  Journal of Physics A: Mathematical and Theoretical 47, 343001 (2014).
  * [25] T Gaskin, GA Pavliotis, M Girolami, Neural parameter calibration for large-scale multi-agent models.  Proceedings of the National Academy of Sciences 120 (2023).
  * [26] National Grid, Transmission Network Shapefiles (https://www.nationalgrid.com/electricity-transmission/network-and-infrastructure/network-route-maps) (2023).
  * [27] SP Energy Networks, Transmission Network GIS Shapefiles (https://www.spenergynetworks.co.uk/pages/utility_map_viewer.aspx) (2023).
  * [28] Scottish and Southern Electricity Networks, Transmission Network GIS Shapefiles (https://www.ssen.co.uk/globalassets/library/connections---useful-documents/network-maps/5-gis-guide-shape-files-v1.pdf) (2023).
  * [29] Department for Business, Energy and Industrial Strategy, Digest of UK Energy Statistics 5: Electricity (2022).
  * [30] Y Kuramoto, Self-entrainment of a population of coupled non-linear oscillators in International Symposium on Mathematical Problems in Theoretical Physics, ed. H Araki.  (Springer-Verlag), pp. 420–422 (1975).
  * [31] G Filatrella, AH Nielsen, NF Pedersen, Analysis of a power grid using a Kuramoto-like model.  The European Physical Journal B 61, 485–491 (2008).
  * [32] M Rohden, A Sorge, M Timme, D Witthaut, Self-Organized Synchronization in Decentralized Power Grids.  Phys. Rev. Lett. 109, 064101 (2012).
  * [33] T Nishikawa, AE Motter, Comparative analysis of existing models for power-grid synchronization.  New Journal of Physics 17, 015012 (2015).
  * [34] YP Choi, Z Li, Synchronization of nonuniform Kuramoto oscillators for power grids with general connectivity and dampings.  Nonlinearity 32, 559–583 (2019).
  * [35] F Basiri, J Casadiego, M Timme, D Witthaut, Inferring power-grid topology in the face of uncertainties.  Phys. Rev. E 98, 012305 (2018).
  * [36] M Timme, Revealing Network Connectivity from Response Dynamics.  Physical Review Letters 98, 224101 (2007).
  * [37] Greater London Authority, 2017 Health Check Report (2017).
  * [38] Greater London Authority, 2015 ward profiles and atlas (2015) Online.
  * [39] A Wilson, A statistical theory of spatial distribution models.  Transportation Research 1, 253–269 (1967).
  * [40] Greater London Authority, Statistical boundary files for London (2011).
  * [41] B Harris, AG Wilson, Equilibrium Values and Dynamics of Attractiveness Terms in Production-Constrained Spatial-Interaction Models.  Environment and Planning A: Economy and Space 10, 371–388 (1978).
  * [42] K Csilléry, MG Blum, OE Gaggiotti, O François, Approximate Bayesian Computation (ABC) in practice.  Trends in Ecology & Evolution 25, 410–418 (2010).
  * [43] E Jennings, M Madigan, astroABC: An Approximate Bayesian Computation Sequential Monte Carlo sampler for cosmological parameter estimation.  Astronomy and Computing 19, 16–22 (2017).
  * [44] M Raissi, P Perdikaris, G Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics 378, 686–707 (2019).
  * [45] Y Gal, Z Ghahramani, Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning.  Proceedings of The 33rd International Conference on Machine Learning 48, 1050–1059 (2016).
  * [46] GA Pavliotis, A Zanoni, A method of moments estimator for interacting particle systems and their mean field limit.  arXiv 2212.00403 [math.NA] (2022).
  * [47] N Zagli, GA Pavliotis, V Lucarini, A Alecio, Dimension reduction of noisy interacting systems.  Phys. Rev. Res. 5, 013078 (2023).
  * [48] DP Kingma, J Ba, Adam: A Method for Stochastic Optimization.  arXiv 1412.6980 [cs.LG] (2014).
  * [49] AM Stuart, Inverse problems: A Bayesian perspective.  Acta Numerica 19, 451–559 (2010).
  * [50] L Riedel, B Herdeanu, H Mack, Y Sevinchan, J Weninger, Utopia: A Comprehensive and Collaborative Modeling Framework for Complex and Evolving Systems.  Journal of Open Source Software 5, 2165 (2020).
  * [51] Y Sevinchan, B Herdeanu, J Traub, dantro: a Python package for handling, transforming, and visualizing hierarchically structured data.  Journal of Open Source Software 5, 2316 (2020).
  * [52] National Grid ESO, What is Frequency? (https://www.nationalgrideso.com/electricity-explained/how-do-we-balance-grid/what-frequency) (2023).

Supporting Information

## Neural networks: Notation and Terminology

A _neural network_ is a sequence of length $L\geq 1$ of concatenated
transformations. Each _layer_ of the net consists of $L_{i}$ _neurons_ ,
connected through a sequence of _weight matrices_
$\mathbf{W}_{i}\leavevmode\nobreak\ \in\leavevmode\nobreak\
\mathbb{R}^{L_{i+1}\times L_{i}}$. Each layer applies the transformation

$\sigma_{i}(\mathbf{W}_{i}\mathbf{x}+\mathbf{b}_{i})$

to the input $\mathbf{x}$ from the previous layer, where
$\mathbf{b}_{i}\in\mathbb{R}^{L_{i+1}}$ is the _bias_ of the $i$-th layer. The
function $\sigma_{i}:\mathbb{R}^{L_{i+1}}\rightarrow\mathbb{R}^{L_{i+1}}$ is
the _activation function_ ; popular choices include the _rectified linear
unit_ (ReLU) activation function $\sigma(x)=\max(x,0)$, and the _sigmoid
activation function_ $\sigma(x)=(1+e^{-x})^{-1}$. A neural net has an _input
layer_ , an _output layer_ , and _hidden layers_ , which are the layers in
between the in- and output layers. If a network only has one hidden layer, we
call it _shallow_ , else we call the neural net _deep_.

[height=7, nodesize=25, nodespacing=12mm, layerspacing=28mm]

; ; ; ;

[count=5, bias=true, title=Input layer, text=$x_{\hiddenlayer}$[count=4,
bias=true, text=$h^{(1)}_{\linklayers}$[count=5, bias=true, title=Hidden
layers, text=$h^{(2)}_{\linklayers}$[count=6, bias=true,
text=$h^{(3)}_{\linklayers}$[count=3, title=Output layer,
text=$\hat{y}_{\linklayers}$

Figure S1: Example of a deep neural network with 3 hidden layers. The inputs
(light blue nodes) are passed through the layers, with links between layers
representing the weight matrices $\mathbf{W}$. Each layer also applies a
_bias_ (yellow nodes), with the network finally producing an output (orange).

(a) Distances $l_{ij}$

(b) Voltages $U_{ij}$

(c) Admittances $Y_{ij}$

(d) Scaled edge weights $a_{ij}\sim|Y_{ij}|U_{ij}^{2}$

Figure S2: Line data statistics for the British power grid. (a) Histogram of
the node distances in the network, with mean, standard deviation, and median
given. (b) Histogram of the line voltages. (c) Histogram of the line
admittances. As is visible, a small number of short lines artifically skew the
distribution. (d) The resulting normalised edge weights. Short edges with an
artificially high admittance are assigned a random value in $[0.9,1]$.

## Constructing the admittance matrix

Here we provide some additional information on the calculation of the network
edge weights from the data. The weights quantify the admittance (inverse
impedance) of the line, that is, how easily the line can transmit electrical
current. The characteristic impedance of a transmission line is given by

$Z_{0}=\sqrt{\dfrac{R+i\Omega L}{i\Omega C}},$ (23)

with $i$ the complex unit, $\Omega=50\ \mathrm{Hz}$ the grid frequency [52],
$L$ the cable inductance, and $C$ its capacitance. We use the following values
for a copper conductor with a cross-section of 1000 mm2, provided by a
standard manufacturer444caledonian-cables.co.uk/products/hv/400kv.shtml of
power grid cables: $R=0.0276$ $\Omega$/km (AC resistance), $L=0.41$ mH/km,
$C=0.150$ $\mu$F/km (assuming a single core). The total impedance along the
line is then given by

$Z(l)=Z_{0}\mathrm{sinh}(\gamma l),$ (24)

where $l$ is the length of the line, and the propagation constant $\gamma$ is
given by

$\gamma=\sqrt{(R+i\Omega L)(i\Omega C)}.$ (25)

In these calculations we neglect the conductance of the dielectric along the
line, which is neglibile at the distances present in the network. The
admittance of the line is then given by

$Y_{ij}=Z_{ij}^{-1},$ (26)

and the total edge weight $a_{ij}$ by

$a_{ij}=|Y_{ij}|U_{ij}^{2}.$ (27)

Long lines typically carry two to four times as many cables as shorter lines.
We account for this by multiplying the admittance of lines over 10 km by a
factor of 2, and those over 80 with a factor of 3 (admittances are additive).
Since the neural network outputs values in $[0,1]$, we scale the edge weights
to this range. However, a small number of short lines artificially skew the
distribution (see figure S2); these are mainly small segments designed to more
accurately capture the geometry of a longer line, though sometimes they
represent short lines to transformers, power stations, etc. To reduce their
impact on the weight distribution, instead of dividing the weights by the
maximum value, we divide by the mean and truncate all weights to $[0,1]$:

$a_{ij}\to\min(1,a_{ij}/\langle a_{ij}\rangle).$ (28)

The scaling factor is absorbed into the coupling coefficient $\kappa$. Those
edges with weight exactly equal to 1 are reassigned a value chosen uniformly
at random in $[0.9,1.0]$. As is visible in fig. S2d, the resulting
distribution of the weights $a_{ij}$ is more uniform than that of the
distances $l_{ij}$ and admittances $Y_{ij}$, since longer lines’ lower
admittance is often compensated by a higher voltage.

These calculations do not account for, among other things, the fact that we
are connecting the nodes with straight lines rather than the real line
trajectory, are not discerning between overground and underground lines, are
assuming a single core per cable, and are not considering differences in
transmission line heights, geometries, and materials. These inaccuracies are
absorbed into the coefficients $\alpha$ and $\beta$, which we tune manually in
order to allow for stable phase-locking.

## Details on the code

The code is uploaded to the Github repository as given in the main text. The
two models relevant to this work are Kuramoto and HarrisWilsonNW.

#### Installation

Detailed installation instructions are given in the repository. First, clone
the repository, install the utopya package and all the required additional
components into a virtual environment, for example via PyPi. In particular,
install pytorch. Enter the virtual environment. Then, from within the project
folder, register the project:

⬇

utopya projects register .

You should get a positive response from the utopya CLI and your project should
appear in the project list when calling:

⬇

utopya projects ls

Note that any changes to the project info file need to be communicated to
utopya by calling the registration command anew. You will then have to
additionally pass the \--exists-action overwrite flag, because a project of
that name already exists. See

⬇

utopya projects register \--help

for more information. Finally, register the project and its models via

⬇

utopya projects register . \--with-models

#### Running the code

To run a model, execute the following command:

⬇

utopya run <model_name>

By default, this runs the model with the settings in the <model_name>_cfg.yml
file. All data and the plots are written to an output directory, typically
located in ~/utopya_output. To run the model with different settings, create a
run_cfg.yml file and pass it to the model like this:

⬇

utopya run <model_name> path/to/run_cfg.yml

This is recommended rather than changing the default settings, because the
defaults are parameters that are known to work and you may wish to fall back
on in the future.

Plots are generated using the plots specified in the <model_name>_plots.yml
file. These too can be updated by creating a custom plot configuration, and
running the model like this:

⬇

utopya run <model_name> path/to/run_cfg.yml \--plots-cfg path/to/plot_cfg.yml

See the Utopia tutorial for more detailed instructions.

All the images in this article can be generated using so-called _configuration
sets_ , which are complete bundles of both run configurations and evaluation
configurations. For example, to generate the predictions on the random network
with $N=1000$ nodes (fig. 4d–e) for the Kuramoto model, you can call

⬇

utopya run Kuramoto \--cfg-set N_1000_example

This will run and evaluate the Kuramoto model with all the settings from the
Kuramoto/cfgs/N_100_example/run.yml and eval.yml configurations.

#### Parameter sweeps

Parameter sweeps are automatically parallelised by utopya, meaning simulation
runs are always run concurrently whenever possible. The data is automatically
stored and loaded into a data tree. To run a sweep, simply add a !sweep tag to
the parameters you wish to sweep over, and specify the values, along with a
default value to be used if no sweep is performed:

⬇

param: !sweep

default: 0

values: [0, 1, 2, 3]

Then in the run configuration, add the following entry:

⬇

perform_sweep: true

Alternatively, call the model with the flag \--run-mode sweep. The
configuration sets used in this work automatically run sweeps whenever needed,
so no adjustment is needed to recreate the plots used in this work.

#### Initialising the neural net

The neural net is controlled from the NeuralNet entry of the configuration:

⬇

NeuralNet:

num_layers: 4

nodes_per_layer:

default: 20

layer_specific:

1: 10

2: 15

biases: # optional; if this entry is omitted no biases are used

default: ~ # default is None (indicted by a tilde in YAML)

layer_specific:

0: default # use pytorch default (xavier uniform)

-1: [-1, 1] # uniform initialisation on a custom interval

activation_funcs:

default: sigmoid

layer_specific: # optional

1: tanh

3:

name: HardTanh # you can also pass a function that takes additional args
and/or kwargs

args:

- -2 # min_value

- +2 # max_value

learning_rate: 0.001 # optional; default is 0.001

optimizer: SGD # optional; default is Adam

num_layers specifies the depth of the net; nodes_per_layer controls the
architecture: provide a default size, and optionally any deviations from the
default under layer_specific. The keys of the layer_specific entry should be
indices of the layer in question. The optional biases entry determines whether
or not biases are to be included in the architecture, and if so how to
initialise them. A default and layer-specific values can again be passed.
Setting an entry to default initialises the values using the pytorch default
initialiser, a Xavier uniform initialisation. Passing a custom interval
instead initialises the biases uniformly at random on that interval, and
passing a tilde ~ (None in YAML) turns the bias off. activation_funcs is a
dictionary specifying the activation functions on each layer, following the
same logic as above. Any pytorch activation function is permissible. If a
function requires additional arguments, these can be passed as in the example
above. Lastly, the optimizer keyword takes any argument allowed in pytorch.
The default optimizer is the Adam optimizer [48] with a learning rate of
$0.02$.

The neural net can be initialised from different initial values in the
parameter space by changing the _random seed_ in the configuration:

⬇

seed: 42

Sweeping over different initialisations is achieved by sweeping over the seed,
as described in the previous section.

#### Training the neural net

The Training entry of the configuration controls the training process:

⬇

Training:

batch_size: 2

device: cpu

true_parameters:

sigma: 0.0

loss_function:

name: MSEloss

# can pass additional args and kwargs here ...

You must specify the noise level to use for the ABM during training; the
default value is $0$. Under the loss_function key you can specify the loss
function to use, and pass any arguments or keyword arguments it may require
using an args or kwargs key. You can use any available pytorch loss function.

The device key sets the training device. The default is the CPU, but you can
also train on the GPU by setting the device to cuda. Note that on Apple
Silicon, the device name is mps. Make sure you have installed the correct
package for your device, and note that, as of writing, some pytorch functions
required for our code (e.g., the trace and hard sigmoid functions) had not yet
been implemented for MPS, hence GPU training on Apple Silicon devices was not
possible.

## Kuramoto model

### Neural Network Architecture

The following is the default configuration for the neural network and training
settings used for the Kuramoto model:

⬇

Kuramoto:

NeuralNet:

num_layers: 5

nodes_per_layer:

default: 20

biases:

default: ~ # No biases

activation_funcs:

default: tanh

layer_specific:

-1: HardSigmoid # hard sigmoid on the last layer

learning_rate: 0.002

optimizer: Adam

Training:

batch_size: 2

loss_function:

name: MSELoss

true_parameters:

sigma: 0

### Training

We train the model using the following operation: let
$\bm{\varphi}(t)=(\varphi_{1}(t),...\varphi_{N}(t))$ be the current phases of
the $N$ nodes,
$\bm{\dot{\varphi}}(t)=(\varphi_{1}(t)-\varphi_{1}(t-1),...,\varphi_{N}(t)-\varphi_{N}(t-1)$
the vector of phase derivatives, $\bm{\omega}=(\omega_{1},...,\omega_{N})$ the
vector of eigenfrequencies, and
$\bm{\Gamma}(t)=(\sin(\varphi_{j}(t)-\varphi_{i}(t))_{ij}$; then in each
iteration of the first-order Kuramoto model ($\beta\neq 0$), we do

$\bm{\varphi}(t+1)=\bm{\varphi}(t)+\dfrac{1}{\beta}\left(\bm{\omega}(t)+\mathrm{diag}\left(\mathbf{\hat{A}}\bm{\Gamma}(t)\right)\right)\mathrm{d}t,$
(29)

where $\mathrm{diag}(\cdot)$ takes the diagonal elements of the matrix. For
the second-order model ($\alpha\neq 0$), we do

$\bm{\varphi}(t+1)=\bm{\varphi}(t)+\dfrac{1}{\alpha}\left[\left(\bm{\omega}(t)+\mathrm{diag}\left(\mathbf{\hat{A}}\bm{\Gamma}(t)\right)-\beta\bm{\dot{\varphi}}(t)\right)\mathrm{d}t+\bm{\dot{\varphi}}(t)\right]\mathrm{d}t.$
(30)

In the first-order case, the initial conditions required are
$\bm{\varphi}(0)$, in the second order case, we need $\bm{\varphi}(0)$ and
$\bm{\varphi}(1)$ (i.e. the initial phases and initial velocities).

### Running the code

Configuration sets to reproduce all numerical experiments except the British
power grid inference are provided; for instance, to produce the predictions
for a random network with $N=100$ nodes, simply run

⬇

utopya run Kuramoto \--cfg-set N_1000_example

## Harris-Wilson model

### Neural Network Architecture

The following is the default configuration for the neural network and training
settings used for the Harris-Wilson model:

⬇

HarrisWilsonNW:

NeuralNet:

num_layers: 2

nodes_per_layer:

default: 20

biases:

default: ~ # No biases

activation_funcs:

default: tanh

layer_specific:

-1: sigmoid # sigmoid on the last layer

learning_rate: 0.002

optimizer: Adam

Training:

batch_size: 2

loss_function:

name: MSELoss

true_parameters:

sigma: 0

### Training

We use the following matrix form of the Harris-Wilson equations to train the
neural net. Let $\mathbf{D}\in\mathbb{R}^{M}$, $\mathbf{O}\in\mathbb{R}^{N}$,
$\mathbf{W}\in\mathbb{R}^{M}$ be the demand vector, origin zone size vector,
and destination zone size vector respectively. The the dynamics are given by

$\mathbf{D}=\mathbf{W}^{\alpha}\odot\left[\left(\mathbf{C}^{\beta}\right)^{T}(\mathbf{O}\odot\mathbf{Z})\right]\in\mathbb{R}^{M},$
(31)

with $\odot$ indicating the Hadamard product, and elementwise exponentiation.
$\mathbf{Z}\in\mathbb{R}^{N}$ is the vector of normalisation constants

$\mathbf{Z}^{-1}=\mathbf{C}^{\beta}\mathbf{W}^{\alpha}.$ (32)

The dynamics then read

$\dot{\mathbf{W}}=\mathbf{W}\odot\epsilon(\mathbf{D}-\kappa\mathbf{W})$ (33)

with given initial conditions $\mathbf{W}(t=0)=\mathbf{W}_{0}$. This
formulation is more conducive to machine learning purposes, since it contains
easily differentiable matrix operations and does not use for-loop iteration.

### Running the code

This is analogous to the Kuramoto case. To reproduce the plot from the main
article, run the following command:

⬇

utopya run HarrisWilsonNW \--cfg-set London_dataset