# Uncovering the non-equilibrium stationary properties in sparse Boolean
networks

Giuseppe Torrisi, Reimer Kühn, Alessia Annibale

###### Abstract

Dynamic processes of interacting units on a network are out of equilibrium in
general. In the case of a directed tree, the dynamic cavity method provides an
efficient tool that characterises the dynamic trajectory of the process for
the linear threshold model. However, because of the computational complexity
of the method, the analysis has been limited to systems where the largest
number of neighbours is small. We devise an efficient implementation of the
dynamic cavity method which substantially reduces the computational complexity
of the method for systems with discrete couplings. Our approach opens up the
possibility to investigate the dynamic properties of networks with fat-tailed
degree distribution. We exploit this new implementation to study properties of
the non-equilibrium steady-state. We extend the dynamical cavity approach to
calculate the pairwise correlations induced by different motifs in the
network. Our results suggest that just two basic motifs of the network are
able to accurately describe the entire statistics of observed correlations.
Finally, we investigate models defined on networks containing bi-directional
interactions. We observe that the stationary state associated with networks
with symmetric or anti-symmetric interactions is biased towards the active or
inactive state respectively, even if independent interaction entries are drawn
from a symmetric distribution. This phenomenon, which can be regarded as a
form of spontaneous symmetry-breaking, is peculiar to systems formulated in
terms of Boolean variables, as opposed to Ising spins.

##  1 Introduction

A broad range of disordered and complex systems can be modelled in terms of
binary units, which update their state according to stochastic Boolean
functions of the neighbouring units. They include spin glasses [1, 2, 3] in
physics, gene-regulatory and immune networks [4, 5, 6, 7, 8] in biology,
artificial neural networks [9, 10, 11, 12] in computer science, agent-based
models [13, 14, 15, 16], models of operational or credit risk [17, 18] in
economics and finance, and a variety of hard combinatorial optimization
problems [19, 20, 21, 22, 23] to name but a few. In recent years there was, in
particular, a surge of interest in systems where each unit interacts only with
a finite set of neighbours, and which are thus defined on finitely connected
networks.

Statistical mechanics provides a rich arsenal of tools to analyse the
equilibrium behaviour of sparse systems [24, 3, 25, 26]. However, the
assumption of equilibrium is not satisfied in systems that are driven,
dissipative or exhibit a degree of asymmetry in the interactions. Here, one
has to resort to tools of non-equilibrium statistical mechanics, which are
however much more challenging. Non-equilibrium dynamics has been analysed in
explicit detail for a class of models where the dynamics of units has an
absorbing state [27, 28, 29, 30] often combined with the simplifying features
that interactions are of pure product form, e.g. in contagion processes [31,
32].

However, the dynamics of a large class of Boolean models does not have an
absorbing state. In the present manuscript, we analyse the stochastic dynamics
of linear threshold models, formulated in terms of $\\{0,1\\}$ variables, and
discuss extensions to threshold models with multi-node interactions and models
with Ising spin $\\{\pm 1\\}$, which do not have absorbing states. In such
models, approximation schemes that are successfully used for dense systems,
such as the heterogeneous mean-field and TAP approaches [23, 33], have been
shown to be ineffective for sparse systems [34, 35]. On the other hand,
generating functional analyses [36] can accurately characterise site averaged
quantities, such as global magnetization and time-lagged correlations [37,
38], in sparse heterogeneous systems. However they require averaging over
graph ensembles, hence they are not able to describe single instances. These
can instead be investigated by the dynamic cavity method [38, 39, 29, 32]
which is particularly effective in the analysis of sparse systems, whenever
short loops are rare. In particular, the dynamic cavity method can potentially
investigate dynamic properties at the level of individual nodes [39]. Single-
node statistics has, for instance, been shown to be highly heterogeneous for
non-equilibrium models with an absorbing state [32, 40, 27]. Little, however,
is known for other models, including spin models and linear threshold models.

In particular, explicit analysis has been feasible only for quite a narrow
choice of networks, namely directed trees and fully asymmetric random graphs,
where bi-directional links are sampled independently [38, 39]. Moreover, both
the dynamic cavity and generating functional methods require evaluating
averages over a configuration space that grows exponentially with the in-
degree of the nodes, which limits the applicability of both methods to
networks with small in-degrees. Due to this computational complexity barrier,
a large class of networks, in particular those characterised by a fat-tailed
degree distribution, are out of reach for analytical characterisation.
Networks with fat-tailed degree distribution are, however, known to appear in
many real-world problems and represent the outcome of generating network
models based on preferential attachment [41, 42, 43, 44, 45, 46]. They are
characterised by large differences in the local environment surrounding nodes,
which are known to lead to heterogeneous node properties at equilibrium [26].
Therefore, the aforementioned exponential complexity in the in-degrees has so
far prevented a detailed study of node heterogeneities, precisely in those
systems where they play a bigger role.

We have recently proposed a method based on dynamic programming that overcomes
this complexity barrier when couplings are chosen randomly from any discrete
set of equidistant values, the simplest example being couplings chosen from
the set $\\{\pm J\\}$, see Ref. [47]. Our algorithm reduces the computational
complexity in the in-degrees of the system from exponential to quadratic, thus
allowing us to characterise the probabilistic evolution of individual node
states in networks with a broad in-degree distribution.

In this manuscript, we extend our analysis of heterogeneous properties of the
linear threshold model in fully asymmetric networks to include pairwise
correlations between node states. Our analysis shows that the statistics of
pairwise correlation is almost completely captured in terms of simple motifs
of length two. This reveals the extent to which, in terms of path distance,
node states significantly influence each other.

In addition, we extend our formalism to investigate systems with multi-node
interactions, using bipartite (or factor) graphs, which provide a more
flexible model for combinatorial control; see Ref. [48]. Multi-node
interactions models can represent generic Boolean functions, which are of
particular interest in the context of gene regulation [4, 49, 50, 51]. We
present our extended formalism for models of indicator variables $0,1$ as well
as models of Ising spins $\pm 1$, for which our formulation in terms of
bipartite graphs leads to (sparse) mixed $p$-spin models.

Finally, we extend our analysis to systems with bi-directional links. For this
class of systems, the cavity method requires following the probability of the
full dynamic trajectory, which cannot be reduced to the product of one-time-
step terms, because bi-directional links introduce retarded self-interactions
of nodes with their past. This leads to an exponential scaling with the time
horizon considered, which effectively restricts the application of the method
to follow the dynamics over a few time-steps only, thus making the
investigation of long time behaviour infeasible. An approximation technique
called the one-time approximation (OTA) [39, 35, 34] has been proposed to
overcome these restrictions.

The OTA factorises the probabilities of dynamic trajectories into one-time-
step probabilities, and it leads to two sets of coupled equations, one for the
marginals and one for the cavity marginals, which become closed under further
assumptions. This makes the problem solvable also for the stationary state.
Both equations are affected by the aforementioned in-degree complexity, hence
the OTA greatly benefits from our dynamic programming method. Within the OTA
framework, different closure schemes have been proposed in the literature [35,
34]. We evaluate the impact that different closure schemes have on the
results, and we use such analysis to motivate our choice. A thorough
inspection of the result leads to a surprising finding: networks with the same
number of positive and negative interactions may sustain a biased distribution
of node activation in presence of bi-directional links. This phenomenon, which
can be regarded as a form of spontaneous symmetry breaking, is peculiar of
systems with $\\{0,1\\}$ representation of states (as opposed to Ising spins
with $\\{\pm 1\\}$ representation).

Our paper is organized as follows. In Sec. 2 we describe the dynamic
programming approach to the linear threshold model on fully asymmetric
networks with fat-tailed degree distribution. Its non-equilibrium stationary
state is characterised in terms of individual nodes statistics in Sec. 3, and
pairwise site correlations in Sec. 4. In Sec. 5 we extend our method both to
models with multi-node interactions, formulated in terms of bipartite
networks, and to networks of Ising spins. In Sec. 6 we include bi-directional
links in the analysis of linear threshold models, using the one-time
approximation. In Sec. 7 we investigate the unexpected occurrence of a form of
spontaneous symmetry breaking in these systems. We summarise and discuss our
results in Sec. 8. Finally, we include two appendices: the first containing
the high noise limit for the asymmetric network model and our multi-node
interactions model, the second containing details concerning the derivation of
OTA schemes. The interested reader can find the code to reproduce the results
shown in this paper at the following link. 111
https://doi.org/10.5281/zenodo.5996772

##  2 Derivation of dynamic cavity equation on fully asymmetric networks

We consider a directed network consisting of $N$ nodes labelled $i=1,\dots N$.
For each edge $(i,j)$, the edge weight $J_{ij}\in\mathbb{R}$ defines the
strength of the interaction carried from node $j$ to node $i$, while
$J_{ij}=0$ if the link $(i,j)$ does not exist. In this section we will
consider fully asymmetric networks, i.e. directed networks such that $\forall
i,j$ if link $(i,j)$ exists, the link in the opposite direction $(j,i)$ does
not exist. We denote $\partial_{i}=\\{j:J_{ij}\neq 0\\}$ the set of
predecessors of node $i$ and $k_{i}^{\mathrm{in}}=|\partial_{i}|$ the in-
degree of node $i$. We will drop the superscript “$\mathrm{in}$” in
$k_{i}^{\mathrm{in}}$ whenever it does not lead to ambiguity. Every node $i$
can be in one of two states, described by a binary state variable
$n_{i}\in\\{0,1\\}$ with $i\in\\{1,\dots N\\}$. Several Boolean models can be
defined on a given graph depending on the choice of the update rule of node
states, a prominent example being the random Boolean network model [4].

In the present section, we consider the linear threshold model defined on a
weighted directed network, which is specified by the interaction matrix
$\bm{J}=(J_{ij})$ and thresholds $\\{\vartheta_{i}\\}_{i=1}^{N}$. A node $i$
evaluates a local field

$h_{i}(\bm{n}_{\partial_{i}})\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\sum_{j}J_{ij}n_{j}\,,$
(1)

and becomes active in the next time-step if the sum of the local field and a
random noise contribution $z_{i}$ exceeds the threshold $\vartheta_{i}$, i.e.,

$n_{i}(t+1)=\Theta\left(h_{i}\left(\bm{n}(t)\right)-\vartheta_{i}-z_{i}(t)\right),$
(2)

where the $z_{i}(t)$ are independent identically distributed random variables,
with cumulative distribution function
$\mathrm{Prob}[z_{i}(t)<x]\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\Phi_{\beta}(x)$,
for all $i,t$, and $T=\beta^{-1}$ is a parameter that characterises the noise
strength. The local field $h_{i}(\bm{n}_{\partial_{i}})$ depends on the states
of nodes that are predecessors of node $i$. We use
$\bm{n}_{\partial_{i}}=\\{n_{j},j\in{\partial_{i}}\\}$ to denote the states of
this set. In the following, we explicitly write the argument of the local
field using the notation $h_{i}(\bm{n}_{\partial_{i}})$. Given the states
$\bm{n}_{\partial_{i}}$, the probability that node $i$ is active at the next
time-step is $\Phi_{\beta}(h_{i}(\bm{n}_{\partial_{i}})-\vartheta_{i})$, the
transition probability to state $n_{i}$, for node $i$, is

$W[n_{i}|h_{i}(\bm{n}_{\partial_{i}}),\vartheta_{i}]=n_{i}\Phi_{\beta}(h_{i}(\bm{n}_{\partial_{i}})-\vartheta_{i})+(1-n_{i})(1-\Phi_{\beta}(h_{i}(\bm{n}_{\partial_{i}})-\vartheta_{i}))\,.$
(3)

Given the joint probability $P(\bm{n}_{\partial_{i}},t)$ of the states
$\bm{n}_{\partial_{i}}$ at time t, the probability that node $i$ is in state
$n_{i}$ at time $t+1$ is given by

$P(n_{i},t+1)=\sum_{\bm{n}_{\partial_{i}}}W[n_{i}|h_{i}(\bm{n}_{\partial_{i}}),\vartheta_{i}]P(\bm{n}_{\partial_{i}},t)\,,$
(4)

where
$\sum_{\bm{n}_{\partial_{i}}}(\cdot)\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\prod_{j\in\partial_{i}}\sum_{n_{j}=\\{0,1\\}}(\cdot)$
indicates the sum over all possible configurations of the microscopic
variables. For finitely coordinated random graphs of the type considered here,
the cavity method [24] can be used to analyse the dynamics of the system. It
is exact on trees and known to become exact for finitely coordinated random
graphs in the thermodynamic limit, as the typical length of any loops in such
systems diverges logarithmically in system size $N$. On the cavity graph where
node $i$, along with the edges connected to it, has been removed, the cavity
approximation assumes

$P^{(i)}(\bm{n}_{\partial_{i}},t)=\prod_{j\in\partial_{i}}P^{(i)}(n_{j},t)\,,$
(5)

which holds on a tree-like structure. The probability
$P^{(i)}(\bm{n}_{\partial_{i}},t)$ of the states of nodes neighbouring $i$
will in general be different from the corresponding probability on the
original graph, i.e. $P(\bm{n}_{\partial_{i}},t)\neq
P^{(i)}(\bm{n}_{\partial_{i}},t)$, but equality does hold on a directed tree
and it holds approximately on a random fully asymmetric graph, in which loops
are rare [38, 39, 50]222Cavity method can be applied also in the general case,
as discussed in [39] and in B. However the evaluation is more challenging..
Using the equivalence of the cavity and non-cavity dynamics, from Equation (4)

$P(n_{i},t+1)=\sum_{\bm{n}_{\partial_{i}}}W[n_{i}|h_{i}(\bm{n}_{\partial_{i}}),\vartheta_{i}]\prod_{j\in\partial_{i}}P(n_{j},t)\,.$
(6)

In particular, the probability of node activation becomes

$P(n_{i}=1,t+1)=\sum_{\bm{n}_{\partial_{i}}}\Phi_{\beta}\left(h_{i}(\bm{n}_{\partial_{i}})-\vartheta_{i}\right)\prod_{j\in\partial_{i}}P(n_{j},t)\,.$
(7)

Let us denote by
$P_{i}(t)\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=P(n_{i}=1,t)$
the activation probability of node $i$ at time $t$, let the angle bracket
$\langle\,\cdot\,\rangle_{\bm{n}_{\partial_{i}},\,t}$ indicates an average
over the states of the predecessors of node $i$, evaluated using their joint
node activation probabilities at time $t$

$\langle(~{}\cdot~{})\rangle_{\bm{n}_{\partial_{i}},\,t}\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\sum_{\bm{n}_{\partial_{i}}}\prod_{j\in\partial_{i}}P_{j}(t)^{n_{j}}[1-P_{j}(t)]^{1-n_{j}}\,(~{}\cdot~{})\,.$
(8)

Combining Equation (8) and Equation (7) a concise expression for the
activation probability $P_{i}(t+1)$ is obtained as

$P_{i}(t+1)=\left\langle\Phi_{\beta}\left(h_{i}(\bm{n}_{\partial_{i}})-\vartheta_{i}\right)\right\rangle_{\bm{n}_{\partial_{i}},\,t}\
.$ (9)

The average in Equation (9) is evaluated over a configuration space containing
$2^{k_{i}}$ elements, as shown in Equation (8), thus the evaluation of the
average hits an (exponential) complexity barrier whenever $k_{i}\gg 1$. In
Ref. [47] we have proposed a new method that uses dynamic programming to
efficiently perform the average in Equation (9) using a polynomial number of
operations. In the following, we recap the main ideas of the method.

#### Dynamic programming

Equation (9) depends on the distribution of states $\mathbf{n}_{\partial_{i}}$
through the local field $h_{i}(\mathbf{n}_{\partial_{i}})$. If couplings are
extracted from a discrete set of values, e.g., $J_{ij}=\\{0,\pm J\\}$, there
is a degeneracy of values of $h_{i}(\bm{n}_{\partial_{i}})$, i.e., several
configurations of $k_{i}$ states $(n_{j_{1}},\dots n_{j_{k_{i}}})$ are
associated with the same value of local field. In fact, for a given
realisation of the interaction terms $\left(J_{ij_{1}}\dots
J_{ij_{k_{i}}}\right)$ the number of possible values of
$h_{i}(\bm{n}_{\partial_{i}})$ is $k_{i}+1$, while the number of
configurations of states is $2^{k_{i}}$. We exploit this degeneracy by
averaging directly over the local fields, rather than over the micro-states of
the set of predecessors, using a dynamic programming approach. We show below
our approach reduces the computational complexity of the problem from
$\mathcal{O}(2^{k_{i}})$ to $\mathcal{O}(k_{i}^{2})$, [47].

Let us consider a node $i$ with in-degree $k_{i}$, and let us use $\\{1,\dots
k_{i}\\}$ to label the indices of the predecessors of node $i$. We first
define the sub-problem

$f_{i}(\ell,\tilde{h})=\left\langle\\!\Phi_{\beta}\Bigg{(}\tilde{h}+\sum_{j=\ell}^{k_{i}}J_{ij}n_{j}-\vartheta_{i}\Bigg{)}\right\rangle_{n_{\ell,\dots,k_{i}},t}\,\mbox{for
}\ell\in\\{1,\dots k_{i}\\},$ (10)

which consists in performing the average over a sub-set of nodes. We will
refer to $\tilde{h}$ as an auxiliary field. At any given $\ell$, the set of
$f_{i}(\ell,\tilde{h})$ of interest would in fact correspond to the set of
averages representing the sub-problems of Equation (9) that remain to be
evaluated after the average over the first $\ell-1$ nodes has been performed.
Quantities in Equation (10) are not performed directly, but through the
recursive relationship

$f_{i}(\ell,\tilde{h})=P_{\ell}(t)\,f_{i}(\ell+1,\tilde{h}+J_{i\ell})+\big{(}1-P_{\ell}(t)\big{)}\,f_{i}(\ell+1,\tilde{h})$
(11)

for $1\leq\ell\leq k_{i}$, with the terminal boundary condition

$f_{i}(k_{i}+1,\tilde{h})=\Phi_{\beta}\left(\tilde{h}-\vartheta_{i}\right)\ .$
(12)

The original average $P_{i}(t+1)$ of Equation (9) that we are ultimately
interested in is, within this backward iteration scheme, obtained as

$P_{i}(t+1)=f_{i}(1,0)\ .$ (13)

In our earlier work, Ref. [47], we have shown that Equation (11) leads to a
dramatic reduction of the computational complexity for binary interactions
$J_{ij}\in\\{0,\pm J\\}$.The computational and memory requirement for the
evaluation of $P_{i}(t+1)$ is seen to scale as $k_{i}^{2}$, for a node with
in-degree $k_{i}$, rather than $2^{k_{i}}$, as it would in a naive evaluation.
This entails that the complexity of our algorithm is
$\mathcal{O}{(\sum_{i}k_{i}^{2})}$.

##  3 Node activation heterogeneity

It has been shown that fully asymmetric graphs do not exhibit a spin-glass
phase in a variety of dynamical models, and the time-dependent local
magnetizations or node activation probabilities are expected to converge to a
fixed point [52, 53, 37, 39]. We exploit this fact to investigate the
trajectory of nodes’ states by iteration of Equation (9). Our dynamic
programming implementation provides an efficient method to evaluate node
activation probabilities. We apply our formalism to inspect the dynamical
properties of networks with fat-tailed degree distributions. The theoretical
evaluations presented in this section would not be feasible without our
dynamic programming implementation.

#### Model Parameters Used in this Study

In the present section, we consider models of synthetic networks in the
configuration model class, i.e. we will study ensembles of maximally random
directed graphs ${\bf J}$ with constrained in-degree and out-degree sequences
of nodes. The in-degree of nodes are drawn from the distribution
$\rho^{\mathrm{in}}(k)=N^{-1}\sum_{i}\delta_{k,k_{i}^{\mathrm{in}}}$. In this
paper, we use the same distribution for the out-degree of nodes
$k_{i}^{\mathrm{out}}$.333If the distributions of -in- and out-degree are
drown from different distributions, a minimal requirement is
$\sum_{i}k_{i}^{\mathrm{in}}=\sum_{i}k_{i}^{\mathrm{out}}$ that guarantees the
conservation of the number of links. The in-degree and out-degree for the same
node are drawn independently. We use $\langle k^{\mathrm{in}}\rangle$ to
denote the mean in-degree. Our synthetic model uses parameters observed in
real-world networks. Specifically, we use parameters extracted from a gene
regulatory network [54], which is characterised by a node in-degree
distribution with a power-law tail $P_{\gamma}(k)\sim\gamma k^{-\gamma-1}$ and
$\gamma=2.81$, which we implement in terms of the discrete fat-tailed
distribution $P_{\gamma}(k)=k^{-\gamma}-(k+1)^{-\gamma}$; it is defined for
positive integers $1\leq k\in\mathbb{N}$ and has the desired power-law
behaviour for large $k$. For every link $(i,j)$, the corresponding interaction
term $J_{ij}$ is sampled from the distribution
$P(J_{ij})=\eta\delta(J_{ij}-J)+(1-\eta)\delta(J_{ij}+J)$ and we use
$\eta=0.621$ as suggested by the data in [54] and $J=1/\sqrt{\langle
k^{\mathrm{in}}\rangle}\approx 0.72$ unless stated otherwise. We set the value
of the threshold to be the same for all sites, $\vartheta_{i}=\vartheta$
$\forall i$. We assume a “thermal” noise distribution with
$\Phi_{\beta}(x)=1/2\left(1+\tanh\frac{\beta x}{2}\right)$.

###  3.1 Comparison between theory and simulation

The dynamic programming algorithm allows to speed up the evaluation of the
node activation probability at time $t+1$ given the knowledge of node
activation probabilities at time $t$. By iterating the procedure over the
time-steps of interest, we solve the dynamics for the nodes activation
probabilities. In this section, we compare the analytical results against
Monte Carlo simulations of the microscopic dynamics. Several results in the
literature indicate that the dynamic cavity method is effective to reproduce
the trajectory of fully asymmetric graphs, see Refs. [39, 55]. While those
results inspect the site average of node activations for Ising spins, we
consider the full distribution of node activations for indicator variables. In
particular, we compare the stationary values of the node activation
probabilities $P_{i}=\lim_{t\rightarrow\infty}P_{i}(t)$ $\forall
i\in\\{1,\dots N\\}$ obtained from simulations with those obtained via the
cavity approach. The stationary activation probabilities are obtained by
running the iterative procedure Equation (9) until convergence. The
convergence is controlled by measuring the difference between activation
probabilities $\epsilon_{t+1}=\mathrm{max}_{i}|P_{i}(t+1)-P_{i}(t)|$ in two
consecutive steps; this criterion is adopted in the entire manuscript. Let us
call $P_{i}=P_{i}(t)$ corresponding to the smallest $t$ satisfying the exit
condition, which we took to be $\epsilon_{t}<10^{-4}$. We test cavity
predictions against simulations of the microscopic dynamics Equation (2). We
evaluate the stationary node activation probability of simulated dynamical
trajectories by taking the sample average of node activations

$P_{i,MC}=\frac{1}{t_{s}}\sum_{t={t_{0}}}^{t_{0}+t_{s}}n_{i}(t)\,,$ (14)

for $t_{s}\gg 1$, and $t_{0}$ a sufficiently long time to allow the dynamics
to relax to stationarity. In the following, we use the general notation
$\Pi(x)$ to denote the probability density function for the realisation of
local quantities $x_{i}$, e.g. $\Pi(P)=N^{-1}\sum_{i}\delta(P_{i}-P)$. In Fig.
1 we show the distribution of the node activation probability as computed via
the cavity method and by estimation from simulations. These are found to be in
excellent agreement, confirming the validity of the cavity method to describe
the stationary state, see Fig. 1. We note that, to reach the resolution
imposed by cavity $\epsilon_{t}<10^{-4}$ through simulations, we need to
simulate long trajectories in Equation (14) with $t_{s}\sim 10^{8}$, which
makes the evaluation from simulation computationally far more expensive than
an evaluation based on the cavity method. For the same realisation of a
network of $N=100,000$ nodes and identical parameter choices, the cavity
implementation took $9\text{\,}\mathrm{s}$, whereas simulations took more than
$11\text{\,}\mathrm{d}$ays on the same machine.

Figure 1: Left: distribution of the node activation probability, as computed
from cavity theoretical calculation (blue solid line) and from Monte Carlo
direct simulation of the microscopic dynamics (dotted green line). Right:
distribution of the differences $\delta_{i}=P_{i,MC}-P_{i}$, evaluated for
each node $i$ of the network. Parameters:
$N=100,000,T/J=0.2,\gamma=2.81,\eta=0.5,\vartheta/J=0.1,J=1$.

#### Time evolution

We investigate the statistical properties of dynamical trajectories through
forward iteration of Equation (9). We initialise the system with random
initial conditions for the node activation probabilities $P_{i}(0)$. We stop
iterations at convergence using the same criterion as above. The time
evolution of the dynamics for each site is shown in Fig. 2. Node activation
probabilities converge after very few time-steps. This property is linked to
the absence of feedback signals in the dynamics, thanks to the directed tree-
like structure of the network. The probability density function of the node
activation probabilities at stationarity is shown in Fig. 3 for different
noise levels $T$. The distribution of the node activation probability is not
uni-modal, indicating the need to inspect the full distribution of node
activation, as measures of central tendency are not representative. This is
also true in the high noise limit, where the problem is directly solvable; see
A.

The multi-modal structure of the node activation distribution can be
rationalised by simple reasoning. Thanks to the fat tailed property of the
network, more than $\mathrm{8}\mathrm{7}\mathrm{\char 37\relax}$ of nodes have
in-degree 1. Let us consider a node $i$ with only one predecessor $j$, let us
call $\rho_{\pm}(P)$ the probability that node $i$ is active, given that
$J_{ij}=\pm J$, and that the node activation probability of that predecessor
is $P$, i.e., $\rho_{\pm}(P)=\left\langle\Phi_{\beta}(\pm
Jn-\vartheta)\right\rangle_{n}$ with $\langle n\rangle_{n}=P$. The explicit
expression of functions $\rho_{\pm}(P)$ becomes

$\begin{split}\rho_{\pm}(P)&=\frac{P}{2}\left[1+\tanh\frac{\beta}{2}\left(\pm
J-\vartheta\right)\right]+\frac{1-P}{2}\left[1-\tanh\frac{\beta\vartheta}{2}\right]\\\
&=\frac{1}{2}\left[1-\tanh\frac{\beta\vartheta}{2}\right]+\frac{P}{2}\left[\tanh\frac{\beta}{2}\left(\pm
J-\vartheta\right)+\tanh\frac{\beta\vartheta}{2}\right]\,.\end{split}$ (15)

The properties of nodes with in-degree 1 can be captured in terms of the map
shown in Fig. 4. We denote $p^{*}$ the fixed point defined as
$p^{*}=\rho_{+}(p^{*})$. The composition of function $\rho_{\pm}$ is able to
describe the trends observed in Fig. 3. In Fig. 5, we highlight the
contribution to the distribution of node activation probability from nodes
with in-degree 1. The values of $\rho_{\pm}(1)$ give the lower and upper
bounds of the node activation probability for nodes with in-degree 1. Even
though nodes with in-degree 1 account for the vast majority of nodes, the
dynamical properties of node activation depends on the contribution of other
nodes as well.

Figure 2: Trajectory of probabilities of node activation up to convergence.
The colour of the pixel located at row $t$ and column $i$ maps the value
$P_{i}(t)$. Colour-bar on the right maps the colours to the value of
activation probability. For visualisation purposes nodes are sorted by the
value of $P$ at the end of iteration, i.e., column locations are sorted
according to their value in the top row. Parameters: $N=200000$,
$\vartheta/J=0.2$, $T/J=0.3,J\approx 0.72$, $\eta=0.621$. Figure 3: Heat-map
of the distribution of node activation probabilities $\Pi(P)$ at different
values of noise parameter $T/J$ (top). Dashed, dotted, and dot-dashed lines
show the graph of different function compositions of $\rho_{\pm}$ at different
values of $T$, see Equation (15) . Vertical dotted lines mark the noise levels
at which the histograms of $\Pi(P)$ are shown in the lower panels. They probe
the regimes at low, intermediate and high noise, corresponding to
$T/J=0.1,0.3,0.6$ from left to right. The central lower panel corresponds to
the setting at the convergence of Fig. 2. Same parameters as in Fig. 2 are
adopted. Figure 4: Graphical solution of the equation
$\rho_{\pm}(P^{*})=P^{*}$. The black line represents the graph of the function
$f(x)=x$. Fixed points $p^{*}_{\pm}$ are obtained by the intersections of the
lines $\rho_{\pm}$ with the bisector. Same parameters as in Fig. 2. Figure 5:
Distribution of the node activation probability and contribution from nodes
with in-degree 1 (orange). Vertical lines indicated the range of possible
values of $P$ for nodes with in-degree 1, as predicted by Fig. 4. Same
parameters as in Fig. 2.

###  3.2 Comparison with the naive mean-field approximation

Solving the dynamics of the linear threshold model for generic couplings is
computationally demanding as already discussed above. In the context of
disordered systems, approximate methods have been proposed, which are
effective for models with dense interactions, such as the dynamic naive mean-
field (nMF) approximation and dynamic TAP equations, for interactions with an
arbitrary degree of symmetry [33], and specifically for fully asymmetric
interactions [56]. However, for diluted systems, dynamic nMF approximations
have been shown to be less effective [57, 35] for systems with both
ferromagnetic and spin-glass interactions, especially at low temperature. In
this section, we investigate the distribution of node activation and we
highlight that even in the region of parameters where the site average of the
node activation probability is well captured by the nMF method, the nMF turns
out to be incapable of reliably capturing the full distribution of node
activation probabilities.

We briefly expose the intuitive idea of the mean-field approximation, and we
refer to Ref. [35] for more details. For a node $i$ let us consider the local
field $h_{i}=\sum_{j}J_{ij}n_{j}$ and the probability associated
$\mathrm{Prob}(h_{i})=\langle\delta(h_{i}-\sum_{j}J_{ij}n_{j})\rangle_{\bm{n}}$.
The core idea of the mean-field approximation is to approximate
$\mathrm{Prob}(h_{i})$ with a delta-function peaked on the mean value $\langle
h_{i}\rangle$, which implies
$\left\langle\Phi_{\beta}\left(h_{i}-\vartheta_{i}\right)\right\rangle\approx\Phi_{\beta}\left(\left\langle
h_{i}\right\rangle-\vartheta_{i}\right)$. This assumption is motivated in the
limit where the number of neighbours of node $i$ diverges where, thanks to the
law of large numbers, one expects $\mathrm{Prob}(h_{i})$ to concentrate — for
suitably normalized $J_{ij}$ — around the mean value. However, in the case of
sparse matrices, the number of neighbours is typically small, therefore, we
expect the nMF to fail. However, our results below indicate that the site
averaged node activation is well captured by the mean-field approximation in a
variety of conditions.

Starting from Equation (2), the associated mean-field equation for fully
asymmetric networks is

$P_{i}(t+1)=\Phi_{\beta}\left(\sum_{j}J_{ij}P_{j}(t)-\vartheta_{i}\right)\,,$
(16)

which represents an approximation of Equation (7). We run Equation (16) for
each node $i$ starting from random initial conditions until convergence. We
apply the heterogeneous mean-field in the same class of networks as above
(fully asymmetric adjacency matrix with discrete coupling strength) and use
our algorithm as ground truth to evaluate the performance of the mean-field
approximation. We compute the site average of the node activation probability
$\langle P\rangle=N^{-1}\sum_{i}P_{i}$ both in the case where the set
$\\{P_{i}\\}_{i=1}^{N}$ is computed by using Equation (16) and Equation (7).
In Fig. 6 we compare the results for the site average of node activation at
different values of bias $\eta$ and threshold $\vartheta$ for two network
models, namely the directed random regular graphs, and networks with fat-
tailed degree distribution. Results indicate the mean-field captures
remarkably well the behaviour of the site average for different choices of
network models and both in the absence and presence of an external threshold.
However, the mean-field approximation captures poorly the heterogeneous site
dependence, and the distribution $\Pi(P)$ presents large differences between
the cavity and mean-field results.

Figure 6:  Comparison between naive mean-field (nMF) and cavity (cav)
solutions at $T/J=0.2$ in random directed networks with $N=100000$ nodes and
regular topology with $k^{\mathrm{in}}=k^{\mathrm{out}}=3$ (left) and
heterogeneous in- and out-degrees, drawn from a fat-tailed distribution with
$\gamma=4$ (right). Upper panels show the average activation probability
$\langle P\rangle=N^{-1}\sum_{i=i}^{N}P_{i}$ versus the bias $\eta$, for
$\vartheta/J=0,0.5$. Results from cavity are shown by triangle markers, while
squares show mean-field results. Lower panels show the distribution $\Pi(P)$
of the site-dependent activation probabilities $\\{P_{i}\\}_{i=1}^{N}$, for
$\vartheta=0$, $\eta=0.7$. Orange lines show results from cavity Equation (7),
blue lines show results from naive mean-field Equation (16).

##  4 Pairwise correlation for the linear threshold model on fully asymmetric
networks

Temporal and pairwise correlations are of great interest for the analysis of
dynamics on networks. For example, correlations are used extensively to
discover patterns of relationships in numerous inverse problems [58, 59, 60,
61]. In the context of spin models on networks, the correlation has been
investigated in the dense regime for asymmetric interactions through mean-
field [56] and an analytical expression for the site correlators can be
formally solved through TAP equations for generic system not at equilibrium
[33]. However, in the latter case, the analytical equations are difficult to
solve. In this section, we introduce a method to evaluate spatio-temporal
correlations in the sparse regime using dynamic cavity techniques.

As shown above, in the cavity approach the evaluation of the one-point
expectations $\langle n_{i}\rangle$ relies on the assumption that the states
associated with the predecessors of a node are independent of one another,
which is true on a tree. Instead, for a given pair of nodes, the states of the
joint predecessors are not independent of one another, even on a tree. The
cavity method has recently been applied in the presence of node correlations
on loopy graphs to investigate percolation, random matrix spectra [62], and
the equilibrium state of Ising spin models [63]. The authors explicitly
include the contribution from short loops, e.g., triangles, which break the
assumption of independence of states the cavity method relies on. Here, we
solve a different problem, namely, we adapt the cavity method to evaluate two-
points expectations. In our context, a new set of motifs that break the
assumption of independence of states arises, which are present even in
directed trees. Therefore, we devise a method that incorporates node
correlations inside the cavity equations.

If the graph admits a giant connected component, every set of nodes belonging
to the giant out-component shares at least one common ancestor, which would
make the evaluation of node correlation infeasible on a large graph. In the
following, we assume correlations to be mainly generated by two network
motifs, where nodes have common predecessors either at distance one or two,
and we test our approximate results against direct simulations.

Let us consider the stationary (time-translation invariant) pairwise connected
correlation (also called covariance)

$C_{ij}^{\mathrm{c}}(\tau)\coloneqq\lim_{t\rightarrow\infty}\,\langle
n_{i}(t+\tau)n_{j}(t)\rangle-P_{i}P_{j}\,,$ (17)

that is defined between any pair of sites $\\{i,j\\}$ and time-lag $\tau$. The
empirical stationary pairwise connected correlation (as computed from Monte-
Carlo simulations) is analogously defined as

$\hat{C}_{ij}^{\mathrm{c}}(\tau)\coloneqq\frac{1}{t_{s}}\sum_{t=t_{0}}^{t_{0}+t_{s}}n_{i}(t+\tau)n_{j}(t)-P_{i,MC}P_{j,MC}\,.$
(18)

for $t_{s}\gg 1$, and $t_{0}$ a sufficiently long time to allow the dynamics
to relax to stationarity. The covariance $C_{ij}^{\mathrm{c}}(\tau)$ depends
both on the sites and on the lag. In the following we restrict ourselves to
two levels of analysis, namely the single-site lagged auto-covariance
$C_{ii}^{\mathrm{c}}(\tau)$, and the covariance at zero time lag,
$C_{ij}^{\mathrm{c}}(0)$. We compare analytical results for those two types of
covariance against simulations.

#### Single-site lagged auto-covariance

For the dynamics described by Equation (2), we first investigate the lagged
auto-covariance $\hat{C}_{ii}^{\mathrm{c}}(\tau)$. We will present and discuss
results in term of the time-lagged local Pearson correlation defined as,
$\rho_{ii}(\tau)\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\hat{C}_{ii}^{\mathrm{c}}(\tau)/\hat{C}_{ii}^{\mathrm{c}}(0)$,
as well as its sample average
$\langle\rho(\tau)\rangle=N^{-1}\sum_{i}\rho_{ii}(\tau)$. Simulation results
show that $\rho_{ii}(\tau)\approx 0$ for $\forall i,\forall\tau>0$ (hence
trivially $\langle\rho(\tau)\rangle\approx 0$), which means that the dynamics
presents no retarded auto-covariance, see Fig. 7. This is consistent with the
theoretical prediction from the cavity formalism on fully asymmetric networks
in the absence of loops, $C_{ii}^{\mathrm{c}}(\tau)=0$ $\forall
i,\forall\tau>0$. Heuristically, non-zero auto-covariance is associated with
feedback loops which send the information back to the starting node with a
delay. However, on a directed tree this never happens.

Figure 7:  Single-site lagged Pearson autocorrelation $\rho_{ii}(\tau)$ for
different nodes $i$ (left) and its site average $\langle\rho(\tau)\rangle$
(right) as a function of the time-lag $\tau$ in a directed random regular
graph with parameters: $N=10000$, $k^{\mathrm{out}}=k^{\mathrm{in}}=2$,
$T/J=0.2,\vartheta=0$.

#### Pairwise connected correlation

Below we compute the analytical pairwise covariance at zero time lag,
$C_{ij}^{c}(0)$, which we refer to as $C_{ij}^{\mathrm{c}}$ for simplicity.
Correlations between two nodes $i,j$ are expected to arise from common
ancestors of the two nodes. Given the earlier result of the absence of single-
site lagged auto-covariance, we expect that if the common ancestor $k$ is not
at the same distance from $i$ and $j$, signals would need a different amount
of time to travel from $k$ to $i$ and $j$, respectively, entailing that
$C_{ij}^{c}(\tau)\neq 0$ only at a non-zero $\tau$ corresponding to the
difference in distances. Here we are interested in the equal-time connected
correlation, $C_{ij}^{c}(0)$, which is thus non-zero only if the common
ancestor $k$ is at the same (finite) distance from nodes $i$ and $j$.
Therefore, network motifs that produce non-zero values of the pairwise equal-
time covariance are required to have their common root node at equal distance
from the leaf nodes; see for example the motifs that are shown in Fig.
8(a)-(b). The (equal-time) connected correlation is obtained from the joint
expectation

$C_{ij}=\left\langle
n_{i}n_{j}\right\rangle=\left\langle\Phi_{\beta}\left(h_{i}(\bm{n}_{\partial_{i}})-\vartheta_{i}\right)\Phi_{\beta}\left(h_{j}(\bm{n}_{\partial_{j}})-\vartheta_{j}\right)\right\rangle_{\bm{n}_{\partial_{i}}\cup\,\bm{n}_{\partial_{j}}}\,,$
(19)

via $C^{\mathrm{c}}_{ij}=C_{ij}-P_{i}P_{j}$, where
$\bm{n}_{\partial_{i}}\cup\,\bm{n}_{\partial_{j}}$ indicates the set of states
associated with the predecessors of nodes $i$ or $j$. In the following we
consider contributions to the correlations coming from the two types of motifs
shown in Fig. 8(a-b). They are (a) a pair of nodes sharing one predecessor,
and (b) a pair of direct descendants from the nodes in motif (a). We show that
those motifs provide the dominant contributions to the $C_{ij}$ distribution.
If we use the standard cavity approximation

$P(\bm{n}_{\partial_{i}}\cup\,\bm{n}_{\partial_{j}},t)\approx\prod_{\ell\in\partial_{i}\cup\,\partial_{j}}P(n_{\ell},t)\,,$
(20)

which assumes that the joint probability of node activations of the
predecessors of nodes $i$ or $j$ factorises over single-node terms, we are
able to take into account the effect of motif (a) (but only motif (a)). This
is because the nodes belonging to the set of predecessors of both $i$ and $j$,
$\\{\ell,\ell\in\partial_{i}\cap\partial_{j}\\}$, appear only once in the
product on the r.h.s. of Equation (20). This then explains why
$C_{ij}^{\mathrm{c}}\neq 0$ in these cases. We evaluate Equation (19) in the
approximation that takes motif (a) into account through dynamic programming in
a similar fashion as done above in the context of the node activation
probability. Let $K=|\partial_{i}\cup\,\partial_{j}|$ be the number of joint
predecessors of nodes $i$ or $j$ and, to ease the notation, let us use
$\\{1,\dots K\\}$ to enumerate the nodes in this set of predecessors. We
define the family of recursive functions
$C_{ij}^{(1)}(\ell,\tilde{h}_{i},\tilde{h}_{j})$ such that $\ell$ represents
the node we perform the average over. If node $\ell$ is active, it gives a
contribution $J_{i\ell}$ and $J_{j\ell}$ to the auxiliary fields
$\tilde{h}_{i}$ and $\tilde{h}_{j}$, respectively. Otherwise the auxiliary
fields do not change. The dynamic programming evaluation for $C_{ij}$ becomes

$\displaystyle C_{ij}$ $\displaystyle=C^{(1)}_{ij}(1,0,0)$ (21) $\displaystyle
C_{ij}^{(1)}(\ell,\tilde{h}_{i},\tilde{h}_{j})$
$\displaystyle=P_{\ell}C_{ij}^{(1)}(\ell+1,\tilde{h}_{i}+J_{i\ell},\tilde{h}_{j}+J_{j\ell})+(1-P_{\ell})C_{ij}^{(1)}(\ell+1,\tilde{h}_{i},\tilde{h}_{j})$
$\displaystyle\mbox{for }1\leq\ell\leq K$ $\displaystyle
C_{ij}^{(1)}(K+1,\tilde{h}_{i},\tilde{h}_{j})$
$\displaystyle=\Phi_{\beta}\left(\tilde{h}_{i}-\vartheta_{i}\right)\Phi_{\beta}\left(\tilde{h}_{j}-\vartheta_{j}\right)\,$

where $P_{\ell}$ is the stationary probability of node $\ell$. From the
connected correlation we obtain the Pearson correlation coefficient

$\rho_{ij}=\frac{C_{ij}^{\mathrm{c}}}{\sqrt{(P_{i}-P_{i}^{2})(P_{j}-P_{j}^{2})}}\,,$
(22)

and we compare the values $\rho_{ij}$ obtained through cavity approximation
with that obtained from the direct evaluation of the microscopic dynamics in
Fig. 8(c). For the reasons explained above, we obtain non-zero values of
$\rho_{ij}$ if and only if nodes $i$ and $j$ have a common predecessor $\ell$,
as in the motif shown in Fig. 8(a). Apart from the trivial peak at 1,
resulting from the diagonal terms of the distribution, the largest absolute
values of $\rho$ are captured nearly perfectly by the motifs illustrated in
Fig. 8(a). On the other hand, such motifs do not capture well the events
giving rise to small correlations.

In order to improve the characterisation of these events, for every pair of
nodes $\\{\ell_{i},\ell_{j}\\}$ that are described by the motif (a), we
investigate the connected correlation associated with the pair of nodes
$\\{i,j\\}$, with $i$ and $j$ a successor of $\ell_{i}$ and $\ell_{j}$
respectively. This motif is illustrated in Fig. 8(b) and consists of a pair of
nodes $\\{i,j\\}$ that has at least one common ancestor $m$ at distance two.
In the following, we use the estimate of the joint expectation for
$\\{\ell_{i},\ell_{j}\\}$ provided by Equation (21) to estimate the joint
expectation of the pair $\\{i,j\\}$. For the two nodes denoted $\ell_{i}$ and
$\ell_{j}$ in the Figure, there is a non-zero correlation,
$\rho_{\ell_{i}\ell_{j}}\neq 0$ as shown above, hence the joint probability of
the states of the predecessors of $i$ or $j$ does not factorise, and Equation
(20) does not hold anymore. For a given pair $\\{i,j\\}$, let us denote the
set of pairs of their predecessors $\\{\\{\ell_{i},\ell_{j}\\}\\}$ which have
a predecessor in common by
$L_{1}=\\{\\{\ell_{i},\ell_{j}\\}:\ell_{i}\in\partial_{i},\ell_{j}\in\partial_{j}\mbox{
and }\partial_{\ell_{i}}\cap\partial_{\ell_{j}}\neq\varnothing\\}$, and let
$L_{0}$ be set of predecessors of $i$ or $j$ which do not appear in $L_{1}$
and are treated as independent. Instead of Equation (20), the marginal
probability of the joint set of neighbours contains the two-points joint
probability of the pairs of states in the set $L_{1}$

$P(\bm{n}_{\partial_{i}}\cup\,\bm{n}_{\partial_{j}},t)\approx\prod_{\ell\in
L_{0}}P(n_{\ell},t)\prod_{\ell_{i},\ell_{j}\in
L_{1}}P(n_{\ell_{i}},n_{\ell_{j}},t)\,.$ (23)

For each pair of sites in $L_{1}$, the joint probability $P(n_{a},n_{b},t)$ is

$P(n_{a},n_{b},t)=\Big{\langle}\left[n_{a}\Phi_{\beta}\left(h_{a}(\bm{n}_{\partial_{a}}))\right)+(1-n_{a})(1-\Phi_{\beta}\left(h_{a}(\bm{n}_{\partial_{a}}))\right)\right]\\\
\left[n_{b}\Phi_{\beta}\left(h_{b}(\bm{n}_{\partial_{b}})\right)+(1-n_{b})(1-\Phi_{\beta}\left(h_{b}(\bm{n}_{\partial_{b}})\right)\right]\Big{\rangle}_{{\bm{n}_{\partial_{a}}\cup\,\bm{n}_{\partial_{b}}},t-1}\,.$
(24)

To simplify the notation, we drop the explicit time dependence in
$P(n_{a},n_{b},t)$; the result can be written in terms of $C_{ab},P_{a},P_{b}$
as

$P(n_{a},n_{b})=(1-n_{a})(1-n_{b})+P_{a}(2n_{a}-1)(1-n_{b})+P_{b}(2n_{b}-1)(1-n_{a})\\\
+C_{ab}(2n_{a}-1)(2n_{b}-1)\,,$ (25)

and $C_{ab}$ is computed from Equation (21). We now evaluate Equation (19)
using Equation (23), which also takes into account motifs (b), in contrast to
Equation (20) which only accounts for motifs (a). The correction applies to
the evaluation of the entries $C_{ij}$ for which a node $m$ exists which is an
ancestor of both nodes $i$ and $j$ through the paths $(i,\ell_{i},m)$ and
$(j,\ell_{j},m)$; see motif in Fig. 8 (b). Let $K_{0}=|L_{0}|$ and
$K_{1}=|L_{1}|$, and let $L_{0}=\\{1,\dots K_{0}\\}$ and
$L_{1}=\\{\\{a_{1},b_{1}\\},\dots\\{a_{K_{1}},b_{K_{1}}\\}\\}$. Equation (21)
is replaced by the dynamic programming relationships

$\displaystyle C_{ij}$ $\displaystyle=C_{ij}^{(1)}(1,0,0)$ (26) $\displaystyle
C_{ij}^{(1)}(\ell,\tilde{h}_{i},\tilde{h}_{j})$
$\displaystyle=P_{\ell}C_{ij}^{(1)}(\ell+1,\tilde{h}_{i}+J_{i\ell},\tilde{h}_{j}+J_{j\ell})+(1-P_{\ell})C_{ij}^{(1)}(\ell+1,\tilde{h}_{i},\tilde{h}_{j})$
$\displaystyle\mbox{for }1\leq\ell\leq K_{0}$ $\displaystyle
C_{ij}^{(1)}(K_{0}+1,\tilde{h}_{i},\tilde{h}_{j})$
$\displaystyle=C^{(2)}_{ij}(1,0,0)$ $\displaystyle
C^{(2)}_{ij}(\ell,\tilde{h}_{i},\tilde{h}_{j})$
$\displaystyle=P(n_{a_{\ell}}=1,n_{b_{\ell}}=1)C_{ij}^{(2)}(\ell+1,\tilde{h}_{i}+J_{ia_{\ell}},\tilde{h}_{j}+J_{jb_{\ell}})$
$\displaystyle~{}~{}+P(n_{a_{\ell}}=1,n_{b_{\ell}}=0)C_{ij}^{(2)}(\ell+1,\tilde{h}_{i}+J_{ia_{\ell}},\tilde{h}_{j})$
$\displaystyle~{}~{}+P(n_{a_{\ell}}=0,n_{b_{\ell}}=1)C_{ij}^{(2)}(\ell+1,\tilde{h}_{i},\tilde{h}_{j}+J_{jb_{\ell}})$
$\displaystyle~{}~{}+P(n_{a_{\ell}}=0,n_{b_{\ell}}=0)C_{ij}^{(2)}(\ell+1,\tilde{h}_{i},\tilde{h}_{j})$
$\displaystyle\mbox{for }1\leq\ell\leq K_{1}$ $\displaystyle
C^{(2)}_{ij}(K_{1}+1,\tilde{h}_{i},\tilde{h}_{j})$
$\displaystyle=\Phi_{\beta}(\tilde{h}_{i}-\vartheta_{i})\Phi_{\beta}(\tilde{h}_{j}-\vartheta_{j})\,,$

where the two point marginals $P(n_{a},n_{b})$ are computed from Equation (25)
using Equation (21) to compute $C_{ab}$. The comparison between simulations
and the analytical calculation is shown in Fig. 8(d). The results indicate
that the motifs in Figs. 8(a) and (b) very accurately capture the salient
properties of pairwise correlations in this system. This finding suggests that
methods for network reconstruction from experimentally measured correlations
may have limited ability to characterise motifs that involve common ancestors
at distance larger than two.

(a)

(b)

(c)

(d)

Figure 8: (a)-(b) Network motifs that contributes to correlation. (c)-(d)
Probability density of the Pearson correlation coefficient $\rho$. Only
elements of the upper triangular matrix (including the diagonal) of $\rho$ are
taken into account for the histogram. Parameters: $k^{\mathrm{in}}=2$.
$N=5000$, $T/J=0.4$, $\vartheta=0$.

##  5 Extensions

###  5.1 Multi-node interaction model

The linear threshold model presented above in Equation (9) describes a dynamic
model consisting of pairwise interactions between nodes. However, approaches
based on pairwise interactions may not always be appropriate [64, 65].
Therefore several models that explicitly consider multi-node interactions
(called high-order models) have been proposed, which can be implemented in
terms of hypergraphs [66], simplicial complexes [67], and bipartite or factor
graphs.

Let us consider a directed bipartite graph $G(V,E)$, where $V$ denotes the set
of nodes and $E$ the set of edges. The set of nodes $V$ is composed of two
disjoint sets $X$ and $Y$, with no edges in $E$ having both endpoints in $X$
or $Y$. Let $|X|=N$, and $|Y|=\alpha N$. We refer to the $N$ elements of $X$
as variable nodes and the $\alpha N$ elements of $Y$ as function nodes. To
define a dynamic process on the bipartite network, we associate a dynamic
variable with every node of $X$. We use $\\{n_{i}\\}_{i=1}^{N}$ and
$\\{\tau_{\mu}\\}_{\mu=1}^{\alpha N}$ to denote the dynamic variables and the
function nodes respectively. Every function node $\mu$ is associated with a
Boolean function $g_{\mu}(\cdot)$, with the function
$g_{\mu}(\bm{n}_{\partial_{\mu}})$ depending on the states of $\mu$’s
predecessors $\bm{n}_{\partial_{\mu}}$. Every variable node $i$ evolves
according to a linear threshold function depending on node $i$’s predecessors
$\bm{\tau}_{\partial_{i}}$. The bipartite graph is defined by the bi-adjacency
matrix ${\bf A}=(A_{\mu i})$ and the interaction matrix ${\bf J}=(J_{i\mu})$,
with entries independent of $\bf{A}$. The matrix ${\bf A}=(A_{\mu i})$ has
dimension $\alpha N\times N$; a value $A_{\mu i}=1$ indicates that a link from
the variable node $i$ to the function node $\mu$ exists, signifying the fact
that the state $n_{i}$ is an argument of the function $g_{\mu}$, and $A_{\mu
i}=0$ otherwise. The interaction matrix ${\bf J}=(J_{i\mu})$ has dimension
$N\times\alpha N$; a non-zero entry $J_{i\mu}$ represents the strength of the
influence of the function $\mu$ on the variable node $i$. For any realisation
of the matrices ${\bf A}$ and ${\bf J}$, we can define $d_{i}({\bf J})$ as the
number of predecessors of node $i$ and $c_{\mu}({\bf A})$ the number of
predecessors of node $\mu$:

$\displaystyle d_{i}({\bf J})=\sum_{\mu}^{\alpha
N}\Theta\left(|J_{i\mu}|\right)\,,\qquad c_{\mu}({\bf A})=\sum_{i}^{N}A_{\mu
i}\,.$ (27)

We consider the dynamics where the state of the function node $\mu$, given by
the binary variable $\tau_{\mu}$, is determined by
$g_{\mu}(\bm{n}_{\partial_{\mu}})$. Conversely, the binary variables $n_{i}$
follow a linear threshold dynamics formulated in terms of the $\tau_{\mu}$,
i.e.

$\displaystyle n_{i}(t+1)$
$\displaystyle=\Theta\left(\sum_{\mu}J_{i\mu}\tau_{\mu}(t)-\vartheta_{i}-z_{i}(t)\right)$
(28) $\displaystyle\tau_{\mu}(t)$
$\displaystyle=g_{\mu}\left(\bm{n}_{\partial_{\mu}}(t)\right)\,.$

where the $z_{i}(t)$ represent independent identically distributed stochastic
noise terms, whose statistics are defined in terms of their cumulative
distribution function $\Phi_{\beta}(z)$ as before. The update rule in Equation
(28) represents a more flexible model than the pairwise interaction model
defined in Equation (2). In Fig. 9 (a) we show an example of a small
monopartite network that captures pairwise interactions as described by
Equation (2): a node $i$ receives the signal from the two nodes $j$ and $k$
through distinct links. Instead, the bipartite graph model shown in Fig 9 (b)
is more flexible, as the function nodes return arbitrary Boolean functions
$g_{\mu}(\cdot),g_{\nu}(\cdot),g_{\lambda}(\cdot)$ of their respective inputs.
These can represent the standard pairwise interactions through function nodes
such as node $\mu$ in Fig. 9(b), which is a function of a single variable
$n_{j}$, but also generic multi-node interactions such as through node
$\lambda$, which depends on two Boolean variables $n_{j}$ and $n_{k}$.

In what follows, we develop the formalism to solve the dynamics of bipartite
graphs of the type represented in Fig. 9 (b), and we show that these systems
benefit from the speedup of the algorithm discussed above. We will
specifically look at an example of a system with multi-node interactions that
is motivated by gene regulations [50].

(a)

(b)

Figure 9:  Snapshot of a network describing a pairwise model; each link
carries the signal from one node only (left). A bipartite graph is used to
represent multi-node interaction, here the node $i$ receives the signal from
nodes $j,k$ through function nodes $\mu,\nu$ and $\lambda$, with node
$\lambda$ representing a Boolean function of two variables (right). Arrows
indicate a link pointing from a function node to a node, square markers join a
function node with the arguments of the function.

In the following we investigate the case
$g_{\mu}(\bm{n})=\prod_{i\in\partial_{\mu}}n_{i}$, which represents an “AND”
gate logic,

$\displaystyle n_{i}(t+1)$
$\displaystyle=\Theta\left(\sum_{\mu}J_{i\mu}\tau_{\mu}(t)-\vartheta_{i}-z_{i}(t)\right)$
(29) $\displaystyle\tau_{\mu}(t)$
$\displaystyle=\prod_{j\in\partial_{\mu}}n_{j}(t)\,.$

Let $P_{i}(t)$ be the probability that state $n_{i}(t)=1$. Let
$\mathcal{P}_{\mu}(t)\coloneqq\mathrm{Prob}\left[\tau_{\mu}(t)=1\right]$ be
the probability that state $\tau_{\mu}(t)=1$. If the directed bipartite graph
does not contain short loops, the dynamic cavity method provides an efficient
approximate solution of the dynamics. In particular the node activation
probabilities are given by the set of expressions

$\displaystyle P_{i}(t+1)$
$\displaystyle=\left\langle\Phi_{\beta}\left(h_{i}(\boldsymbol{\tau}_{\partial_{i}})-\vartheta_{i}\right)\right\rangle_{\bm{\tau}_{\partial_{i}}}$
(30) $\displaystyle\mathcal{P}_{\mu}(t)$
$\displaystyle=\prod_{j\in\partial_{\mu}}P_{j}(t)$ (31)

with $h_{i}(\bm{\tau}_{\partial_{i}})=\sum_{\mu}J_{i\mu}\tau_{\mu}$. Equation
(30) is of the same form as Equation (9) and thus the dynamic programming
approach can speed-up the evaluation of the average. Let us consider a node
$i$ with in-degree $d_{i}$, and let $\nu\in\\{1,\dots d_{i}\\}$ be the indices
of the predecessors of node $i$. We define a function $f_{i}(\nu,\tilde{h})$
recursively

$\displaystyle f_{i}(\nu,\tilde{h})$
$\displaystyle=\mathcal{P}_{\nu}(t)f_{i}(\nu+1,\tilde{h}+J_{i\nu})+(1-\mathcal{P}_{\nu}(t))f_{i}(\nu+1,\tilde{h})$
$\displaystyle\mbox{for }1\leq\nu\leq d_{i}$ (32) $\displaystyle
f_{i}(d_{i}+1,\tilde{h})$
$\displaystyle=\Phi_{\beta}\left(\tilde{h}-\vartheta\right)\ .$

In terms of this recursion, we have

$P_{i}(t+1)=f_{i}(1,0)\ .$ (33)

In analogy with the pairwise model, we assume that
$J_{i\mu}\in\\{-J_{i},0,J_{i}\\}$ with $J_{i}$ a real positive number that may
depend on the variable node. Under this condition, the dynamic programming
algorithm provides a speedup of the evaluation of node activation
probabilities.

#### Node activation probability with multi-node interaction

As for the linear threshold model, we inspect the stationary state of the node
activation probabilities at different values of the noise intensity; results
are shown in Fig. 10. In our investigation, we use the same degree
distribution of the in- and out-degree as in the monopartite network, both for
the variable and function nodes. The variable $\eta$ refers to the fraction of
positive terms in the interaction matrix $\bm{J}$ in analogy to the
monopartite case.

Dominant peaks of the node activation probability can be rationalised in terms
of the discrete stochastic map of Equation (15). This should not come as a
surprise: given our choice of the degree distribution, the majority of nodes
only have a pairwise interaction with a single predecessor node. Therefore the
same argument used to derive the theoretical lines of Equation (15), which
describes the temperature dependence of some of the prominent peaks in the
distribution of node activation probabilities, still holds. However, the
distribution shown in Fig. 10 is substantially different from the distribution
obtained in the linear threshold model at zero thresholds, which was analysed
in Ref. [47], highlighting contributions due to multi-node interactions. The
effect of multi-node interactions is particularly easy to analyse in the high
noise regime, which is discussed in A.

Figure 10:  Heat-map of the distribution of node activation probabilities
$\Pi(P)$ at different values of noise parameter $T$ resulting from the
microscopic dynamics detailed in Equation (29) (top). Dashed, dotted, and dot-
dashed lines show the graph of different function compositions of $\rho_{\pm}$
at different values of $T$, see Equation (15). Vertical dotted lines mark the
noise levels at which the histograms of $\Pi(P)$ are shown in the bottom
panels. They probe the regimes at low, intermediate, and high noise
corresponding to $T/J=0.1,0.4,0.8$, from left to right. Parameters:
$N=200000$, $\alpha=1$, $\gamma=2.81$, $\vartheta=0$.

###  5.2 Any Boolean function can be written in terms of AND gates

We have demonstrated above, how our linear programming method can be extended
to analyse the stochastic dynamics of multi-node interaction models, using a
formalisation in terms of bipartite graphs in which all function nodes are
logical AND gates, i.e., we have
$\tau_{\mu}=g_{\mu}(\bm{n}_{\partial_{\mu}})=\prod_{j\in\partial_{\mu}}n_{j}$
for all $\mu$. We now show that the bipartite model with AND gates is
universal in the sense that a local field of the form
$h_{i}=\sum_{\mu}J_{i\mu}\tau_{\mu}$ as it appears in Equation (29) can
represent arbitrary functions of $k_{i}=|\partial_{i}|$ Boolean variables. To
see this, consider an arbitrary real function
$g:\\{0,1\\}^{k}\rightarrow\mathbb{R}$ of $k$ Boolean variables. It can always
be written using an expansion of the form

$g(\bm{n})=a^{(0)}+\sum_{i}^{k}a^{(1)}_{i}n_{i}+\sum_{i<j}^{k}a^{(2)}_{ij}n_{i}n_{j}+\sum_{i<j<\ell}^{k}a^{(3)}_{ij\ell}n_{i}n_{j}n_{\ell}+\dots+a^{(k)}\prod_{\ell=1}^{k}n_{\ell}\
,$ (34)

in which terms are arranged by increasing order, with suitable coefficients
$a^{0},\bm{a}^{(1)},\bm{a}^{(2)},\dots,a^{(k)}$, which can be found
recursively by considering values of $g(\bm{n})$ successively for
configurations $\bm{n}$, with non-zero components on subsets of
$\\{1,\dots,k\\}$ of increasing size. The local field $h_{i}$ in Equation (29)
is exactly of this form if the $\mu\in\partial_{i}$, i.e. the $\mu$ for which
$J_{i\mu}\neq 0$, are arranged in order of increasing $|\partial_{\mu}|$.

###  5.3 Spin dynamics

Let $\sigma_{i}$, $i\in\\{1,\dots N\\}$, denote a set of Ising spin variables,
taking values in $\\{-1,1\\}$, located on the $N$ nodes of a complex network,
and let $J_{ij}$ denote the strength of the interaction from node $j$ to node
$i$. The Ising spin dynamics is defined as

$\sigma_{i}(t+1)=\mathrm{sgn}\left[\sum_{j}^{N}J_{ij}\sigma_{j}(t)-\vartheta_{i}-z_{i}(t)\right],$
(35)

where $\vartheta_{i}$ is a threshold, and $z_{i}(t)$ is an independent
identically distributed random variable defined as above with
$\mathrm{Prob}[z<x]=\Phi_{\beta}(x)$.444Note that, in analogy with the rest of
the manuscript, the variables $\vartheta_{i}$ represent a threshold and not an
external field. The average of $\mathrm{sgn}(x-z)$ over the noise distribution
can be expressed in terms of the cumulative distribution function of the noise
as

$\tilde{\Phi}_{\beta}(x)\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\langle\mathrm{sgn}(x-z)\rangle_{z}=2\Phi_{\beta}(x)-1\,.$
(36)

Let us denote by
$m_{i}(t)\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\sum_{\sigma}P_{i}(\sigma,t)\sigma=\langle\sigma_{i}\rangle_{\sigma_{i},t}$
the time dependent average of the Ising spin at node $i$. For a fully
asymmetric graph, the cavity method provides an expression for the
magnetisation in a similar form as above. Performing the average in (35) gives

$m_{i}(t+1)=\left\langle\tilde{\Phi}_{\beta}\left(h_{i}(\bm{\sigma}_{\partial_{i}})-\vartheta_{i}\right)\right\rangle_{\bm{\sigma}_{\partial_{i}},\,t}\
,$ (37)

in which $h_{i}(\bm{\sigma}_{\partial_{i}})=\sum_{j}J_{ij}\sigma_{j}$ and
$\bm{\sigma}_{\partial_{i}}$ plays a role analogous to Equation (8) for spin
variables. In analogy to the linear threshold model, Equation (35) can be
evaluated efficiently through the recursion

$\displaystyle m_{i}(t+1)$ $\displaystyle=f_{i}(1,0)$ (38) $\displaystyle
f_{i}(\ell,\tilde{h})$
$\displaystyle=\dfrac{1+m_{\ell}(t)}{2}f_{i}(\ell+1,\tilde{h}+J_{i\ell})+\dfrac{1-m_{\ell}(t)}{2}f_{i}(\ell+1,\tilde{h}-J_{i\ell})$
$\displaystyle\mbox{for }1\leq\ell\leq k_{i}$ $\displaystyle
f_{i}(k_{i}+1,\tilde{h})$
$\displaystyle=\tilde{\Phi}_{\beta}\left(\tilde{h}-\vartheta_{i}\right)\,.$

####  Mixed $p$-spin model

We briefly specify the expression for the p-spin model, whose equilibrium
properties have long been investigated in the literature [68]. Let us consider
a system defined on a bipartite network with state variables
$\sigma_{i}\in\\{\pm 1\\}$ and $\tau_{\mu}\in\\{\pm 1\\}$. The update rule is

$\displaystyle\sigma_{i}(t+1)$
$\displaystyle=\mbox{sign}\left(\sum_{\mu}J_{i,\mu}\tau_{\mu}(t)+z_{i}(t)-\vartheta_{i}\right)\,$
(39) $\displaystyle\tau_{\mu}(t)$
$\displaystyle=\prod_{j\in\partial_{\mu}}\sigma_{j}(t)\,.$

Let $m_{\mu}(t)=\left\langle\tau_{\mu}\right\rangle_{\tau_{\mu},t}$, and
$m_{i}(t)=\left\langle\sigma_{i}\right\rangle_{\sigma_{i},t}$ denote the time
dependent averages of the $\tau_{\mu}$ and the $\sigma_{i}$, respectively. The
time evolution of the local magnetizations $m_{i}(t)$ is given by

$\displaystyle m_{i}(t+1)$
$\displaystyle=\left\langle\tilde{\Phi}_{\beta}\left(h_{i}(\boldsymbol{\tau}_{\partial_{i}})-\vartheta_{i}\right)\right\rangle_{\bm{\tau}_{\partial_{i}},t}\,,$
(40) $\displaystyle m_{\mu}(t)$
$\displaystyle=\prod_{j\in\partial_{\mu}}m_{j}(t)\,,$

with
$h_{i}(\bm{\tau}_{\partial_{i}})=\sum_{\mu\in{\partial_{i}}}J_{i\mu}\tau_{\mu}$.
Equation (40) can be efficiently evaluated using the recursive procedure below

$\displaystyle m_{i}(t+1)$ $\displaystyle=f_{i}(1,0)$ (41) $\displaystyle
m_{i}(\nu,\tilde{h})$
$\displaystyle=\frac{1+m_{\nu}(t)}{2}f_{i}(\nu+1,\tilde{h}+J_{i\nu})+\frac{1-m_{\nu}(t)}{2}f_{i}(\nu+1,\tilde{h}-J_{i\nu})$
$\displaystyle\mbox{for }1\leq\nu\leq d_{i}$ $\displaystyle
m_{i}(d_{i}+1,\tilde{h})$
$\displaystyle=\tilde{\Phi}_{\beta}\left(\tilde{h}-\vartheta_{i}\right)$

##  6 Beyond fully asymmetric networks

Feedback loops are associated with multi-stability and oscillatory behaviour
of dynamical systems [69, 70]. These two phenomena are often observed also in
biological systems. Indeed, feedback mechanisms represent an essential
ingredient of biological networks, e.g., in the context of gene regulation
[71, 72, 73]. However, it is not yet clear what happens in the case where a
combination of many of these patterns is present in a network. The presence of
bi-directional links with opposite signs of the interactions (such as in
predator-prey systems) have been shown to produce oscillatory behaviour in the
context of continuous linear dynamics [74], but little is known in the context
of discrete state dynamics, e.g., for the linear threshold model. The dynamics
investigated in Sect. 3 is defined on tree-like fully asymmetric network
models, where feedback loops are absent.

In this section, we consider the linear threshold model of Equation (2) in the
presence of bi-directional links. Any bi-directional link can be considered as
a de-facto short loop of length two of the network. As expected, the presence
of these short loops makes the use of the dynamic cavity method more
challenging. In the cavity graph, where a node and its links have been
removed, the presence of a short loop causes a retarded self-interactions that
make the dynamics of the cavity variables non-Markovian; see Ref. [39, 38] and
B. An approximation scheme, called the one-time approximation, has been
proposed with the intent to reduce the computational complexity created by
this form of non-Markovianity [39, 75]. The approximation factorises the
memory kernel to include a one-time-step memory term. It exactly describes the
full dynamics of fully asymmetric networks and the equilibrium steady-state of
systems with symmetric interactions in the replica symmetric phase.
Furthermore, the one-time approximation has been shown to effectively capture
the stationary state also of systems with partially symmetric links in the
high noise regime [75, 55]. However, the literature has so far focused on the
agreement of macroscopic observables, such as the mean activation, while the
distribution of individual node activations has generally been overlooked. In
this section, we investigate the heterogeneous node activation in networks
with bi-directional links using the OTA and compare the results with
simulations. Among other things, our investigation also leads to an unexpected
observation: the stationary state of networks with symmetric and anti-
symmetric interactions are biased towards the active and the inactive states,
respectively.

###  6.1 The one-time-step approximation

In this section we will only provide brief account of the OTA, leaving a more
extensive discussion to the B. Let us consider the model Equation (2) in the
case where the network of interactions contains bi-directional links. For any
node $i$ in the graph, the probability of the state at time $t+1$ depends on
the states of node $i$ neighbours, $\bm{n}_{\partial_{i}}$, at the previous
time-step $t$ through the transition probability
$W[n_{i}|h_{i}(\bm{n}_{\partial_{i}}),\vartheta_{i}]$, see Equation (4). If
links are bi-directional, the states of nodes $\bm{n}_{\partial_{i}}$ at time
$t$ are not independent, because the state of every node $j\in\partial_{i}$
depends on the past history states of node $i$ through the terms
$J_{ji}n_{i}(s)$ for $s\in\\{t-1,t-3\dots\\}$.555Bi-directional links are the
only origin of loops present on a tree, thus every node depends on its own
history at even time lag. In order to describe the dynamics of the nodes
$j\in\partial_{i}$ on the cavity graph from which node $i$ and connections to
it are removed, it is useful to introduce the time-dependent cavity thresholds
$\vartheta_{j}^{(i)}(s)\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\vartheta_{j}-J_{ji}n_{i}(s)$
for $s=0,\dots t-1$ and
$\vartheta_{j}^{(i)}(0)\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\vartheta_{j}$,
which encodes the trajectory of node $i$ in terms of threshold terms, and the
local field in the cavity graph

$h_{j}^{(i)}\left(\bm{n}^{s}_{\partial_{j}}\right)\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\sum_{\ell\in\partial_{j}\setminus{\\{i\\}}}J_{j\ell}n_{\ell}^{s}\,,$
(42)

where here and in what follows we have denoted $n_{j}^{s}=n_{j}(s)$ to make
the notation more compact. Let us call $\bm{n}_{j}^{0,\dots,t}$ the trajectory
of node $j$ from time $0$ to $t$. For every bi-directional link $(i,j)$, the
dynamic cavity method provides a recursive expression for the conditional
probability
$P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t}|\bm{\vartheta}^{(i)0,\dots,t-1}_{j}\right)$
of the trajectory of node $j$ in terms of the time-dependent thresholds
$\bm{\vartheta}^{(i)0,\dots,t}_{j}=(\vartheta_{j}^{(i)}(0),\dots,\vartheta_{j}^{(i)}(t))$,
as discussed in [39, 38] and in B. The OTA assumes that, for every
$j\in\partial_{i}$, the probability of the trajectory factorises into one-time
conditional probabilities

$P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t}|\bm{\vartheta}^{(i),0,\dots,t-1}_{j}\right)\approx
P^{(i)}_{j}(n_{j}^{0})\prod_{s=0}^{t-1}P^{(i)}_{j}(n_{j}^{s+1}|\vartheta^{(i),s}_{j})\,.$
(43)

The factorisation procedure proposed in Equation (43) is not sufficient to
fully specify the cavity and marginal probability, but an additional closure
condition is required, as shown in B. In B we re-investigate the derivation of
the OTA, we show that different choices of the closure conditions lead to the
different versions of the OTA known in the literature [34, 35], and we also
propose a new version of the closure condition. It has been shown that, for
symmetric couplings, the version of Ref. [35] admits the equilibrium state as
a solution, while the version of Ref. [34] has not been designed with that
goal. We compare the performance of OTA under different closure conditions,
including the accuracy in characterising the individual node activation
probabilities, as well as the computational complexity of their
implementation. We show the version of Ref. [35] outperforms the version of
Ref. [34] at equilibrium, but the situation is reversed in the non-equilibrium
model we investigate. However, our results indicates the version of Ref. [34]
is more consistent through the different network symmetries investigated.
Based on our assessment, we expect the version of Ref. [34] to provide the
best trade-off between accuracy and ease of implementation to tackle the non-
equilibrium stationary state. The closed-form expressions for the single-time
marginal probability of the OTA we will use are 666In the following, we drop
the time dependence from the cavity thresholds as this simplifies notation and
does not lead to any ambiguity.

$P_{i}\left(n_{i}^{t}|\vartheta_{i}\right)=\sum_{n_{i}^{t-2}}\sum_{\bm{n}^{t-1}_{\partial_{i}}}W\left[n_{i}^{t}|h_{i}\left(\bm{n}_{\partial_{i}}^{t-1}\right),\vartheta_{i}\right]\left[\prod_{j\in\partial_{i}}P_{j}^{(i)}\left(n_{j}^{t-1}|\vartheta_{j}-J_{ji}n_{i}^{t-2}\right)\right]P_{i}(n_{i}^{t-2}|\vartheta_{i})\,,$
(44)

and for the cavity marginal,

$P_{j}^{(i)}\left(n_{j}^{t-1}|\vartheta_{j}^{(i)}\right)=\sum_{n_{j}^{t-3}}\sum_{\bm{n}^{t-2}_{\partial_{j}\setminus
i}}W\left[n_{j}^{t-1}|h_{j}^{(i)}\left(\bm{n}^{t-2}_{\partial_{j}}\right),\vartheta_{j}^{(i)}\right]\left[\prod_{\ell\in\partial_{j}\setminus
i}P_{\ell}^{(j)}\left(n_{\ell}^{t-2}|\vartheta_{\ell}-J_{\ell
j}n_{j}^{t-3}\right)\right]P_{j}(n_{j}^{t-3}|\vartheta_{j})\,.$ (45)

The possible values of $\vartheta_{j}^{(i)}$ that appear in Equation (45) are
$\vartheta_{j}^{(i)}=\\{\vartheta_{j},\vartheta_{j}-J_{ji}\\}$. To solve the
dynamics over the trajectory of length $t$, the initial conditions for the
marginal and cavity probabilities $P_{j}(n_{j}^{0}|\vartheta_{j})$,
$P_{j}^{(i)}(n_{j}^{0}|\vartheta_{j})$,
$P_{j}^{(i)}(n_{j}^{0}|\vartheta_{j}-J_{ji})$ are set for $\forall j$ and
$\forall i:J_{ij}\neq 0$, and $P_{j}(n_{j}^{1}|\vartheta_{j})$ is computed
from Equation (6). To evaluate the trajectory of node activation, Equation
(45) and Equation (44) are propagated over the time-steps of interests
$0,\dots t$. At every time iteration, both Equation (45) and Equation (44)
would face a complexity barrier created by a large in-degree of node $i$ and
$j$ respectively, and their evaluation benefit from our dynamic programming
implementation. To this end, we note that the term
$\sum_{\bm{n}^{t-1}_{\partial_{i}}}W[\ldots][\prod_{j\in\partial_{i}}\ldots]$
in Equation (44) has the same structure as Equation (8) and can be evaluated
using dynamic programming as explained in Sec. 3.1. This also applies to the
corresponding average in Equation (45).

In this section, we apply the OTA to networks with fat-tailed degree
distributions and bi-directional links. We will investigate three different
degrees of symmetry of the interaction matrix, namely antisymmetric, symmetric
and uncorrelated interaction matrix $J$, and we highlight how the symmetry
affects the performance of the OTA in the regime of low temperature.

###  6.2 Comparison between OTA and simulations

#### Parameters of the model

We first generate an undirected network in the configuration model class
consisting of $N$ nodes. In analogy with above, degrees of the nodes are
sampled from the fat-tailed distribution
$P_{\gamma}(k)=k^{-\gamma}-(k+1)^{-\gamma}$. We then assign interaction terms
to the two directions of bi-directional links. In particular, interaction
terms may or may not have opposite sign in the two directions, i.e.
$J_{ij}J_{ji}=\pm J^{2}$. In the following, we investigate the correlation of
the sign in bi-directional links and we focus on three major symmetries:
symmetric ($J_{ij}=J_{ji}$), uncorrelated ($\langle J_{ij}J_{ji}\rangle=0$),
and antisymmetric ($J_{ij}=-J_{ji}$) interaction matrix $J$. In all settings,
the number of positive and negative interactions are statistically the same.
We iterate the OTA equations until convergence using the same criterion we
adopt on fully asymmetric networks, as we have discussed in Sect. 3.1.

#### Results

In order to compare the performance of OTA against Monte Carlo simulation, we
define $P_{MC}$ and $P_{OTA}$ the probability of node activation derived from
microscopic Monte Carlo simulation of the dynamics and from Eqs. (44)(45)
respectively at stationarity. In Fig. 11 we compare $P_{MC}$ and $P_{OTA}$ for
bi-directional interaction matrices with the three different types of
symmetries detailed above. We recall that in the symmetric setting, the
stationary state of the dynamics is an equilibrium state, and in the absence
of replica symmetry breaking the equilibrium state is described by belief
propagation [24].

Figure 11: Comparison of the one-time approximation (OTA) and Monte Carlo (MC)
simulations at stationarity in networks with bi-directional links and weights
chosen, from left to right, as uncorrelated, symmetric and anti-symmetric. Top
panels show scatter plots of single-node activation probabilities computed
from OTA versus MC results. Bottom panels show the distribution of node
activation probabilities obtained through MC and OTA. Simulations are obtained
through an average of over 500 trajectories. Parameters are $N=10000$,
$\gamma=4$, $T/J=0.5,\vartheta=0$, $t_{s}=10^{4}$ steps.

To evaluate the accuracy of the OTA, we consider the mean square difference of
$P_{OTA}$ to $P_{MC}$

$\mathrm{error}=\sqrt{\frac{1}{N}\sum_{i}\left(P_{i,\mathrm{MC}}-P_{i,OTA}\right)^{2}}\,.$
(46)

Our results suggest that the OTA is able to capture the statistics of node
activation in the stationary state remarkably well, even in the regime of
moderately low temperatures when bi-directional links are present, both for
the symmetric and the uncorrelated case. However, at still lower temperature
we expect the results from OTA to be less accurate, as observed in earlier
work for spin models [75, 35].

Our earlier analysis on fully-asymmetric networks revealed that spatial
correlations are short-ranged and auto-covariance is zero at any finite time
lag. On the other hand, the presence of correlations in the weights is
associated with time-lagged auto-correlation at non-zero lag, as shown in Fig.
12. Simulation results indicate that the antisymmetric and symmetric settings
exhibit the slowest decrease of the time-lagged auto-covariance, while the
uncorrelated setting exhibits the fastest. Nevertheless, the auto-covariance
decays quite consistently in few steps even for the symmetric and
antisymmetric interaction models, justifying the use of OTA to investigate the
stationary state. This observation suggests that a different factorisation
that includes more time-steps than in Equation (43) may provide a better
characterisation of the symmetric and antisymmetric cases.

Finally, Fig. 11 also reveals that the stationary distribution of node
activations for symmetric and antisymmetric interaction is skewed, even though
we are in the regime with zero thresholds, $\vartheta_{i}=0$ $\forall i$, and
unbiased distribution of couplings, i.e., $\langle J_{ij}\rangle=0$. We refer
to this setting as “unbiased parameter choice” in the following. This
behaviour is confirmed in Fig. 12, where mean node activation are shown at
different levels of the noise. In the following, we discuss this unexpected
result and we provide an intuitive explanation.

Figure 12: Left: Site average of the lagged auto-correlation
$\langle\rho(\tau)\rangle$ for a network with bi-directional links. Three
choices of the interaction matrix $J$ are considered: symmetric (circle),
uncorrelated (triangles), and antisymmetric (squares) interaction matrices.
The correlation is computed from Monte Carlo simulated trajectories of the
dynamics. Right: site average of the node activation at different temperatures
$\langle P\rangle=N^{-1}\sum_{i}P_{i}$ for symmetric (blue, upper curve),
uncorrelated (orange, middle curve), and antisymmetric (green, lower curve)
interaction matrices. Squares indicate the result obtained through Monte Carlo
simulations, circles are obtained through the one-time approximation. Network
and simulation parameters are the same as for Fig. 11.

##  7 Spontaneous symmetry breaking of the linear threshold model for unbiased
parameter choices

We have observed that, for an unbiased parameter choice, the linear model
defined in terms of $\\{0,1\\}$ variables exhibits a skewed distribution of
node activations in the presence of bi-directional links. This is somewhat
surprising, as in spin systems, the average magnetisation is well known to be
zero in the case of an u updatenbiased parameter choice. In this section we
compare the two models and we provide a simple example that illustrates how
such a phenomenology might.

The dynamics described in terms of indicator variables $n_{i}\in\\{0,1\\}$ can
be reformulated in terms of spin variables $\sigma_{i}\in\\{\pm 1\\}$. Recall
that the dynamics of the linear threshold model of Equation (2), defined in
terms of $N$ indicator variables, is given by

$n_{i}(t+1)=\Theta\left(\sum_{j}J_{ij}n_{j}(t)-\vartheta_{i}-z_{i}(t)\right)\,.$
(47)

Rewriting the analogous dynamics of a model of $N$ interacting spins as
defined in Equation (35) as

$\sigma_{i}(t+1)=\mathrm{sgn}\left(\sum_{j}J^{\prime}_{ij}\sigma_{j}(t)-\vartheta_{i}^{\prime}-z_{i}(t)\right)\,,$
(48)

with interaction matrix $J^{\prime}$ and thresholds
$\\{\vartheta^{\prime}_{i}\\}_{i=1}^{N}$, one can map one model onto the other
by setting

$J^{\prime}_{ij}=J_{ij}/2,\qquad\vartheta_{i}^{\prime}=\vartheta_{i}-\frac{1}{2}\sum_{\ell}J_{i\ell},\quad\forall
i,j\,.$ (49)

Unless the mapping of Equation (49) is implemented exactly, the trajectories
of the Boolean linear threshold model and the spin model will be different. In
analogy with the above, we use $m_{i}=\langle\sigma_{i}\rangle$ and
$P_{i}=\langle n_{i}\rangle$ to denote the stationary averages of spins and
indicator variables respectively. The phase diagram for spin systems has been
extensively studied at equilibrium, and it has been shown [24] that for an
unbiased parameter choice the mean magnetisation is zero, i.e.,
$N^{-1}\sum_{i}m_{i}=0$, even at low temperatures.

In a system with indicator variables instead, both Monte Carlo simulations and
an analysis using OTA show that the system exhibits spontaneous symmetry
breaking $N^{-1}\sum_{i}P_{i}\neq 1/2$ when the parameters are unbiased, as
shown in Fig. 12. In particular, at stationarity the average site activation
$N^{-1}\sum_{i}P_{i}$ is larger than $1/2$ for symmetric interaction matrices
and smaller than $1/2$ in the case of antisymmetric interaction matrices.

We present an intuitive argument that explains how the biased distribution of
node activations emerges in a model with an unbiased choice of parameters. To
this end, let us consider a mini-network composed of just two nodes $i,j$
connected by a bi-directional link, i.e. $J_{ij}J_{ji}\neq 0$. A key mechanism
in bi-directional links is the dependence of a node state on its own past,
i.e. node $i$ interacts with its neighbours $\partial_{i}$ which in turn
interact with $i$. The probability of the state of node $i$ at time $t+2$ can
be written in terms of the probability of the state of node $i$ at time $t$
using the Markov property of the microscopic dynamics,

$\displaystyle P_{i}(n_{i}^{t+2})$ $\displaystyle=\left\langle
W(n_{i}^{t+2}|J_{ij}n_{j})\right\rangle_{n_{j}\sim P(n_{j}^{t+1})}$ (50)
$\displaystyle P_{j}(n_{j}^{t+1})$ $\displaystyle=\left\langle
W(n_{j}^{t+1}|J_{ji}n_{i})\right\rangle_{n_{i}\sim P(n_{i}^{t})}\,.$

Considering the cases $J_{ij},J_{ji}\in\\{\pm J\\}$, we note that $\langle
W(1|Jn)\rangle_{n\sim P}=\rho_{+}(P)$ and $\langle W(1|-Jn)\rangle_{n\sim
P}=\rho_{-}(P)$ with $\rho_{\pm}(P)$ defined earlier in Equation (15). In
accordance with the unbiased parameters choice, the external threshold is set
to zero, $\vartheta_{i}=\vartheta_{j}=0$.

#### Symmetric network

In a symmetric network it holds $J_{ij}=J_{ji}$. Thus only the two following
cases can be realised:

  * •

ferromagnetic coupling $J_{ij}=J_{ji}=J$, leading to
$P_{i}(t+2)=\rho_{+}(\rho_{+}(P_{i}(t)))$.

  * •

antiferromagneric coupling $J_{ij}=J_{ji}=-J$, leading to
$P_{i}(t+2)=\rho_{-}(\rho_{-}(P_{i}(t)))$

In the stationary state, the node activation probability is given by the fixed
point of the map $p^{+}=\rho_{+}(p^{+})$ in the case of ferromagnetic
couplings, or $p^{-}=\rho_{-}(p^{-})$ in the case of antiferromagnetic
couplings. We now exploit the above observation by considering an unbiased
network of isolated dimers, allowing us to use the above result obtained for
individual dimers. The site average of the node activation probability gives

$\langle P\rangle=\frac{p^{+}+p^{-}}{2}=\left(2-\frac{1}{2}\tanh^{2}(\beta
J/2)\right)^{-1}>1/2\,.$ (51)

While this result is derived for a special class of networks, we expect that
for a generic unbiased symmetric network the site average activation is larger
than $1/2$, as shown by Fig. 12 where Monte Carlo simulation shows that
$\langle P\rangle>1/2$ at stationarity. Indeed, from Equation (51) this simple
argument predicts that at $T=0$, $P=2/3$ which is comparable with the
simulation result in Fig. 12 (see blue markers).

#### Antisymmetric network

An antisymmetric network defined by $J_{ij}=-J_{ji}$ is always unbiased. Let
$J_{ij}=J$ (then $J_{ji}=-J)$, Equation (50) gives

  * •

$P_{i}(t+2)=\rho_{+}(\rho_{-}(P_{i}(t)))$

  * •

$P_{j}(t+2)=\rho_{-}(\rho_{+}(P_{j}(t)))$.

In the stationary state, the node activation of node $i$ satisfies
$p_{i}=\rho_{+}(\rho_{-}(p_{i}))$, and the node activation of node $j$
satisfies $p_{j}=\rho_{-}(\rho_{+}(p_{j}))$. Again, we consider a network of
isolated dimers. The site average of the node activation probability gives

$\langle P\rangle=\frac{p_{i}+p_{j}}{2}=\left(2+\frac{1}{2}\tanh^{2}(\beta
J/2)\right)^{-1}<1/2\,,$ (52)

which suggests that in an unbiased antisymmetric network, the site average
activation is smaller than $1/2$ and it is confirmed in Fig. 12. From Equation
(52), this simple argument predicts that at $T=0$, $P=2/5$ which is comparable
with the result in Fig. 12 (see green markers).

The result of Equation (51) is an example where the dynamics of Ising spin and
indicator variables leads to macroscopically different stationary states.
While a spin model defined for unbiased parameter choices presents
$N^{-1}\sum_{i}m_{i}=0$ at equilibrium, for indicator variable we have
rationalised that in general $2\sum_{i}P_{i}-1\neq 0$.

##  8 Summary and discussion

In this paper we have solved the stochastic dynamics of Boolean linear
threshold models for a class of networks characterised by fat-tailed degree
distribution. On fully asymmetric networks, the distribution of the node
activation probabilities in the stationary state is not unimodal in general
but displays a rich structure. Salient features of the distribution can in
part be rationalised in terms of discrete stochastic maps resulting from the
underlying network structure and dynamics.

We have adapted the dynamic cavity method to study equal-time pairwise
correlations in the stationary state. For every pair of nodes, we have
computed the correlation produced by common ancestors at a distance of one or
two from the pair. We demonstrate that these two motifs capture most of the
statistics of values, indicating that pairwise correlation is effectively a
local property of the network.

Finally, we have investigated the individual node activation statistics in the
stationary state for two additional classes of models, namely for models with
multi-node interactions and for linear threshold models defined on networks
with bi-directional links. For models with multi-node interactions, we have
adopted a bipartite graph representation in which factor nodes represent “AND”
gates, and we have shown that this construction is universal in the sense that
it allows representing any Boolean function. We have previously used models of
this type to describe the regulatory mechanism of transcription factors (TFs),
and we have so far focused on the percolation problem [50, 51] as a way to
identify structural properties that gene regulatory networks need to exhibit
in order to support multi-cellular life. The present manuscript provides an
additional perspective that incorporates the role of noise and of inhibition,
which were outside the focus of our previous analysis. In the high noise
limit, we observe that TFs having a large in-degree are more robust to
spontaneous noise-induced activation, which suggests a potential use as noise
filters. On the other hand, results from percolation theory suggest that the
average in-degree of TF complexes should be small to allow sustaining a
multiplicity of attractors. Hence, the combination of those results may hint
at a trade-off between these antagonistic conditions being relevant for gene
regulatory networks.

To study systems with bi-directional links, we apply an approximation
procedure, called OTA, to analyse the stationary state. The OTA assumes a
factorisation of the memory kernel which in some sense amounts to a Markov
assumption. The reason why the OTA is found to be effective in describing the
non-equilibrium stationary state may thus be very well be related to the fact
that the lagged auto-correlation is found to be short-ranged.

The dynamics of the linear threshold model with bi-directional links yields a
somewhat surprising result: a network with an unbiased distribution of $\pm J$
interactions sustains a biased distribution of node activation probability,
both for symmetric and anti-symmetric interaction matrices. We provide a
heuristic argument that helps to explain this phenomenon. Symmetric and anti-
symmetric interactions are known to appear in biological networks as a
mechanism for positive and negative feedback mechanisms respectively [76]. The
present finding is of potential interest in the context of biological networks
and may help explain sustained activation or quiescence of nodes. Also, it
demonstrates once more how the intuition built about the Ising spins model may
be misleading if one investigates models defined on indicator variables [77].

##  9 Acknowledgements

GT is supported by the EPSRC Centre for Doctoral Training in Cross-
Disciplinary Approaches to Non-Equilibrium Systems (CANES EP/L015854/1).

References

## References

  * [1] S. F. Edwards and P. W. Anderson.  Theory of Spin Glasses.  J. Phys. F, 5:965–974, 1975.
  * [2] D. Sherrington and S. Kirkpatrick.  Solvable Model of a Spin-Glass.  Phys. Rev. Lett., 35:1792–1796, 1975.
  * [3] L. Viana and A. J Bray.  Phase diagrams for dilute spin glasses.  jpc, 18(15):3037, 1985.
  * [4] S. A Kauffman.  Metabolic stability and epigenesis in randomly constructed genetic nets.  J. Theor. Biol, 22(3):437–467, 1969.
  * [5] B. Derrida and Y. Pomeau.  Random networks of automata: a simple annealed approximation.  Europhys. Lett., 1(2):45, 1986.
  * [6] G. Parisi.  A simple model for the immune network.  Proc. Natl. Acad. Sci., U. S. A, 87(1):429–433, 1990.
  * [7] E. Agliari, A. Barra, A. Galluzzi, F. Guerra, and F. Moauro.  Multitasking associative networks.  Phys. Rev. Lett., 109:268101, Dec 2012.
  * [8] E. Agliari, A. Annibale, A. Barra, A. C. C Coolen, and D. Tantari.  Immune networks: Multitasking capabilities near saturation.  J. Phys. A Math. Theor., 46(41):415003, 2013.
  * [9] J. J. Hopfield.  Neural Networks and Physical Systems with Emergent Collective Computational Abilities.  Proc. Natl. Acad. Sci., U. S. A, 79:2554–2558, 1982.  Reprinted in [78].
  * [10] J. Hertz, A. Krogh, and R. G. Palmer.  Introduction to the Theory of Neural Computation.  CRC Press, 2018.
  * [11] D. Amit, H. Gutfreund, and H. Sompolinsky.  Statistical Mechanics of Neural Networks Near Saturation.  Annals of Physics, 173:30–67, 1987.
  * [12] P. Sollich, D. Tantari, A. Annibale, and A. Barra.  Extensive parallel processing on scale-free networks.  Phys. Rev. Lett., 113(23):238106, 2014.
  * [13] D. Challet and Y. C. Zhang.  Emergence of Cooperation and Organization in An Evolutionary Game.  Physica A, 246:407–418, 1997.
  * [14] A. C. C. Coolen.  The Mathematical Theory of Minority Games–Statistical Mechanics of Interacting Agents.  Oxford University Press, Oxford, 2005.
  * [15] G. Iori.  Avalanche Dynamics and Trading Friction Effects on Stock Market Returns.  Int. J. Mod. Phys. C, 10:1149–1162, 1999.
  * [16] S. Bornholdt.  Expectation Bubbles in a Spin Model of Market Intermittency from Frustration Across Scales.  Int. J. Mod. Phys. C, 12:667–674, 2001.
  * [17] K. Anand and R. Kühn.  Phase Transitions in Operational Risk.  Phys. Rev. E, 75:016111, 2007.
  * [18] J. P. L. Hatchett and R. Kühn.  Effects of Economic Interactions on Credit Risk.  J. Phys. A, 39:2231–2251, 2006.
  * [19] M. Weigt and A. K Hartmann.  Number of guards needed by a museum: A phase transition in vertex covering of random graphs.  Phys. Rev. Lett., 84(26):6118, 2000.
  * [20] S. Cocco and R.Monasson.  Statistical physics analysis of the computational complexity of solving random satisfiability problems using backtrack algorithms.  Eur. Phys. J. B, 22(4):505–531, 2001.
  * [21] O. C. Martin, R. Monasson, and R. Zecchina.  Statistical mechanics methods and phase transitions in optimization problems.  Theoretical computer science, 265(1-2):3–67, 2001.
  * [22] S. Franz, M. Leone, F. Ricci-Tersenghi, and R. Zecchina.  Exact Solutions for Diluted Spin Glasses and Optimization Problems.  Phys. Rev. Lett., 87:127209, 2001.
  * [23] M. Mézard, G. Parisi, and R. Zecchina.  Analytic and algorithmic solution of random satisfiability problems.  Sci., 297(5582):812–815, 2002.
  * [24] M.Mézard and G.Parisi.  The bethe lattice spin glass revisited.  Eur. Phys. J. B, 20(2):217–233, 2001.
  * [25] A.C.C. Coolen.  Replica methods for loopy sparse random graphs.  In J. Phys. Conf. Ser., volume 699, page 012022, 2016.
  * [26] A. Annibale, A. C. C. Coolen, and G. Bianconi.  Network Resilience Against Intelligent Attacks Constrained by the Degree-Dependent Node Removal Cost.  J. Phys. A, 43:395001, 2010.
  * [27] F. Altarelli, A. Braunstein, L. Dall’Asta, and R. Zecchina.  Optimizing spread dynamics on graphs by message passing.  J. Stat. Mech.: Theory Exp, 2013(09):P09011, 2013.
  * [28] A. Y. Lokhov, M. Mézard, and Zdeborová.  Dynamic message-passing equations for models with unidirectional dynamics.  Phys. Rev. E, 91(1):012811, 2015.
  * [29] P. Paga and R. Kühn.  Contagion in an Interacting Economy.  JSTAT, 03:P03008, 2015.
  * [30] B. Li and D. Saad.  Impact of presymptomatic transmission on epidemic spreading in contact networks: A dynamic message-passing analysis.  Phys. Rev. E, 103(5):052303, 2021.
  * [31] B. Karrer and M. E. J.Newman.  Message Passing Approach for General Epidemic Models.  Phys. Rev. E, 82:016101, 2010.
  * [32] A. Y. Lokhov, M. Mézard, H. Ohta, and L. Zdeborová.  Inferring the origin of an epidemic with a dynamic message-passing algorithm.  Phys. Rev. E, 90(1):012801, 2014.
  * [33] Y. Roudi and J. Hertz.  Dynamical tap equations for non-equilibrium ising spin glasses.  J. Stat. Mech., 2011(03):P03031, 2011.
  * [34] P. Zhang.  Inference of kinetic Ising model on sparse graphs.  J. Stat. Phys, 148(3):502–512, 2012.
  * [35] E Aurell and H. Mahmoudi.  Dynamic mean-field and cavity methods for diluted ising systems.  Phys. Rev. E, 85(3):031119, 2012.
  * [36] C. De Dominicis.  Dynamics as Substitute for Replicas in Systems with Quenched Random Impurities.  Phys. Rev. B, 18:4913–4919, 1978.
  * [37] J. P. L. Hatchett, B. Wemmenhove, I. Pérez Castillo, T. Nikoletopoulos, N. S. Skantzos, and A. C. C. Coolen.  Parallel Dynamics of Disordered Ising Spin Systems on Finitely Connected Random Graphs.  J. Phys. A, 37:6201–6220, 2004.
  * [38] K. Mimura and A. C. C. Coolen.  Parallel dynamics of disordered ising spin systems on finitely connected directed random graphs with arbitrary degree distributions.  J. Phys. A, 42(41):415001, 2009.
  * [39] I. Neri and D. Bollé.  The cavity approach to parallel dynamics of ising spins on a graph.  J. Stat. Mech.: Theory Exp, 2009(08):P08009, 2009.
  * [40] R. Kühn and T. Rogers.  Heterogeneous Micro-Structure of Percolation in Sparse Networks.  Europhys. Lett., 118:68003, 2017.
  * [41] J. Leskovec and A. Krevl.  SNAP Datasets: Stanford Large Network Dataset Collection.  http://snap.stanford.edu/data, 2014.
  * [42] R. Albert and A.-L. Barabási.  Statistical Mechanics of Complex Networks.  Rev. Mod. Phys., 74:47–97, 2002.
  * [43] S. N. Dorogovtsev and J. F. F Mendes.  Evolution of Networks: from Biological Networks to the Internet and WWW.  Oxford University Press, Oxford, 2003.
  * [44] M. E. J. Newman.  Networks: an Introduction, 2nd Ed. Oxford Univ. Press, Oxford, 2018.
  * [45] E. Estrada.  The Structure of Complex Networks: Theory and Applications.  Oxford University Press, Oxford, 2011.
  * [46] V. Latora, V. Nicosia, and G. Russo.  Complex Networks: Principles, Methods and Applications.  Cambridge University Press, Cambridge, 2017.
  * [47] G. Torrisi, A. Annibale, and R. Kühn.  Overcoming the complexity barrier of the dynamic message-passing method in networks with fat-tailed degree distributions.  Phys. Rev. E, 104(4):045313, 2021.
  * [48] T. Fink and R. Hannam.  Boolean composition restricts biological logics.  arXiv preprint arXiv:2109.12551, 2021.
  * [49] A. J Gates, R. B Correia, X. Wang, and L. M Rocha.  The effective graph reveals redundancy, canalization, and control pathways in biochemical regulation and signaling.  Proc. Natl. Acad. Sci., 118(12), 2021.
  * [50] G. Torrisi, R. Kühn, and A. Annibale.  Percolation in the Gene Regulatory Network.  JSTAT, 083501:31p, 2020.
  * [51] R. Hannam, R. Kühn, and A. Annibale.  Percolation in Bipartite Boolean Networks and its Role in Sustaining Life.  J. Phys. A, 52:334002, 2019.
  * [52] JA Hertz, G Grinstein, and SA Solla.  Memory networks with asymmetric bonds.  In AIP Conference Proceedings, volume 151, pages 212–218. American Institute of Physics, 1986.
  * [53] H Gutfreund, JD Reger, and AP Young.  The nature of attractors in an asymmetric spin glass with deterministic dynamics.  Journal of Physics A: Mathematical and General, 21(12):2775, 1988\.
  * [54] H. Han, J. Cho, S. Lee, A. Yun, H. Kim, et al.  Trrust v2: an expanded reference database of human and mouse transcriptional regulatory interactions.  Nucleic Acids Res., 46(D1):D380–D386, 2018.
  * [55] G. Del Ferraro and E. Aurell.  Dynamic message-passing approach for kinetic spin models with reversible dynamics.  Phys. Rev. E, 92(1):010102, 2015.
  * [56] M. Mézard and J. Sakellariou.  Exact mean-field inference in asymmetric kinetic ising systems.  J. Stat. Mech., 2011(07):L07001, 2011.
  * [57] J. P Gleeson.  High-accuracy approximation of binary-state dynamics on networks.  Phys Rev. Let., 107(6):068701, 2011.
  * [58] B. Deplancke and N. Gheldof.  Gene regulatory networks: methods and protocols.  Springer, 2012.
  * [59] A. F. Bitbol, R. S Dwyer, L. J Colwell, and N. Wingreen.  Inferring interaction partners from protein sequences.  Proc. Natl. Acad. Sci. U.S.A., 113(43):12180–12185, 2016.
  * [60] F. Morcos, A. Pagnani, B. Lunt, A. Bertolino, D. S Marks, C. Sander, R. Zecchina, J. N Onuchic, T. Hwa, and M. Weigt.  Direct-coupling analysis of residue coevolution captures native contacts across many protein families.  Proc. Natl. Acad. Sci. U.S.A., 108(49):E1293–E1301, 2011.
  * [61] A. Cavagna, I. Giardina, F. Ginelli, T. Mora, D. Piovani, R. Tavarone, and A. M Walczak.  Dynamical maximum entropy approach to flocking.  Phys. Rev.E, 89(4):042707, 2014.
  * [62] George T Cantwell and Mark EJ Newman.  Message passing on networks with loops.  Proceedings of the National Academy of Sciences, 116(47):23398–23403, 2019.
  * [63] A. Kirkley, G. T Cantwell, and MEJ Newman.  Belief propagation for networks with loops.  Science Advances, 7(17):eabf1211, 2021.
  * [64] R. Lambiotte, M. Rosvall, and I. Scholtes.  From networks to optimal higher-order models of complex systems.  Nat. Phys., 15(4):313–320, 2019.
  * [65] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri.  Networks beyond pairwise interactions: Structure and dynamics.  Physics Reports, 2020.
  * [66] G. Ghoshal, V. Zlatić, G. Caldarelli, and M. EJ Newman.  Random hypergraphs and their applications.  Physical Review E, 79(6):066118, 2009.
  * [67] G. Petri and A. Barrat.  Simplicial activity driven model.  Physical review letters, 121(22):228301, 2018.
  * [68] T. R. Kirkpatrick and D. Thirumalai.  Dynamics of the Structural Glass Transition and the $p$-Spin Interaction Spin-Glass Model.  Phys. Rev. Lett., 58:2091–2094, 1987.
  * [69] Thomas R.  On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations.  In Numerical methods in the study of critical phenomena, pages 180–193. Springer, 1981.
  * [70] E. Plahte, T. Mestl, and S. W Omholt.  Feedback loops, stability and multistationarity in dynamical systems.  J. Biol. Syst, 3(02):409–413, 1995.
  * [71] M. Freeman.  Feedback control of intercellular signalling in development.  Nature, 408(6810):313–319, 2000.
  * [72] H. Zhang, Y. Chen, and Y. Chen.  Noise propagation in gene regulation networks involving interlinked positive and negative feedback loops.  PloS one, 7(12):e51840, 2012.
  * [73] S. Chakravarty and A. Csikász-Nagy.  Systematic analysis of noise reduction properties of coupled and isolated feed-forward loops.  PLoS Comput. Biol., 17(12):e1009622, 2021.
  * [74] A. M. Mambuca, C. Cammarota, and I. Neri.  Dynamical systems on large networks with predator-prey interactions are stable and exhibit oscillations.  Phys. Rev. E, 105:014305, Jan 2022.
  * [75] E. Aurell and H. Mahmoudi.  A message-passing scheme for non-equilibrium stationary states.  J. Stat. Mech., 2011(04):P04014, 2011.
  * [76] L. Cardelli, A. Csikász-Nagy, N. Dalchau, M. Tribastone, and M. Tschaikowski.  Noise reduction in complex biological switches.  Sci. Rep., 6(1):1–12, 2016.
  * [77] C. Campajola, F. Lillo, P. Mazzarisi, and D. Tantari.  On the equivalence between the kinetic ising model and discrete autoregressive processes.  Journal of Statistical Mechanics: Theory and Experiment, 2021(3):033412, 2021.
  * [78] J. A. Anderson and E. Rosenfeld, editors.  Neurocomputing: Foundations of Research.  MIT Press, Cambridge, 1988.

##  Appendix A High noise regime

In the high noise ($\beta\rightarrow 0$ limit), the expansion of
$\Phi_{\beta}(x)$ in Equation (9) to first order in $\beta x$ 777Note that
$\Phi_{\beta}(x)=\Phi_{1}(\beta x)$, both for thermal noise, and for Gaussian
noise, if we equate $\beta=\sigma$ in the Gaussian case gives

$P_{i}(t+1)\approx\frac{1}{2}+\beta\Phi_{1}^{\prime}(0)\left(\sum_{j}J_{ij}P_{j}(t)-\vartheta_{i}\right)\,.$
(53)

Using $P_{j}(t)=\frac{1}{2}+\mathcal{O}(\beta)$, we obtain that

$P_{i}(t+1)=\frac{1}{2}+\beta\Phi_{1}^{\prime}(0)\left(\frac{1}{2}\sum_{j}^{N}J_{ij}-\vartheta_{i}\right)+\mathcal{O}(\beta^{2})\,,$
(54)

which is independent of the time-step $t$. The result of Equation (54) states
that in the large noise limit, the node activation of node $i$, $P_{i}$, is
solely determined by the imbalance of the interactions with $i$’s
predecessors, i.e. by $\sum_{j}J_{ij}$. The mean node activation probability
over sites $\langle P\rangle=\sum_{i}P_{i}/N$ reduces to

$\langle
P\rangle=\frac{1}{2}+\beta\Phi_{1}^{\prime}(0)\left(\frac{1}{2N}\sum_{i,j}J_{ij}-\frac{1}{N}\sum_{i}\vartheta_{i}\right)+\mathcal{O}(\beta^{2})\,,$
(55)

where the $\mathcal{O}(\beta)$ correction depends only on the macroscopic
quantities $\sum_{i,j}J_{ij}/N$ and $\sum_{i}\vartheta_{i}/N$. This result is
compatible with the expansion around the paramagnetic phase done in ref. [39].

Similarly, the high noise approximation for the AND gate dynamics discussed in
Equation (31) gives

$P_{i}\approx\frac{1}{2}+\beta\Phi_{1}^{\prime}(0)\left[\sum_{\mu}J_{i\mu}\mathcal{P}_{\mu}-\vartheta_{i}\right]\,,$
(56)

where

$\mathcal{P}_{\mu}=\left(\frac{1}{2}\right)^{c_{\mu}}+\beta\Phi_{1}^{\prime}(0)\sum_{j\in\partial_{\mu}}\sum_{\nu}J_{j\nu}\left(\frac{1}{2}\right)^{c_{\nu}}+\mathcal{O}(\beta^{2})\,.$
(57)

In this case, the average node activation probability is therefore given by

$\langle
P\rangle=\frac{1}{2}+\beta\Phi_{1}^{\prime}(0)\left[\frac{1}{N}\sum_{i,\mu}J_{i\mu}\left(\frac{1}{2}\right)^{c_{\mu}}-\frac{1}{N}\sum_{i}\vartheta_{i}\right]+\mathcal{O}(\beta^{2})\
.$ (58)

In the high noise limit, the average activation probability of factor nodes
converges to

$\langle\mathcal{P}\rangle=\sum_{c_{\mu}}P(c_{\mu})\left(\frac{1}{2}\right)^{c_{\mu}}\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=G^{\mathrm{C}}_{0}(1/2)\
,\quad\mbox{as }\beta\to 0.$ (59)

with $G^{\mathrm{C}}_{0}(x)$ the generating function of the in-degree of
factor nodes. This demonstrates, among other things, that the signalling
through AND gates with large in-degrees is effectively suppressed in the high
noise limit. Thus, factor nodes with large in-degree may operate as noise
filters.

##  Appendix B Derivation of the One-Time Approximation

#### Dynamic cavity trajectory with bi-directional links

For completeness, we include in this appendix a brief derivation of the OTA
approximation from the dynamic cavity equations, mainly following Ref. [39,
35]. We also compare different closure schemes of the OTA which have been
proposed in the literature and propose a new version as an alternative.

Let us define $P_{i}\left(\bm{n}_{i}^{0,\dots,t}|\vartheta_{i}\right)$ the
probability of the trajectory
$\bm{n}_{i}^{0,\dots,t}\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=(n_{i}(0)\dots
n_{i}(t))$ of node $i$ over the time span $0,\dots,t$, and
$P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t-1}|\bm{\vartheta}_{j}^{(i)0,\dots,t-2}\right)$
the conditional probability of the trajectory
$\bm{n}_{j}^{0,\dots,t-1}\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=(n_{j}(0)\dots
n_{j}(t-1))$ of node $j$, given the history of node $i$, in the cavity graph
where node $i$ and its links have been removed, so that the past values of
node $i$ act as external thresholds. Their contribution can be combined with
those of the constant threshold $\vartheta_{j}$, to define time-dependent
thresholds
$\bm{\vartheta}_{j}^{(i)0,\dots,t-1}\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=(\vartheta_{j}^{(i)}(0)\dots\vartheta_{j}^{(i)}(t-1))$
with

$\vartheta_{j}^{(i)}(s)=\vartheta_{j}-J_{ji}n_{i}(s)\,,$ (60)

where $J_{ji}n_{i}^{s}$ is the signal received by $j$ from $i$ as a retarded
interaction. Above, we have used the notation with time superscript to
indicate the collection of variables over the trajectory. In what follows, to
ease the notation, we will also use superscripts for single time quantities,
e.g. $n_{i}(s)=n_{i}^{s}$ and $\vartheta_{j}^{(i)}(s)=\vartheta_{j}^{(i),s}$.

The probability of a single-node trajectory does not factorise in the time-
steps, but one needs to solve the set of equations, see Ref. [39, 75]

$P_{i}\left(\bm{n}_{i}^{0,\dots,t}|\vartheta_{i}\right)=\sum_{\bm{n}^{0}_{\partial_{i}}}\dots\sum_{\bm{n}^{t-1}_{\partial_{i}}}\left[\prod_{s=0}^{t-1}W\left[n_{i}^{s+1}|h_{i}(\bm{n}_{\partial_{i}}^{s}),\vartheta_{i}\right]\right]\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\prod_{j\in\partial_{i}}P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t-1}|\bm{\vartheta}_{j}-J_{ji}\bm{n}_{i}^{0,\dots,t-2}\right)P_{i}(n_{i}^{0})\,,$
(61)
$P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t-1}|\bm{\vartheta}^{(i)0,\dots,t-2}_{j}\right)=\sum_{\bm{n}^{0}_{\partial_{j}\setminus
i}}\dots\sum_{\bm{n}^{t-2}_{\partial_{j}\setminus
i}}\left[\prod_{s=0}^{t-2}W\left[n_{j}^{s+1}|h_{j}^{(i)}(\bm{n}_{\partial_{j}}^{s}),\vartheta_{j}^{(i),s}\right]\right]\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\prod_{\ell\in\partial_{j}\setminus
i}P_{\ell}^{(j)}\left(\bm{n}_{\ell}^{0,\dots,t-2}|\bm{\vartheta}_{\ell}-J_{\ell
j}\bm{n}_{j}^{0,\dots,t-3}\right)P_{j}(n_{j}^{0})\,,$ (62)

where $\bm{\vartheta}_{j}=(\vartheta_{j},\dots,\vartheta_{j})$ is a constant
vector with identical entries $\vartheta_{j}$. Note that for any node
$j\in\partial_{i}$, i.e. such that $J_{ij}\neq 0$, one has $J_{ji}\neq 0$ only
if the link $(i,j)$ is bidirectional, so the retarded self-interaction terms
in Equation (60) are absent for unidirectional links. The transition
probability
$W\left[n_{j}^{s+1}|h_{j}^{(i)}(\bm{n}_{\partial_{j}}^{s}),\vartheta_{j}^{(i),s}\right]$
indicates the conditional probability of the state $n_{j}^{s+1}$ at time $s+1$
given the values of the cavity local field
$h_{j}^{(i)}(\bm{n}_{\partial_{j}}^{s})$,

$h_{j}^{(i)}(\bm{n}_{\partial_{j}}^{s})\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\sum_{\ell\in\partial_{j}\setminus{i}}J_{j\ell}n_{\ell}^{s}\,,$
(63)

and the cavity threshold $\vartheta_{j}^{(i),s}$. The detailed function form
of
$W\left[n_{j}^{s+1}|h_{j}^{(i)}(\bm{n}_{\partial_{j}}^{s}),\vartheta_{j}^{(i),s}\right]$
depends on the choice of the noise distribution. In analogy with the above we
adopt a thermal noise model for which

$W\left[n_{j}^{s+1}|h_{j}^{(i)}(\bm{n}_{\partial_{j}}^{s}),\vartheta_{j}^{(i),s}\right]=\frac{1}{2}\left\\{1+(2n_{j}^{s+1}-1)\tanh\left[\frac{\beta}{2}\left(h_{j}^{(i)}(\bm{n}_{\partial_{j}}^{s})-\vartheta_{j}^{(i),s}\right)\right]\right\\}\,.$
(64)

In Equation (61) and Equation (62) the probability of the state of node $i$
depends on the states of neighbours $\partial_{i}$, which in turn depend
parametrically on the node $i$ through the cavity thresholds. This feedback
effect is responsible for the non-Markovian structure of Equation (61) and
Equation (62), entailing that the computational complexity for a trajectory of
time length $t$ grows as $2^{t}$, making the evaluation feasible only for
relatively few time-steps.

From now on, we focus on the single site cavity probability and we discuss the
quantity
$P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t}|\bm{\vartheta}^{(i),0,\dots,t-1}_{j}\right)$
for the trajectory up to time $t$ and not up to time $t-1$ as done in Equation
(62) to simplify the notation of the lagged time dependence. A pragmatic
approach is to assume that the single site cavity probability trajectory
$P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t}|\bm{\vartheta}^{(i),0,\dots,t-1}_{j}\right)$
factorises in a Markovian fashion under the one-time approximation (OTA)[39]

$P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t}|\bm{\vartheta}^{(i),0,\dots,t-1}_{j}\right)\approx
P^{(i)}_{j}(n_{j}^{0})\prod_{s=0}^{t-1}P^{(i)}_{j}(n_{j}^{s+1}|\vartheta^{(i),s}_{j})\,.$
(65)

Substituting Equation (65) into Equation (62)

$P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t}|\bm{\vartheta}^{(i),0,\dots,t-1}_{j}\right)=\left[\sum_{\bm{n}^{t-1}_{\partial_{j}\setminus
i}}W\left[n_{j}^{t}|h_{j}^{(i)}(\bm{n}_{\partial_{j}}^{t-1}),\vartheta_{j}^{(i),t-1}\right]\prod_{\ell\in\partial_{j}\setminus
i}P_{\ell}^{(j)}\left(n_{\ell}^{t-1}|\vartheta_{\ell}-J_{\ell
j}n_{j}^{t-2}\right)\right]\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times
P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t-1}|\bm{\vartheta}^{(i)0,\dots,t-2}_{j}\right)\,,$
(66)

where we have grouped together the terms of Equation (65) at the previous
time-steps

$P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t-1}|\bm{\vartheta}^{(i)0,\dots,t-2}_{j}\right)=\sum_{\bm{n}^{0}_{\partial_{j}\setminus
i}}\dots\sum_{\bm{n}^{t-2}_{\partial_{j}\setminus
i}}\left[\prod_{s=0}^{t-2}W\left[n_{j}^{s+1}|h_{j}^{(i)}(\bm{n}_{\partial_{j}}^{s}),\vartheta_{j}^{(i),s}\right]\right]\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\prod_{\ell\in\partial_{j}\setminus
i}\left[\prod_{s=0}^{t-3}P_{\ell}^{(j)}\left(n_{\ell}^{s+1}|\vartheta_{\ell}-J_{\ell
j}n_{j}^{s}\right)P_{\ell}(n_{\ell}^{0})\right]P_{j}(n_{j}^{0})\,.$ (67)

Applying the OTA factorization to the rightmost factor in Equation (66), i.e.
Equation (67), one obtains

$P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t}|\bm{\vartheta}^{(i),0,\dots,t-1}_{j}\right)=\left[\sum_{\bm{n}^{t-1}_{\partial_{j}\setminus
i}}W\left[n_{j}^{t}|h_{j}^{(i)}(\bm{n}_{\partial_{j}}^{t-1}),\vartheta_{j}^{(i),t-1}\right]\prod_{\ell\in\partial_{j}\setminus
i}P_{\ell}^{(j)}\left(n_{\ell}^{t-1}|\vartheta_{\ell}-J_{\ell
j}n_{j}^{t-2}\right)\right]\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\left[\prod_{s=0}^{t-2}P^{(i)}_{j}(n_{j}^{s+1}|\vartheta^{(i),s}_{j})\right]P_{j}(n_{j}^{0})\,.$
(68)

This expression is particularly suitable to compute the marginal over the
trajectory of node $j$, i.e.,

$P_{j}^{(i)}\left(n_{j}^{t}|\bm{\vartheta}_{j}^{(i),0,\dots,t-1}\right)=\sum_{\bm{n}_{j}^{0\dots
t-1}}P_{j}^{(i)}\left(\bm{n}_{j}^{0,\dots,t}|\bm{\vartheta}^{(i),0,\dots,t-1}_{j}\right)\,.$
(69)

Thanks to the factorization of the last term in Equation (68), the sum over
$\bm{n}_{j}^{0\dots t-3,t-1}$ is trivial, as the term in the first square
brackets does not depend on these variables, and conditional probabilities are
normalised to one.

$P_{j}^{(i)}\left(n_{j}^{t}|\vartheta_{j}^{(i),t-1},\vartheta_{j}^{(i),t-3}\right)=\\\
\sum_{n_{j}^{t-2}}\left[\sum_{\bm{n}^{t-1}_{\partial_{j}\setminus
i}}W\left[n_{j}^{t}|h_{j}^{(i)}(\bm{n}_{\partial_{j}}^{t-1}),\vartheta_{j}^{(i),t-1}\right]\prod_{\ell\in\partial_{j}\setminus
i}P_{\ell}^{(j)}\left(n_{\ell}^{t-1}|\vartheta_{\ell}-J_{\ell
j}n_{j}^{t-2}\right)\right]P^{(i)}_{j}(n_{j}^{t-2}|\vartheta^{(i),t-3}_{j})$
(70)

In Equation (70) the probability of node $j$ depends on the two cavity
thresholds $\vartheta_{j}^{(i),t-1},\vartheta_{j}^{(i),t-3}$. This expression
is directly derived from the OTA factorisation Equation (65). The one-step
object $P_{j}^{(i)}\left(n_{j}^{t}|\vartheta_{j}^{(i),t-1}\right)$ can be
obtained imposing a “closure” condition that enforces the Markovian behaviour,
and different choices could be made. In [39, 75] the authors assume that, in
the stationary state, the cavity threshold does not have an explicit time
dependence $\vartheta_{j}^{(i),s}\approx\vartheta_{j}^{(i)}$ $\forall s$. In
[34] instead, the cavity term is approximated by the non-cavity marginal
$P^{(i)}_{j}(n_{j}^{t-2}|\vartheta^{(i),t-3}_{j})\approx
P_{j}(n_{j}^{t-2}|\vartheta_{j})$. Another approach is to truncate the
retarded dependence by taking an average over the state $n_{i}^{t-3}$, which
we detail below. Remembering from Equation (60) that the
$\vartheta_{j}^{(i),s}=\vartheta_{j}-J_{ji}n_{i}^{s}$, then

$P(\vartheta_{j}^{(i),t-3})=\sum_{n_{i}^{t-3}}\delta_{\vartheta_{j}^{(i),t-3},\vartheta_{j}-J_{ji}n_{i}^{t-3}}P(n_{i}^{t-3}|\vartheta_{i})\,,$
(71)

with $\delta$ indicating the Kronecker delta. The one-time cavity marginal
appearing in Equation (65) is given by

$P_{j}^{(i)}\left(n_{j}^{t}|\vartheta_{j}^{(i),t-1}\right)=\sum_{n_{i}^{t-3}}P_{j}^{(i)}\left(n_{j}^{t}|\vartheta_{j}^{(i),t-1},\vartheta_{j}-J_{ji}n_{i}^{t-3}\right)P_{i}(n_{i}^{t-3}|\vartheta_{i})\,.$
(72)

Substituting Equation (70) in Equation (72), the expression for the one-step
probability becomes

$P_{j}^{(i)}\left(n_{j}^{t}|\vartheta_{j}^{(i),t-1}\right)=\sum_{n_{j}^{t-2}}\left[\sum_{\bm{n}^{t-1}_{\partial_{j}\setminus
i}}W\left[n_{j}^{t}|h_{j}^{(i)}(\bm{n}_{\partial_{j}}^{t-1}),\vartheta_{j}^{(i),t-1}\right]\prod_{\ell\in\partial_{j}\setminus
i}P_{\ell}^{(j)}\left(n_{\ell}^{t-1}|\vartheta_{\ell}-J_{\ell
j}n_{j}^{t-2}\right)\right]P^{(i)}_{j}(n_{j}^{t-2}|\vartheta_{j})$ (73)

with

$P^{(i)}_{j}(n_{j}^{t-2}|\vartheta_{j})=\sum_{n_{i}^{t-3}}P^{(i)}_{j}(n_{j}^{t-2}|\vartheta_{j}-J_{ji}n_{i}^{t-3})P_{i}(n_{i}^{t-3}|\vartheta_{i})$
(74)

We compare the different OTA closure conditions in Equation (70), which we
denote:

  * •

C.1 corresponding to the assumption
$P^{(i)}_{j}(n_{j}^{t-2}|\vartheta^{(i),t-3}_{j})\approx
P^{(i)}_{j}(n_{j}^{t-2}|\vartheta^{(i),t-1}_{j})$, from Ref. [35],

  * •

C.2 corresponding to the assumption
$P^{(i)}_{j}(n_{j}^{t-2}|\vartheta^{(i),t-3}_{j})\approx
P_{j}(n_{j}^{t-2}|\vartheta_{j})$, from Ref. [34],

  * •

C.3 corresponding to the assumption
$P^{(i)}_{j}(n_{j}^{t-2}|\vartheta^{(i),t-3}_{j})\approx
P^{(i)}_{j}(n_{j}^{t-2}|\vartheta_{j})$ using the definition in Equation (74).

We evaluate the distribution of the stationary single-node activation
probability for the OTA expressions corresponding to closures C.1, C.2, and
C.3. We benchmark the theory with Monte Carlo simulations for symmetric and
antisymmetric networks, and we show in Fig. 13 the distribution of the
distance between theory and simulation for the single-site activation
probabilities. Our results indicate that in the symmetric case closure C.1
outperforms both C.2 and C.3, which is expected since the expression
associated with C.1 admits the equilibrium solution of belief propagation, as
discussed in Ref. [39]. In the antisymmetric case instead, the relative
performances of the methods are changed, with closure C.1 providing the worst
result roughly by a factor of two compared to closures C.2 and C.3. We use the
mean square distance between OTA closures and simulation to quantify the error
in Tab. 1. Our results indicate that closures C.2 and C.3 are less dependent
on the network symmetry, while C.1 presents important variations between the
two cases examined. From a computational point of view, closures C.1 and C.2
are comparable and they are simpler to perform compared to C.3, since in
closure C.3, the Eq. (74) needs to be performed for every link and for every
time-step of interest. Hence we adopt closure C.2 since it is simpler to
implement than C.3, and our results suggest it gives the same performances as
closure C.3.

Figure 13: Comparison of different closures in the OTA. C.1 is the closure of Ref. [35]. C.2 is the closure of Ref. [34], C.3 is the closure of Equation (74). The distribution of the difference between the probability obtained using OTA and the probability obtained using Monte Carlo simulations. Network with symmetric (left) and antisymmetric interactions (right). The network belongs to the random regular graph class with $k^{\mathrm{in}}=k^{\mathrm{out}}=3$. Simulations are obtained through an average of over 500 trajectories. Parameters are $N=10000$, $T/J=1$, and $\vartheta=0$, $t_{s}=10^{4}$ steps. Table 1: Error, as defined in Equation (46), of different closure schemes for symmetric, antisymmetric and uncorrelated interactions. Same parameters as in Fig. 13. | C.1 | C.2 | C.3
---|---|---|---
symmetric | $2.3\text{\times}{10}^{-04}$ | $6.0\text{\times}{10}^{-04}$ | $6.2\text{\times}{10}^{-04}$
antisymmetric | $8.6\text{\times}{10}^{-04}$ | $5.5\text{\times}{10}^{-04}$ | $5.3\text{\times}{10}^{-04}$
uncorrelated | $3.8\text{\times}{10}^{-04}$ | $3.4\text{\times}{10}^{-04}$ | $3.4\text{\times}{10}^{-04}$