# A unifying framework for mean-field theories of asymmetric kinetic Ising
systems

Miguel Aguilera<EMAIL_ADDRESS>IAS-Research Center for Life, Mind, and
Society, Department of Logic and Philosophy of Science, University of the
Basque Country, Donostia, Spain.  Department of Informatics & Sussex
Neuroscience, University of Sussex, Brighton, UK.  ISAAC Lab, Aragón Institute
of Engineering Research, University of Zaragoza, Zaragoza, Spain.  S. Amin
Moosavi  Graduate School of Informatics, Kyoto University, Kyoto, Japan.
Present address: Department of Neuroscience, Brown University, USA  Hideaki
Shimazaki  Graduate School of Informatics, Kyoto University, Kyoto, Japan.
Present address: Center for Human Nature, Artificial Intelligence, and
Neuroscience (CHAIN), Hokkaido University, Japan

###### Abstract

Kinetic Ising models are powerful tools for studying the non-equilibrium
dynamics of complex systems. As their behavior is not tractable for large
networks, many mean-field methods have been proposed for their analysis, each
based on unique assumptions about the system’s temporal evolution. This
disparity of approaches makes it challenging to systematically advance mean-
field methods beyond previous contributions. Here, we propose a unifying
framework for mean-field theories of asymmetric kinetic Ising systems from an
information geometry perspective. The framework is built on Plefka expansions
of a system around a simplified model obtained by an orthogonal projection to
a sub-manifold of tractable probability distributions. This view not only
unifies previous methods but also allows us to develop novel methods that, in
contrast with traditional approaches, preserve the system’s correlations. We
show that these new methods can outperform previous ones in predicting and
assessing network properties near maximally fluctuating regimes.

††preprint: APS/123-QED

## Introduction

Advances in high-throughput data acquisition technologies for very large
biological and social systems are providing unprecedented possibilities to
investigate their complex, non-equilibrium dynamics. For example, optical
recordings from genetically modified neural populations make it possible to
simultaneously monitor activities of the whole neural network of behaving C.
elegans [1] and zebrafish [2], as well as thousands of neurons in the mouse
visual cortex [3]. Such networks generally exhibit out-of-equilibrium dynamics
[4], and are often found to self-organize near critical regimes at which their
fluctuations are maximized [5, 6]. Evolution of such systems cannot be
faithfully captured by methods assuming an asymptotic equilibrium state.
Therefore, in general, there is a pressing demand for mathematical tools to
study the dynamics of large-scale, non-equilibrium complex systems and to
analyze high-dimensional data sets recorded from them.

The kinetic Ising model with asymmetric couplings is a prototypical model for
studying such non-equilibrium dynamics in biological [7, 8] and social systems
[9]. It is described as a discrete-time Markov chain of interacting binary
units, resembling the nonlinear dynamics of recurrently connected neurons. The
model exhibits non-equilibrium behavior when couplings are asymmetric or when
model parameters are subject to rapid changes, ruling out quasi-static
processes. These conditions induce a time reversal asymmetry in dynamical
trajectories, leading to positive entropy production (the second law of
thermodynamics) as revealed by the fluctuation theorem [10, 11, 12, 13, 14,
15] (see [16, 17] for reviews). This time-asymmetry is characteristic of non-
equilibrium systems as it can only be displayed by systems in which energy
dissipation takes place [18]. In the case of symmetric connections and static
parameters, the model converges to an equilibrium stationary state. Therefore,
it is a generalization of its equilibrium counterpart known as the
(equilibrium) Ising model [19].

The forward Ising problem refers to calculating statistical properties of the
model, such as mean activation rates (mean magnetizations of spins) and
correlations, given the parameters of the model. In contrast, inference of the
model parameters from data is called the inverse Ising problem [20]. In this
regard, kinetic Ising models [21, 22] and their equilibrium counterparts [23,
24, 25] have become popular tools for modeling and analyzing biological and
social systems. In addition, they capture memory retrieval dynamics in
classical associative networks. Namely, they are equivalent to the Boltzmann
machine, extensively used in machine learning applications [20].
Unfortunately, exact solutions of the forward and inverse problems often
become computationally too expensive due to the combinatorial explosion of
possible patterns in large networks or the high volume of data, and
applications of exact or sampling-based methods are limited in practice to
around a hundred of neurons [25, 26, 5]. Therefore, analytical approximation
methods are necessary for analysing large systems. In this endeavour, mean-
field methods have emerged as powerful tools to track down otherwise
intractable statistical quantities.

The standard mean-field approximations to study equilibrium Ising models are
the classical naive mean-field (nMF) and the more accurate Thouless-Anderson-
Palmer (TAP) approximations [27]. These methods have also been employed to
solve the inverse Ising problem [28, 29, 30, 31]. Plefka demonstrated that the
nMF and TAP approximations for the equilibrium model can be derived using the
power series expansion of the free energy around a model of independent spins,
a method which is now referred to as the Plefka expansion [32]. This expansion
up to the first and second orders leads to the nMF and TAP mean-field
approximations respectively. The Plefka expansion was later formalized by
Tanaka and others in the framework of information geometry [33, 34, 35, 36,
37].

In non-equilibrium networks, however, the free energy is not directly defined,
and therefore it is not obvious how to apply the Plefka expansions. Kappen and
Spanjers [38] proposed an information geometric approach to mean-field
solutions of the asymmetric Ising model with asynchronous dynamics. They
showed that their second-order approximation for an asymmetric model in the
stationary state is equivalent to the TAP approximation for equilibrium
models. Later, Roudi and Hertz derived TAP equations for non-stationary states
using a Legendre transformation of the generating functional of the set of
trajectories of the model [39]. Another study by Roudi and Hertz extended
mean-field equations to provide expressions for the non-stationary delayed
correlations assuming the presence of equal-time correlations at the previous
step [40]. Yet another interesting method proposed by Mézard and Sakellariou
approximates the local fields by a Gaussian distribution according to the
central limit theorem, yielding more accurate results for fully asymmetric
networks [41]. This method was later extended to include correlations at the
previous time-step, improving the results for symmetric couplings [42]. More
recently, Bachschmid-Romano et al. extended the path-integral methods in [39]
with Gaussian effective fields [43], not only recovering [41] for fully
asymmetric networks but also proposing a method that better approximates mean
rate dynamics by conserving autocorrelations of units. Although many choices
of mean-field methods are available, the diversity of methods and assumptions
makes it challenging to advance systematically over previous contributions.

Here, we propose a unified approach for mean-field approximations of the Ising
model. While our method is applicable to symmetric and equilibrium models, we
focus for generality on asymmetric kinetic Ising models. Our approach is
defined as a family of Plefka expansions in an information geometric space.
This approach allows us to unify and relate existing mean-field methods and
provide expressions for other statistics of the systems such as pairwise
correlations. Furthermore, our approach can be extended beyond classical mean-
field assumptions to propose novel approximations. Here, we introduce an
approximation based on a pairwise model that better captures network
correlations, and we show that it outperforms existing approximations of
kinetic Ising models near a point of maximum fluctuations. We also provide a
data-driven method to reconstruct and test if a system is near a phase
transition by combining the forward and inverse Ising problems, and
demonstrate that the proposed pairwise model more accurately estimates the
system’s fluctuations and its sensitivity to parameter changes. These results
confirm that our unified framework is a useful tool to develop methods to
analyze large-scale, non-equilibrium biological and social dynamics operating
near critical regimes. In addition, since the methods are directly applicable
to Boltzmann machine learning, the geometrical framework introduced here is
relevant in machine learning applications.

The paper is organized as follows. First, we introduce the kinetic Ising model
and its statistical properties of interest. Second, we introduce our framework
for the Plefka approximation methods from a geometric perspective. To explain
how it works, we derive the classical naive and TAP mean-field approximations
under the proposed framework. Third, we show that our approach can unify other
known mean-field approximation methods. We then propose a novel pairwise
approximation under this framework. Finally, we compare different mean-field
approximations in solving the forward and inverse Ising problems, as well as
in performing the data-driven assessment of the system’s sensitivity. The last
section is devoted to discussion.

## Results

### The kinetic Ising model

The kinetic Ising model is the least structured statistical model containing
delayed pairwise interactions between its binary components (i.e., a maximum
caliber model [44]). The system consists of $N$ interacting binary variables
(down or up of Ising spins or inactive or active of neural units)
$s_{i,t}\in\\{-1,+1\\},i=1,2,..,N$, evolving in discrete time steps $t$ with
parallel dynamics. Given the configuration of spins at $t-1$,
$\mathbf{s}_{t-1}=\\{s_{1,t-1},s_{2,t-1},\dots,s_{N,t-1}\\}$, spins
$\mathbf{s}_{t}$ at time $t$ are conditionally independent random variables,
updated as a discrete-time Markov chain, following

$\displaystyle P\left(\mathbf{s}_{t}|\mathbf{s}_{t-1}\right)=$
$\displaystyle\prod_{i}\frac{\mathrm{e}^{s_{i,t}h_{i,t}}}{2\cosh h_{i,t}},$
(1) $\displaystyle h_{i,t}=$ $\displaystyle H_{i}+\sum_{j}J_{ij}s_{j,t-1}.$
(2)

The parameters $\mathbf{H}=\\{H_{i}\\}$ and $\mathbf{J}=\\{J_{ij}\\}$
represent local external fields at each spin and couplings between pairs of
spins respectively. When the couplings are asymmetric (i.e, $J_{ij}\neq
J_{ji}$), the system is away from equilibrium because the process is
irreversible with respect to time. Given the probability mass function of the
previous state $P(\mathbf{s}_{t-1})$, the distribution of the current state
is:

$P\left(\mathbf{s}_{t}\right)=\sum_{\mathbf{s}_{t-1}}P\left(\mathbf{s}_{t}|\mathbf{s}_{t-1}\right)P\left(\mathbf{s}_{t-1}\right).$
(3)

This marginal distribution $P(\mathbf{s}_{t})$ is not factorized (except at
$\mathbf{J}=\mathbf{0}$), but it rather exhibits a complex statistical
structure, generally containing higher-order spin interactions. We can apply
this equation recursively, e.g. decomposing $P\left(\mathbf{s}_{t-1}\right)$
in terms of the distribution $P(\mathbf{s}_{t-2})$, and trace the evolution of
the system from the initial distribution $P(\mathbf{s}_{0})$.

In this article, we use variants of the Plefka expansion to calculate some
statistical properties of the system. Namely, we investigate the average
activation rates $\mathbf{m}_{t}$, correlations between pairs of units
(covariance function) $\mathbf{C}_{t}$, and delayed correlations
$\mathbf{D}_{t}$ given by

$\displaystyle m_{i,t}=$
$\displaystyle\sum_{\mathbf{s}_{t}}s_{i,t}P(\mathbf{s}_{t}),$ (4)
$\displaystyle C_{ik,t}=$
$\displaystyle\sum_{\mathbf{s}_{t}}s_{i,t}s_{k,t}P(\mathbf{s}_{t})-m_{i,t}m_{k,t},$
(5) $\displaystyle D_{il,t}=$
$\displaystyle\sum_{\mathbf{s}_{t},\mathbf{s}_{t-1}}s_{i,t}s_{l,t-1}P(\mathbf{s}_{t},\mathbf{s}_{t-1})-m_{i,t}m_{l,t-1}.$
(6)

Note that $\mathbf{m}_{t}$ and $\mathbf{D}_{t}$ are sufficient statistics of
the kinetic Ising model. Therefore, we will use them in solving the inverse
Ising problem (see Methods). We additionally consider the equal-time
correlations $\mathbf{C}_{t}$ as they are commonly used to describe the
systems, and are investigated by some of the mean-field approximations in the
literature [40]. Calculation of these expectation values is analytically
intractable and computationally very expensive for large networks, due to the
combinatorial explosion of the number of possible states. To reduce this
computational cost, we approximate the marginal probability distributions (Eq.
3) by the Plefka expansion method that utilizes an alternative, tractable
distribution.

### Geometrical approach to mean-field approximation

Information geometry [37, 45, 46] provides clear geometrical understanding of
information-theoretic measures and probabilistic models [47, 48, 15]. Using
the language of information geometry, we introduce our method for mean-field
approximations of kinetic Ising systems.

Let $\mathcal{P}_{t}$ be the manifold of probability distributions at time $t$
obtained from Eq. 3. Each point on the manifold corresponds with a set of
parameter values. The manifold $\mathcal{P}_{t}$ contains sub-manifolds
$\mathcal{Q}_{t}$ of probability distributions with analytically tractable
statistical properties (See Fig. 1). We use this tractable manifold, i.e., a
reference model, to approximate a target point
$P(\mathbf{s}_{t}|\mathbf{H},\mathbf{J})$ in the manifold $\mathcal{P}_{t}$
and its statistical properties $\mathbf{m}_{t},\mathbf{C}_{t},\mathbf{D}_{t}$.

The simplest sub-manifold $\mathcal{Q}_{t}$ is the manifold of independent
models, used in classical mean-field approximations to compute average
activation rates. Each point on this sub-manifold corresponds to a
distribution

$Q\left(\mathbf{s}_{t}|\mathbf{\Theta}_{t}\right)=\prod_{i}\frac{\mathrm{e}^{s_{i,t}\Theta_{i,t}}}{2\cosh\Theta_{i,t}},$
(7)

where $\bm{\Theta}_{t}=\\{\Theta_{i,t}\\}$ is the vector of parameters that
represents a point in $\mathcal{Q}_{t}$. This distribution does not include
couplings between units, and its average activation rate is immediately given
as $m_{i,t}=\tanh\Theta_{i,t}$.

Our first goal is to find the average activation rates of the target
distribution $P(\mathbf{s}_{t}|\mathbf{H},\mathbf{J})$. It turns out that they
can be obtained from the independent model $Q(\mathbf{s}_{t}|\bm{\Theta}_{t})$
that minimizes the following Kullback-Leibler (KL) divergence from
$P(\mathbf{s}_{t})$:

$\displaystyle D(P||Q)=$
$\displaystyle\sum\limits_{\mathbf{s}_{t}}P(\mathbf{s}_{t})\log\frac{P(\mathbf{s}_{t})}{Q(\mathbf{s}_{t})}.$
(8)

The independent model $Q(\mathbf{s}_{t}|\bm{\Theta}_{t}^{*})(\equiv
Q^{*}(\mathbf{s}_{t}))$ that minimizes the KL divergence has activation rates
$\mathbf{m}_{t}$ identical to those of the target distribution
$P(\mathbf{s}_{t})$ [38] because the minimizing points $\Theta^{*}_{i,t}$
satisfy (for $i=1,\ldots,N$)

$\displaystyle\frac{\partial
D(P||Q)}{\partial\Theta_{i,t}}\bigg{\rvert}_{\bm{\Theta}_{t}=\bm{\Theta}_{t}^{*}}=$
$\displaystyle-\sum\limits_{\mathbf{s}_{t}}s_{i,t}P(\mathbf{s}_{t})+\tanh\Theta^{*}_{i,t}$
$\displaystyle=$ $\displaystyle m_{i,t}^{Q^{*}}-m_{i,t}^{P}=0,$ (9)

where $m_{i,t}^{P}$ and $m_{i,t}^{Q^{*}}$ are respectively expectation values
of $s_{i,t}$ by $P(\mathbf{s}_{t})$ and
$Q(\mathbf{s}_{t}|\mathbf{\Theta}_{t}^{*})$. As these values are equal, for
the rest of the paper we will drop their superscripts and just write $m_{i,t}$
for simplicity. Note that the result of this approximation is indifferent of
the system’s correlations. Later in the paper we will consider approximations
that take into account pairwise correlations.

From an information geometric point of view, given $\mathbf{m}_{t}$ (or
$\mathbf{\Theta}_{t}^{*}$), we may consider a family of points defined as a
linear mixture of $P(\mathbf{s}_{t})$ and
$Q(\mathbf{s}_{t}|\mathbf{\Theta}_{t}^{*})$ for which $\mathbf{m}_{t}$ is kept
constant (the dashed line $\mathcal{A}$ in Fig. 1). This is known as an
m-geodesic, and it is orthogonal to the e-flat manifold $\mathcal{Q}_{t}$,
constituting an m-projection to this manifold [47, 37]. Thus, the previous
search of $Q(\mathbf{s}_{t}|\mathbf{\Theta}_{t}^{*})$ given
$P(\mathbf{s}_{t}|\mathbf{H},\mathbf{J})$ is equivalent to finding the
orthogonal projection point from $P(\mathbf{s}_{t}|\mathbf{H},\mathbf{J})$ to
the manifold $\mathcal{Q}_{t}$ of independent models [36, 37].

Figure 1: A geometric view of the approximations based on Plefka expansions.
The point $P(\mathbf{s}_{t})$ is the marginal distribution of a kinetic Ising
model at time $t$. The submanifold $\mathcal{Q}_{t}$ is a set of tractable
distributions, for example a manifold of independent models. The points in
$\mathcal{A}$ correspond to a m-geodesic, that is a linear mixture of
$P(\mathbf{s}_{t})$ and $Q^{*}(\mathbf{s}_{t})$ on $\mathcal{Q}_{t}$, where
for independent $\mathcal{Q}_{t}$ all points on $\mathcal{A}$ share the same
mean values $\mathbf{m}_{t}$. Geometrically, $\mathcal{A}$ constitutes the
m-projection from $P(\mathbf{s}_{t})$ to $\mathcal{Q}_{t}$, defining
$Q^{*}(\mathbf{s}_{t})$ as the closest point in the submanifold
$\mathcal{Q}_{t}$ to the point $P(\mathbf{s}_{t})$ [47]. The Plefka expansion
is defined by expanding an $\alpha$-dependent distribution
$P_{\alpha}(\mathbf{s}_{t})$ that satisfies
$P_{\alpha=0}(\mathbf{s}_{t})=Q^{*}(\mathbf{s}_{t})$ and
$P_{\alpha=1}(\mathbf{s}_{t})=P(\mathbf{s}_{t})$.

### The Plefka expansion

Although the m-projection provides the exact and unique average activation
rates, its calculation in practice requires the complete distribution
$P(\mathbf{s}_{t})$. In the Plefka expansion, we relax the constraints of the
m-projection, and introduce another set of more tractable distributions that
passes only through $P(\mathbf{s}_{t}|\mathbf{H},\mathbf{J})$ and
$Q(\mathbf{s}_{t}|\mathbf{\Theta}_{t}^{*})$ (the solid line in Fig. 1). This
distribution is defined using a new conditional distribution introducing a
parameter $\alpha$ that connects a distribution on the manifold
$\mathcal{Q}_{t}$ with the original distribution $P(\mathbf{s}_{t})$:

$\displaystyle P_{\alpha}(s_{i,t}|\mathbf{s}_{t-1})=$
$\displaystyle\frac{\mathrm{e}^{s_{i,t}{h}_{i,t}(\alpha)}}{2\cosh{h}_{i,t}(\alpha)},$
(10) $\displaystyle h_{i,t}(\alpha)=$
$\displaystyle(1-\alpha)\Theta_{i,t}+\alpha(H_{i}+\sum_{j}J_{ij}s_{j,t-1}).$
(11)

At $\alpha=0$,
$P_{\alpha=0}(\mathbf{s}_{t}|\mathbf{s}_{t-1})=Q(\mathbf{s}_{t}|\mathbf{\Theta}_{t})$,
and $\alpha=1$ leads to
$P_{\alpha=1}(\mathbf{s}_{t}|\mathbf{s}_{t-1})=P(\mathbf{s}_{t}|\mathbf{s}_{t-1})$.
Using this alternative conditional distribution
$P_{\alpha}(s_{i,t}|\mathbf{s}_{t-1})$, we construct an approximate marginal
distribution $P_{\alpha}(\mathbf{s}_{t})$. Consequently, expectation values
with respect to $P_{\alpha}(\mathbf{s}_{t})$ are functions of $\alpha$. We
thus write the statistics of the approximate system as
$\mathbf{m}_{t}(\alpha)$, $\mathbf{C}_{t}(\alpha)$, and
$\mathbf{D}_{t}(\alpha)$.

The Plefka expansions of these statistics are defined as the Taylor series
expansion of these functions around $\alpha=0$. In the case of the mean
activation rate, the expansion up to the $n$th-order leads to:

$\mathbf{m}_{t}(\alpha)=\mathbf{m}_{t}(\alpha=0)+\sum_{k=1}^{n}\frac{\alpha^{k}}{k!}\frac{\partial^{k}\mathbf{m}_{t}(\alpha=0)}{\partial\alpha^{k}}+\mathcal{O}(\alpha^{(n+1)}),$
(12)

where $\mathcal{O}(\alpha^{(n+1)})$ stands for the residual error of the
approximation of order $n+1$ and higher. For the $n$th-order approximation, we
neglect the residual terms as
$\mathcal{O}(\alpha^{(n+1)})\big{\rvert}_{\alpha=1}\approx 0$. Note that all
coefficients of expansion are functions of $\mathbf{\Theta}_{t}$. The mean-
field approximation is computed by setting $\alpha=1$ and finding the value of
$\mathbf{\Theta}_{t}^{*}$ that satisfies Eq. 12. Note that since the original
marginal distribution is recovered at $\alpha=1$, the equality of Eq. 9 holds:
$\mathbf{m}_{t}(\alpha=1)=\mathbf{m}_{t}(\alpha=0)$. Then, we have

$\sum_{k=1}^{n}\frac{1}{k!}\frac{\partial^{k}\mathbf{m}_{t}(\alpha=0)}{\partial\alpha^{k}}=0,$
(13)

which should be solved with respect to the parameters $\mathbf{\Theta}_{t}$.
Since we neglected the terms higher than the $n$-th order, the solution may
not lead to the exact projection, $Q(\mathbf{s}_{t}|\mathbf{\Theta}_{t}^{*})$.
In this study, we investigate the first ($n=1$) and second ($n=2$) order
approximations. Moreover we can apply the same expansion to approximate the
correlations $\mathbf{C}_{t}$ and $\mathbf{D}_{t}$, using Eq. 10.

What is the difference between this approach and other mean-field methods?
Conventionally, naive mean-field approximations are obtained by minimizing
$D(Q||P)$ as opposed to $D(P||Q)$ (Eq. 8) [49, 36]. This approach is typically
used in variational inference to construct a tractable approximate posterior
in machine learning problems. Following the Bogolyubov inequality, minimizing
this divergence is equivalent to minimizing the variational free energy.
Geometrically, it comprises an e-projection of
$P(\mathbf{s}_{t}|\mathbf{H},\mathbf{J})$ to the sub-manifold
$\mathcal{Q}_{t}$, which does not result in
$Q(\mathbf{s}_{t}|\mathbf{\Theta}_{t}^{*})$. Namely, minimizing $D(Q||P)$, as
well as minimization of other $\alpha$-divergences except for $D(P||Q)$,
introduces a bias in the estimation of the mean-field approximation [36, 37].
In contrast, if we consider the m-projection point that minimizes $D(P||Q)$,
we can approximate the exact value of $\mathbf{m}_{t}$ using Eq. 12 up to an
arbitrary order.

In the subsequent sections we show that different approximations of the
marginal distribution $P(\mathbf{s}_{t})$ in Eq. 3 can be constructed by
replacing $P(s_{i,\tau}|\mathbf{s}_{\tau-1})$ with
$P_{\alpha}(s_{i,\tau}|\mathbf{s}_{\tau-1})$ for different pairs $i,\tau$
(here we will explore the cases of $\tau=t$ and $\tau=t-1$). More generally,
we show in Supplementary Note 1 that this framework can be extended to a
marginal path of arbitrary length $k$,
$P(\mathbf{s}_{t-k+1},\dots,\mathbf{s}_{t})$. In addition, we are not
restricted to manifolds of independent models. The independent model is
adopted as a reference model to approximate the average activation rate, but
one can also more accurately approximate correlations using this method. In
this vein, we can extend our framework to use reference manifolds
$\mathcal{Q}_{t-k+1:t}$ (of models
$Q(\mathbf{s}_{t-k+1},\dots,\mathbf{s}_{t})$) that include interactions, e.g.,
pairwise couplings between elements at two different time points, to more
accurately approximate the delayed correlations (see Supplementary Note 1). By
systematically defining these reference distributions, we will provide a
unified approach to derive Plefka approximations of $\mathbf{m}_{t}$,
$\mathbf{C}_{t}$, and $\mathbf{D}_{t}$, including the one that utilizes a
pairwise structure.

### Plefka[$t-1,t$]: Expansion around independent models at times $t-1$ and
$t$

Before elaborating different mean-field approximations, we demonstrate our
method by deriving the known results of the classical nMF and TAP
approximations for the kinetic Ising model [38, 39]. In order to derive these
classical mean-field equations, we make a Plefka expansion around the points
$\mathbf{\Theta}^{*}_{t}$ and $\mathbf{\Theta}^{*}_{t-1}$ that are
respectively obtained by orthogonal projection to the independent manifolds
$\mathcal{Q}_{t}$ and $\mathcal{Q}_{t-1}$, computed as in Eq. 9. Note that
assuming an approximation where previous distributions (e.g., $t-2,t-3,\dots$)
are also independent yields exactly the same result. In this way, we derive
the nMF and TAP equations of a model defined by a marginal probability
distribution $P_{\alpha}^{[t-1:t]}$. Using Eqs. (3) and (10), we write

$P_{\alpha}^{[t-1:t]}(\mathbf{s}_{t})=\sum_{\begin{subarray}{c}\mathbf{s}_{t-1}\\\
\mathbf{s}_{t-2}\end{subarray}}P_{\alpha}(\mathbf{s}_{t}|\mathbf{s}_{t-1})P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$
(14)

where $P_{\alpha=0}^{[t-1:t]}(\mathbf{s}_{t})=Q(\mathbf{s}_{t})$ and the
original distribution is recovered for
$P_{\alpha=1}^{[t-1:t]}(\mathbf{s}_{t})=P(\mathbf{s}_{t})$.

Following Eq. 13, for the first order approximation we have $\frac{\partial
m_{i,t}(\alpha=0)}{\partial\alpha}=0$. Since the derivative of the first order
moment is

$\displaystyle\left.\frac{\partial m_{i,t}(\alpha=0)}{\partial\alpha}\right.=$
$\displaystyle(1-m_{i,t}^{2})(-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}),$
(15)

by solving the equation, we find $\Theta^{*}_{i,t}\approx
H_{i}+\sum_{j}J_{ij}m_{j,t-1}$ that leads to the naive mean-field
approximation:

$m_{i,t}\approx\tanh[H_{i}+\sum_{j}J_{ij}m_{j,t-1}].$ (16)

We apply the same expansion to approximate the correlations, expanding
$C_{ik,t}(\alpha)$ and $D_{il,t}(\alpha)$ around $\alpha=0$ up to the first
order using $\Theta_{i,t}=\Theta^{*}_{i,t}$. Then we obtain

$\displaystyle C_{ik,t}\approx$ $\displaystyle 0,\quad i\neq k,$ (17)
$\displaystyle D_{il,t}\approx$ $\displaystyle
J_{il}(1-m_{i,t}^{2})(1-m_{l,t-1}^{2}).$ (18)

Detailed calculations are presented in Supplementary Note 2.

To obtain the second order approximation, we need to solve $\frac{\partial
m_{i}(\alpha=0)}{\partial\alpha}+\frac{1}{2}\frac{\partial^{2}m_{i}(\alpha=0)}{\partial\alpha^{2}}=0$
from Eq. 13. Here the second order derivative is given as

$\displaystyle\frac{\partial^{2}m_{i,t}(\alpha=0)}{\partial\alpha^{2}}$
$\displaystyle\approx-2m_{i,t}(1-m_{i,t}^{2})\sum_{j}J_{ij}^{2}(1-m_{j,t-1}^{2}),$
(19)

where terms of the order higher than quadratic were neglected (see
Supplementary Note 2 for further details). From these equations, we find
$\Theta^{*}_{i,t}\approx
H_{i}+\sum_{j}J_{ij}m_{j,t-1}-m_{i,t}\sum_{j}J_{ij}^{2}(1-m_{j,t-1}^{2})$
leading to the TAP equation:

$m_{i,t}\approx\tanh[H_{i}+\sum_{j}J_{ij}m_{j,t-1}-m_{i,t}\sum_{j}J_{ij}^{2}(1-m_{j,t-1}^{2})].$
(20)

Having $\Theta^{*}_{i,t}$, we can incorporate TAP approximations of the
correlations by expanding $C_{ik,t}(\alpha)$ and $D_{il,t}(\alpha)$ (see
Supplementary Note 2 for details) as:

$\displaystyle C_{ik,t}\approx$
$\displaystyle(1-m_{i,t}^{2})(1-m_{k,t}^{2})\sum\limits_{j}J_{ij}J_{kj}(1-m_{j,t-1}^{2}),$
$\displaystyle\quad i\neq k,$ (21) $\displaystyle D_{il,t}\approx$
$\displaystyle
J_{il}(1-m_{i,t}^{2})(1-m_{l,t-1}^{2})(1+2J_{il}m_{i,t}m_{l,t-1}).$ (22)

In these approximations, Eqs. 16 and 20 of activation rates $\mathbf{m}_{t}$
correspond to the classical nMF and TAP equations of the kinetic Ising model
[38, 39]. The mean-field equations for the equal-time and delayed correlations
(Eqs. 17, 18 and 21, 22) are novel contributions from applying the Plefka
expansion to correlations.

We note that, using the equations above, we can compute the approximate
statistical properties of the system at $t$
($\mathbf{m}_{t},\mathbf{C}_{t},\mathbf{D}_{t}$) from $\mathbf{m}_{t-1}$.
Therefore, the system evolution is described by recursively computing
$\mathbf{m}_{t}$ from an initial state $\mathbf{m}_{0}$ (for both transient
and stationary dynamics), although approximation errors accumulate over the
iterations. After we introduce a unified view of mean-field approximations in
the subsequent sections, we will numerically examine approximation errors of
these various methods in predicting statistical structure of the system.

### Generalization of mean-field approximations

In the previous section, we described a Plefka approximation that uses a model
containing independent units at time $t-1$ and $t$ to construct a marginal
probability distribution $P_{\alpha}^{[t-1:t]}(\mathbf{s}_{t})$. This is,
however, not the only possible choice of approximation. As we mentioned above,
other approximations have been introduced in the literature. In [40],
expressions are provided for the non-stationary delayed correlations
$\mathbf{D}_{t}$ as a function of $\mathbf{C}_{t-1}$. In [41], an
approximation is derived by assuming that units at state $\mathbf{s}_{t-1}$
are independent while correlations of $\mathbf{s}_{t}$ are preserved.

Figure 2: Unified mean-field framework. Original model (A) and family of
generalized Plefka expansions (B-E). Gray lines represent connections that are
proportional to $\alpha$ and thus removed in the approximated model to perform
the Plefka expansions, while solid black lines are conserved and dashed lines
are free parameters. Plefka[$t-1,t$] (B) retrieves the classical mean-field
equations [38, 39]. Plefka[$t$] (C) results in a novel method which preserves
correlations of the system at $t-1$, incorporating equations similar to [40].
Plefka[$t-1$] (D) assumes independent activity at t-1, and in its first order
approximation reproduces the results in [41]. Plefka2[$t$] (E) represents a
novel pairwise approximation which performs better in approximating
correlations.

In the following sections, we show that various approximation methods,
including those mentioned above, can be unified as Plefka expansions. Each
method of the approximation corresponds to a specific choice of the sub-
manifold $\mathcal{Q}_{t}$ at each time step. Fig. 2 show the corresponding
sub-manifolds $\mathcal{Q}_{t-1:t}$ of possible approximations, where gray
lines represent interactions that are affected by $\alpha$ in the Plefka
expansion. The mean-field approximations in the previous section were obtained
by using the model represented in Fig. 2B, where the couplings at time $t-1$
and $t$ are affected by $\alpha$. Below, we present systematic applications of
the Plefka expansions around other reference models in order to approximate
the original distribution (Fig. 2C-E). By doing so, we not only unify the
previously reported mean-field approximations but also provide novel solutions
that can provide more precise approximations than known methods.

### Plefka[$t$]: Expansion around an independent model at time $t$

For the Plefka[$t-1,t$] approximation, explained above, the system becomes
independent for $\alpha=0$ at $t$ as well as $t-1$. This leads to
approximations of $\mathbf{m}_{t},\mathbf{C}_{t},\mathbf{D}_{t}$ being
specified by $\mathbf{m}_{t-1}$, while being independent of $\mathbf{C}_{t-1}$
and $\mathbf{D}_{t-1}$. In [40], the authors describe a mean-field
approximation by performing new expansion over the classical nMF and TAP
equations that takes into account previous correlations $\mathbf{C}_{t-1}$.
Here, our framework allows us to obtain similar results by considering only a
Plefka expansion over manifold $\mathcal{Q}_{t}$ while assuming that we know
the properties of $P(\mathbf{s}_{t-1})$ (Fig. 2C). Therefore, we denote this
approximation as $P^{[t]}_{\alpha}$ and consider

$P_{\alpha}^{[t]}(\mathbf{s}_{t})=\sum_{\mathbf{s}_{t-1}}P_{\alpha}(\mathbf{s}_{t}|\mathbf{s}_{t-1})P(\mathbf{s}_{t-1}).$
(23)

In Supplementary Note 3 we derive the equations for this approximation. For
the first order, we obtain

$\displaystyle m_{i,t}\approx$
$\displaystyle\tanh[H_{i}+\sum_{j}J_{ij}m_{j,t-1}],$ (24) $\displaystyle
C_{ik,t}\approx$ $\displaystyle 0,\quad i\neq k,$ (25) $\displaystyle
D_{il,t}\approx$
$\displaystyle(1-m_{i,t}^{2})\sum\limits_{j}J_{ij}C_{jl,t-1}.$ (26)

Note that Eqs. 24 and 25 are the same as the nMF Plefka[$t-1,t$] equations.
Eq. 26 includes $\mathbf{C}_{t-1}$, being exactly the same result obtained in
[40, Eq. 4]. The second order approximations leads to:

$\displaystyle m_{i,t}\approx$
$\displaystyle\tanh[H_{i}+\sum_{j}J_{ij}m_{j,t-1}-m_{i,t}\sum_{jl}J_{ij}J_{il}C_{jl,t-1}],$
(27) $\displaystyle C_{ik,t}\approx$
$\displaystyle(1-m_{i,t}^{2})(1-m_{k,t}^{2})\sum_{jl}J_{ij}J_{kl}C_{jl,t-1},\quad
i\neq k,$ (28) $\displaystyle D_{il,t}\approx$
$\displaystyle(1-m_{i,t}^{2})\sum\limits_{j}J_{ij}C_{jl,t-1}(1+2J_{il}m_{i,t}m_{l,t-1}).$
(29)

All update rules include the effect of $\mathbf{C}_{t-1}$. We can see that if
we use the covariance matrix of the independent model at $t-1$, we recover the
results of $P_{\alpha}^{[t-1:t]}$ in the previous section. In contrast with
[40], we provide a novel approximation method that depends on previous
correlations using a single expansion (instead of two subsequent expansions),
and additionally present approximated equal-time correlations.

### Plefka[$t-1$]: Expansion around an independent model at time $t-1$

In [41], a mean-field method is proposed by approximating the effective field
$\mathbf{h}_{t}$ as the sum of a large number of independent spins,
approximated by a Gaussian distribution, yielding exact results for fully
asymmetric networks in the thermodynamic limit. In our framework, we describe
this approximation as an expansion around the projection point from
$P(\mathbf{s}_{t-1})$ to the submanifold $\mathcal{Q}_{t-1}$, using a model
where only $\mathbf{s}_{t-1}$ are independent (Fig. 2D). In this case (see
Supplementary Note 4), the effective field $\mathbf{h}_{t}$ at the submanifold
is a sum of independent terms, which for large $N$ yields a Gaussian
distribution.

By defining

$P_{\alpha}^{[t-1]}(\mathbf{s}_{t},\mathbf{s}_{t-1})=\sum_{\mathbf{s}_{t-2}}P(\mathbf{s}_{t}|\mathbf{s}_{t-1})P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$
(30)

we see that now the expansion is defined for the marginal distribution of the
path $\\{\mathbf{s}_{t-1},\mathbf{s}_{t}\\}$ (see Supplementary Note 1). The
first order equations for this method are

$\displaystyle m_{i,t}\approx$
$\displaystyle\int\mathrm{D}_{x}\tanh[H_{i}+\sum_{j}J_{ij}m_{j,t-1}+x\sqrt{\Delta_{i,t}}],$
(31) $\displaystyle C_{ik,t}\approx$
$\displaystyle\int\mathrm{D}_{xy}^{\rho_{ik}}\tanh[H_{i}+\sum_{j}J_{ij}m_{j,t-1}+x\sqrt{\Delta_{i,t}}]$
$\displaystyle\cdot\tanh[H_{k}+\sum_{l}J_{kl}m_{l,t-1}+y\sqrt{\Delta_{j,t}}]$
$\displaystyle-m_{i,t}m_{k,t},\quad i\neq k,$ (32) $\displaystyle
D_{il,t}\approx$
$\displaystyle\sum\limits_{j}J_{ij}C_{jl,t-1}\int\mathrm{D}_{x}\big{(}1$
$\displaystyle-\tanh^{2}[H_{i}+\sum_{j}J_{ij}m_{j,t-1}+x\sqrt{\Delta_{i,t}}]\big{)}.$
(33)

Here we use
$\mathrm{D}_{x}=\frac{\mathrm{d}x}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}x^{2}\right)$,
$\mathrm{D}_{xy}^{\rho_{ik}}=\frac{\mathrm{d}x\mathrm{d}y}{2\pi\sqrt{1-\rho_{ik}^{2}}}\exp\left(-\frac{1}{2}\frac{(x^{2}+y^{2})-2\rho_{ik}xy}{1-\rho_{ik}^{2}}\right)$,
$\Delta_{i,t}=\sum_{j}J_{ij}^{2}(1-m_{j,t-1}^{2})$ and
$\rho_{xy}=\sum_{j}J_{ij}J_{kj}(1-m_{j,t-1}^{2})/\sqrt{\Delta_{i,t}\Delta_{j,t}}$.
Derivations are described in Supplementary Note 4. These results are exactly
the same as those presented for $\mathbf{m}_{t},\mathbf{D}_{t}$ in [41],
adding an additional expression for $\mathbf{C}_{t}$. For this approximation,
we do not consider the second order equations since they are computationally
much more expensive than the other approximations.

### Plefka2[$t$]: Expansion around a pairwise model

The proposed framework is also a powerful tool to develop novel Plefka
expansions. To make the expansions more accurately approximate target
statistics, we can consider a reference manifold composed of multiple time
steps while maintaining some of the parameters in the system (see
Supplementary Note 1). Motivated by this idea, here we propose new methods
that directly approximate pairwise activities of the units by choosing a
reference manifold that preserves a coupling term.

Let us first consider the joint probability of any arbitrary pair of units at
time $t-1$ and $t$ to compute the delayed correlations (Fig. 2E, left).
Namely, we consider the joint probability of spins $s_{i,t}$ and $s_{l,t-1}$:

$\displaystyle P(s_{i,t},s_{l,t-1})=$
$\displaystyle\sum_{\begin{subarray}{c}\mathbf{s}_{\setminus l,t-1}\\\
\mathbf{s}_{t-2}\end{subarray}}P(s_{i,t}|\mathbf{s}_{t-1})P(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$
(34)

with $\mathbf{s}_{\setminus l,t-1}$ containing all elements of
$\mathbf{s}_{t-1}$ except $s_{l,t-1}$. As a reference manifold
$\mathcal{Q}_{t-1:t}$, we consider the dependency among only the units $i$ and
$l$:

$\displaystyle Q(s_{i,t},s_{l,t-1})=$ $\displaystyle
Q(s_{i,t}|s_{l,t-1})Q(s_{l,t-1})$ $\displaystyle=$
$\displaystyle\frac{\mathrm{e}^{s_{i,t}\theta_{i,t}(s_{l,t-1})}}{2\cosh\theta_{i,t}(s_{l,t-1})}\frac{\mathrm{e}^{s_{l,t-1}\Theta_{l,t-1}}}{2\cosh\Theta_{l,t-1}},$
(35)

where $\theta_{i,t}(s_{l,t-1})=\Theta_{i,t}+\Delta_{il,t}s_{l,t-1}$.
Orthogonal projection to $\mathcal{Q}_{t}$ is equivalent to minimizing the KL
divergence $D(P||Q)$ with respect to the parameters:

$\displaystyle\frac{\partial
D(P||Q)}{\partial\Theta_{i,t}}\bigg{\rvert}_{\begin{subarray}{l}\bm{\theta}_{t}=\bm{\theta}_{t}^{*}\\\
\bm{\Theta}_{t-1}=\bm{\Theta}_{t-1}^{*}\end{subarray}}=$ $\displaystyle
m_{i,t}^{Q^{*}}-m_{i,t}^{P}=0,$ (36) $\displaystyle\frac{\partial
D(P||Q)}{\partial\Theta_{l,t-1}}\bigg{\rvert}_{\begin{subarray}{l}\bm{\theta}_{t}=\bm{\theta}_{t}^{*}\\\
\bm{\Theta}_{t-1}=\bm{\Theta}_{t-1}^{*}\end{subarray}}=$ $\displaystyle
m_{l,t-1}^{Q^{*}}-m_{l,t-1}^{P}=0,$ (37) $\displaystyle\frac{\partial
D(P||Q)}{\partial\Delta_{il,t}}\bigg{\rvert}_{\begin{subarray}{l}\bm{\theta}_{t}=\bm{\theta}_{t}^{*}\\\
\bm{\Theta}_{t-1}=\bm{\Theta}_{t-1}^{*}\end{subarray}}=$ $\displaystyle\langle
s_{i,t}s_{l,t-1}\rangle_{(Q^{*}\cdot P)}-\langle s_{i,t}s_{l,t-1}\rangle_{P}$
$\displaystyle=0.$ (38)

with

$\langle x\rangle_{(Q^{*}\cdot P)}=\sum_{\begin{subarray}{c}s_{i,t}\\\
s_{l,t-1}\end{subarray}}x\,Q(s_{i,t}|s_{l,t-1},\bm{\theta}_{t}^{*},\bm{\Theta}_{t-1}^{*})P(s_{l,t-1}).$
(39)

As in the previous approximations, $P(s_{i,t},s_{l,t-1})$ is connected to
$Q(s_{i,t},s_{l,t-1}|\bm{\theta}_{t}^{*},\bm{\Theta}_{t-1}^{*})$ through an
$\alpha$-dependent probability

$\displaystyle P_{\alpha}^{2[t]}(s_{i,t},s_{l,t-1})=$
$\displaystyle\sum_{\begin{subarray}{c}\mathbf{s}_{\setminus l,t-1}\\\
\mathbf{s}_{t-2}\end{subarray}}P_{\alpha}(s_{i,t}|\mathbf{s}_{t-1})P_{\alpha}(s_{l,t-1}|\mathbf{s}_{t-2})$
$\displaystyle\cdot P(\mathbf{s}_{\setminus
l,t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$ (40)

with conditional probabilities given by

$\displaystyle P_{\alpha}(s_{i,t}|\mathbf{s}_{t-1})=$
$\displaystyle\frac{\mathrm{e}^{s_{i,t}h_{i,t}(\alpha)}}{2\cosh
h_{i,t}(\alpha)},$ (41) $\displaystyle h_{i,t}(\alpha)=$
$\displaystyle(1-\alpha)\theta_{i,t}(s_{l,t-1})$
$\displaystyle+\alpha(H_{i}+\sum_{j}J_{ij}s_{j,t-1}),$ $\displaystyle
P_{\alpha}(s_{l,t-1}|\mathbf{s}_{t-2})=$
$\displaystyle\frac{\mathrm{e}^{s_{l,t-1}{h}_{l,t-1}(\alpha)}}{2\cosh{h}_{l,t-1}(\alpha)},$
(42) $\displaystyle h_{l,t-1}(\alpha)=$
$\displaystyle(1-\alpha)\Theta_{l,t-1}+\alpha(H_{l}+\sum_{n}J_{ln}s_{n,t-2}).$

As in the cases above, we can calculate the equations for the first and second
order approximations (see Supplementary Note 5). Here, for the second order
approximation (which is more accurate than the first order) we have that:

$\displaystyle\theta_{i,t}^{*}(s_{l,t-1})\approx$ $\displaystyle
H_{i}+\sum_{j}J_{ij}m_{j,t-1}+J_{il}(s_{l,t-1}-m_{l,t-1})$
$\displaystyle+\big{(}\sum_{j\neq
l,n}J_{ij}J_{ln}D_{jn,t-1}\big{)}(s_{l,t-1}-m_{l,t-1})$
$\displaystyle-\tanh\theta_{i,t}^{*}(s_{l,t-1})\sum_{jn\neq
l}J_{ij}J_{in}C_{jn,t-1},$ (43) $\displaystyle\Theta_{l,t-1}^{*}\approx$
$\displaystyle H_{l}+\sum_{n}J_{ln}m_{n,t-2}$ $\displaystyle-
m_{l,t-1}\sum_{mn}J_{lm}J_{ln}C_{mn,t-2},$ (44)

which directly leads to calculation of means and delayed correlations as:

$\displaystyle m_{i,t}\approx$
$\displaystyle\sum_{s_{l,t-1}}\tanh\theta_{i,t}^{*}(s_{l,t-1})Q^{*}(s_{l,t-1}),$
(45) $\displaystyle m_{l,t-1}\approx$ $\displaystyle\tanh\Theta_{l,t-1}^{*},$
(46) $\displaystyle D_{il,t}\approx$
$\displaystyle\sum_{s_{l,t-1}}\tanh\theta_{i,t}^{*}(s_{l,t-1})s_{l,t-1}Q^{*}(s_{l,t-1})$
$\displaystyle-m_{i,t}m_{l,t-1}.$ (47)

These results are related to previous work [43] that included autocorrelations
as one of the constraints to derive the Plefka approximation. Instead, here we
provide a Plefka approximation that includes delayed correlations between any
pair of units.

To compute the above approximations, we need to know $\mathbf{C}_{t-1}$ and
$\mathbf{C}_{t-2}$. Here, we provide similar pairwise Plefka approximations
for the pairwise distribution at time $t$, $P(s_{i,t},s_{k,t})$. Since
$s_{i,t},s_{k,t}$ are conditionally independent, we can construct a model in
which first $s_{k,t}$ is computed from $\mathbf{s}_{t-1}$, and then $s_{i,t}$
is computed conditioned on $s_{k,t},\mathbf{s}_{t-1}$ (Fig. 2E, right):

$\displaystyle Q(s_{i,t},s_{k,t})=$ $\displaystyle
Q(s_{i,t}|s_{k,t})Q(s_{k,t}),$ (48) $\displaystyle
P_{\alpha}^{2[t]}(s_{i,t},s_{k,t})=$
$\displaystyle\sum_{\mathbf{s}_{t-1}}P_{\alpha}(s_{i,t}|s_{k,t},\mathbf{s}_{t-1})$
(49) $\displaystyle\cdot
P_{\alpha}(s_{k,t}|\mathbf{s}_{t-1})P(\mathbf{s}_{t-1}),$ (50)

with conditional probabilities given by

$\displaystyle P_{\alpha}(s_{i,t}|s_{k,t},\mathbf{s}_{t-1})=$
$\displaystyle\frac{\mathrm{e}^{s_{i,t}h_{i,t}(\alpha)}}{2\cosh
h_{i,t}(\alpha)},$ (51) $\displaystyle h_{i,t}(\alpha)=$
$\displaystyle(1-\alpha)\theta_{i,t}(s_{k,t})$
$\displaystyle+\alpha(H_{i}+\sum_{j}J_{ij}s_{j,t-1}),$ $\displaystyle
P_{\alpha}(s_{k,t}|\mathbf{s}_{t-1})=$
$\displaystyle\frac{\mathrm{e}^{s_{k,t}{h}_{k,t}(\alpha)}}{2\cosh{h}_{k,t}(\alpha)},$
(52) $\displaystyle h_{k,t}(\alpha)=$
$\displaystyle(1-\alpha)\Theta_{k,t}+\alpha(H_{k}+\sum_{l}J_{kl}s_{l,t-1}).$

Note that $\theta_{i,t}$ is a function of $s_{k,t}$ that accounts for equal-
time correlations between $s_{i,t}$ and $s_{k,t}$. Computed similarly to
delayed correlations, the second order approximation yields (see Supplementary
Note 5):

$\displaystyle\theta_{i,t}^{*}(s_{k,t})\approx$ $\displaystyle
H_{i}+\sum_{j}J_{ij}m_{j,t-1}$
$\displaystyle+(\sum_{jl}J_{ij}J_{kl}C_{jl,t-1})(s_{k,t}-m_{k,t})$
$\displaystyle-\tanh\theta_{i,t}^{*}(s_{k,t})\sum_{jl}J_{ij}J_{il}C_{jl,t-1},$
(53) $\displaystyle\Theta_{k,t}^{*}\approx$ $\displaystyle
H_{k}+\sum_{l}J_{kl}m_{l,t-1}$ $\displaystyle-
m_{k,t}\sum_{jl}J_{kj}J_{kl}C_{jl,t-1}.$ (54)

Using these equations, approximate equal-time correlations are given as

$\displaystyle C_{ik,t}\approx$
$\displaystyle\sum_{s_{k,t}}\tanh\theta_{i,t}^{*}(s_{k,t})s_{k,t}Q^{*}(s_{k,t})-m_{i,t}m_{k,t}.$
(55)

Note that the approximation of equal-time correlations may not be symmetric
for $C_{ik,t}$ and $C_{ki,t}$. In the results of this paper we use the average
of the two.

Figure 3: Forward Ising problem. Top: Evolution of average activation rates
(magnetizations) (A), equal-time correlations (B) and delayed correlations (C)
found by different mean-field methods for $\beta=\beta_{c}$. Middle:
Comparison of the activation rates (D), equal-time correlations (E) and
delayed correlations (F) found by the different Plefka approximations
(ordinate, $p$ superscript) with the original values (abscissa, $o$
superscript) for $\beta=\beta_{c}$ and $t=128$. Black lines represent the
identity line. Bottom: Average squared error of the magnetizations
$\epsilon_{\mathbf{m}}=\langle\langle(m^{o}_{i,t}-m^{p}_{i,t})^{2}\rangle_{i}\rangle_{t}$
(G), equal-time correlations
$\epsilon_{\mathbf{C}}=\langle\langle(C^{o}_{ik,t}-C^{p}_{ik,t})^{2}\rangle_{ik}\rangle_{t}$
(H), and delayed correlations
$\epsilon_{\mathbf{D}}=\langle\langle(D^{o}_{ik,t}-D^{p}_{ik,t})^{2}\rangle_{il}\rangle_{t}$
(I) for 21 values of $\beta$ in the range $[0.7\beta_{c},1.3\beta_{c}]$.

### Comparison of the different approximations

In the subsequent sections, we compare the family of Plefka approximation
methods described above by testing their performance in the forward and
inverse Ising problems. More specifically, we compare the second order
approximations of Plefka[$t-1,t$] and Plefka[$t$], the first order
approximation of Plefka[$t-1$], and the second order pairwise approximation of
Plefka2[$t$]. We define an Ising model as an asymmetric version of the kinetic
Sherrington-Kirkpatrick (SK) model, setting its parameters around the
equivalent of a ferromagnetic phase transition in the equilibrium SK model.
External fields $H_{i}$ are sampled from independent uniform distributions
$\mathcal{U}(-\beta H_{0},\beta H_{0})$, $H_{0}=0.5$, whereas coupling terms
$J_{ij}$ are sampled from independent Gaussian distributions
$\mathcal{N}(\beta\frac{J_{0}}{N},\beta^{2}\frac{J_{\sigma}^{2}}{N})$,
$J_{0}=1,J_{\sigma}=0.1$, where $\beta$ is a scaling parameter (i.e., an
inverse temperature).

Generally, mean-field methods are suitable for approximating properties of
systems with small fluctuations. However, there is evidence that many
biological systems operate in critical, highly fluctuating regimes [6, 5]. In
order to examine different approximations in such a biologically plausible yet
challenging situation, we select the model parameters around a phase
transition point displaying large fluctuations.

To find such conditions, we employed path integral methods to solve the
asymmetric SK model (Supplementary Note 6). We find that the stationary
solution of the asymmetric model displays for our choice of parameters a non-
equilibrium analogue of a critical point for a ferromagnetic phase transition,
which takes place at $\beta_{c}\approx 1.1108$ in thermodynamic limit (see
Supplementary Note 6, Supplementary Fig. 1). Note that the uniformly
distributed bias terms $\mathbf{H}$ shifts the phase transition point from
$\beta=1$ obtained at $\mathbf{H}=\mathbf{0}$. By simulation of the finite
size system, we confirmed that the maximum fluctuations in the model are found
near the theoretical $\beta_{c}$, which shows maximal covariance values (see
Supplementary Note 6, Supplementary Fig. 2).

Fluctuations of a system are generally expected to be maximized at a critical
phase transition [19]. In addition, entropy production (a signature of time
irreversibility) has been suggested as an indicator of phase transitions. For
example, it presents a peak at the transition point of a continuous phase
transition in a non-equilibrium Curie-Weiss Ising model with oscillatory field
[50] and some instances of mean-field majority vote models [51, 52]. We found
that the entropy production of the kinetic Ising system is also maximized
around $\beta_{c}$ (discussed later, see also Methods for its derivation).

### Forward Ising problem

We examine the performance of the different Plefka expansions in predicting
the evolution of an asymmetric SK model of size $N=512$ with random
$\mathbf{H}$ and $\mathbf{J}$. To study the non-stationary transient dynamics
of the model, we start from $\mathbf{s}_{0}=\mathbf{1}$ (all elements set to 1
at $t=0$) and recursively update its state for $T=128$ steps. We repeated this
stochastic simulation for $R=10^{6}$ trials for $21$ values of $\beta$ in the
range $[0.7\beta_{c},1.3\beta_{c}]$, except for the reconstruction of the
phase transition where we used $R=10^{5}$ and $201$ values of $\beta$ in the
same range. Using the $R$ samples, we computed the statistical moments and
cumulants of the system, $\mathbf{m}_{t}$, $\mathbf{C}_{t}$, and
$\mathbf{D}_{t}$ at each time step. We then computed their averages over the
system units, i.e., $\langle m_{i,t}\rangle_{i}$, $\langle
C_{ik,t}\rangle_{ik}$ and $\langle D_{il,t}\rangle_{il}$, where the angle
bracket denotes average over indices of its subscript.

The black solid lines in Fig. 3A,B,C display non-stationary dynamics of these
averaged statistics from $t=0,\dots,128$, simulated by the original model at
$\beta=\beta_{c}$. In comparison, color lines display these statistics
predicted by the family of Plefka approximations that are recursively computed
using the obtained equations, starting from the initial state
$\mathbf{m}_{0}=\mathbf{1}$, $\mathbf{C}_{0}=\mathbf{0}$ and
$\mathbf{D}_{0}=\mathbf{0}$. We observe that although the recursive
application of all the approximation methods provides good predictions for the
transient dynamics of the mean activation rates $\mathbf{m}_{t}$ until its
convergence (Fig. 3A), the predictions using Plefka[$t$] and especially the
proposed Plefka2[$t$] approximations are closer to the true dynamics than the
others. Evolution of the mean equal-time and time-delayed correlations
$\mathbf{C}_{t},\mathbf{D}_{t}$ is precisely captured only by our new method
Plefka2[$t$]. In contrast, Plefka[$t$] overestimates correlations while
Plefka[$t-1$] and Plefka[$t-1,t$] underestimate correlations.

Performance of the methods in predicting individual activation rates and
correlations are displayed in Fig. 3D, E, F by comparing vectors
$\mathbf{m}_{t}$, $\mathbf{C}_{t}$ and $\mathbf{D}_{t}$ at the last time step
($t=128$) of the original model ($o$ superscript) and those of the Plefka
approximations ($p$ superscript). For activation rates $\mathbf{m}_{t}$, the
proposed Plefka2[$t$] and Plefka[$t$] perform slightly better than the others
(see also Fig. 3A). While being overestimated by Plefka[$t$], underestimated
moderately by Plefka[$t-1$] and significantly by Plefka[$t-1,t$], equal-time
and time-delayed correlations $\mathbf{C}_{t},\mathbf{D}_{t}$ are best
predicted by Plefka2[$t$] (Fig. 3E,F).

The above results are obtained at the critical $\beta=\beta_{c}$, intuitively
the most challenging point for mean-field approximations. In order to further
show that our novel approximation Plefka2[$t$] systematically outperforms the
others in a wider parameter range, we repeated the analysis for different
inverse temperatures $\beta$ (the same random parameters are applied for all
$\beta$). Fig. 3G,H,U respectively show the averaged squared errors (averaged
over time and units) of the activation rates $\epsilon_{\mathbf{m}}$, equal-
time correlations $\epsilon_{\mathbf{C}}$ and delayed correlations
$\epsilon_{\mathbf{D}}$ between the original model and approximations,
averaged over units and time for values of $\beta$ in the range
$[0.7\beta_{c},1.3\beta_{c}]$. Figs. 3G,H,I show that Plefka2[$t$] outperforms
the other methods in computing $\mathbf{m}_{t},\mathbf{C}_{t},\mathbf{D}_{t}$
(with the exception of a certain region of $\beta>\beta_{c}$ in which
Plefka[$t$] is slightly better), yielding consistently a low error bound for
all values of $\beta$. Errors of these approximations are smaller when the
system is away from $\beta_{c}$.

Figure 4: Inverse Ising problem. Top: Inferred external fields (A) and
couplings (B) found by different mean-field models plotted versus the real
ones for $\beta=\beta_{c}$. Black lines represent the identity line. Bottom:
Average squared error of inferred external fields
$\epsilon_{\mathbf{H}}=\langle(H^{o}_{i}-H^{p}_{i})^{2}\rangle_{i}$ (C) and
couplings
$\epsilon_{\mathbf{J}}=\langle(J^{o}_{ij}-J^{p}_{ij})^{2}\rangle_{ij}$ (D) for
21 values of $\beta$ in the range $[0.7\beta_{c},1.3\beta_{c}]$.

### Inverse Ising problem

We apply the approximation methods to the inverse Ising problem by using the
data generated above for the trajectory of $T=128$ steps and $R=10^{6}$ trials
to infer the parameters of the model, $\mathbf{H}$ and $\mathbf{J}$. The model
parameters are estimated by the Boltzmann learning method under the maximum
likelihood principle: $\mathbf{H}$ and $\mathbf{J}$ are updated to minimize
the differences between the average rates $\mathbf{m}_{t}$ or delayed
correlations $\mathbf{D}_{t}$ of the original data and the model
approximations, which can significantly reduce computational time (see
Methods). While Boltzmann learning requires to compute the likelihood of every
point in a trajectory and every trial ($RT$ calculations) each iteration, we
can estimate the gradient at each iteration in a one-shot computation by
applying the Plefka approximations (Methods). At $\beta=\beta_{c}$ (Fig.
4A,B), we observe that the classical Plefka[$t-1,t$] approximation adds
significant offset values to the fields $\mathbf{H}$ and couplings
$\mathbf{J}$. In contrast, Plefka[$t$], Plefka[$t-1$] and Plefka2[$t$] are all
precise in estimating the values of $\mathbf{H}$ and $\mathbf{J}$.

Fig. 4C and D show the mean squared error $\epsilon_{\mathbf{H}}$,
$\epsilon_{\mathbf{J}}$ for bias terms and couplings between the original
model and the inferred values for different $\beta$. In this case, errors are
large in the estimation of $\mathbf{J}$ for Plefka[$t-1,t$]. In comparison,
Plefka[$t$], Plefka[$t-1$] and Plefka2[$t$] work equally well even in the high
fluctuation regime ($\beta\approx\beta_{c}$). Since the inverse Ising problem
is solved by applying approximation one single time step (per iteration), it
is not as challenging as the forward problem that can accumulate errors by
recursively applying the approximations. Therefore, different approximations
other than the classical mean-field Plefka[$t-1,t$] perform equally well in
this case.

### Phase transition reconstruction

We have shown how different methods perform in computing the behavior of the
system (forward problem) and inferring the parameters of a given network from
its activation data (inverse problem). Combining the two, we can ask how well
the methods explored here can reconstruct the behavior of a system from data,
potentially exploring behaviors under different conditions than the recorded
data.

First, in Fig. 5.A-C we examine how the different approximation methods
approximate fluctuations (equal-time and time-delayed covariances) and the
entropy production (see Methods) at $t=128$ after solving the forward problem
by recursively applying the approximations for the 128 steps. As we mentioned
above, the asymmetric SK model explored here presents maximum fluctuations and
maximum entropy production around $\beta=\beta_{c}$ (Supplementary Note 6,
Supplementary Fig. 2). However, we see that Plefka[$t-1,t$] and Plefka[$t-1$]
cannot reproduce the behavior of correlations $\mathbf{C}_{t}$ and
$\mathbf{D}_{t}$ of the original SK model around the transition point.
Plefka[$t$] and Plefka2[$t$] show much better performance in capturing the
behavior of $\mathbf{C}_{t}$ and $\mathbf{D}_{t}$ in the phase transition,
although Plefka[$t$] overestimates both correlations. Additionally, all the
methods capture the phase transition in entropy production, though Plefka[$t$]
overestimates its value around $\beta_{c}$ and Plefka2[$t$] is more precise
than the other methods.

Next, we combine the forward and inverse Ising problem and try to reproduce
the transition in the asymmetric SK model in the models inferred from the
data. We first obtain the values of $\mathbf{H},\mathbf{J}$ by solving the
inverse problem from the data sampled at $\beta=\beta_{c}$, and next we solve
again the forward problem with those estimated parameters re-scaled by a new
inverse temperature $\tilde{\beta}$. The results for the correlations (Fig.
5.D-E) show that in this case Plefka[$t-1,t$] works badly, not being able to
capture the transition. Plefka[$t-1$] shows similar performances as in the
forward problem, and Plefka[$t$] and Plefka2[$t$] have a similar behavior,
underestimating fluctuations slightly. When we analyze entropy production of
the system (Fig. 5.F), we find that Plefka2[$t$] exhibits better performance
with a high precision, with Plefka[$t-1$] slightly overestimating it,
Plefka[$t$] underestimating it, and Plefka[$t-1,t$] not capturing the phase
transition. Overall, the results above suggest that Plefka2[$t$] is better
suited to identify non-equilibrium phase transitions in models reconstructed
from experimental data.

Figure 5: Reconstructing phase transition of kinetic Ising systems. Top:
Average of the Ising model’s equal-time correlations (A), delayed correlations
(B), and entropy production (shown as an exponential for better presentation
of its maximum) (C), at the last step $t=128$ found by different mean-field
methods for $\beta=\beta_{c}$. Bottom (D, E, F): The same as above using the
reconstructed network $\mathbf{H},\mathbf{J}$ by solving the inverse Ising
problem at $\beta=\beta_{c}$ and multiplying a fictitious inverse temperature
$\tilde{\beta}$ to the estimated parameters. The stars are marked at the
values of $\tilde{\beta}$ that yield maximum fluctuations or maximum entropy
production.

## Discussion

We have proposed a framework that unifies different mean-field approximations
of the evolving statistical properties of non-equilibrium Ising models. This
allows us to derive approximations premised on specific assumptions about the
correlation structure of the system previously proposed in the literature.
Furthermore, using our framework we derive a new approximation (Plefka2[$t$])
using atypical assumptions for mean-field methods, i.e. the maintenance of
pairwise correlations in the system. This new pairwise approximation
outperforms existing ones for approximating the behavior of an asymmetric SK
model near the non-equilibrium equivalent of a ferromagnetic phase transition
(see Supplementary Note 6), where classical mean-field approximations face
problems. This shows that the proposed methods are useful tools to analyze
large-scale, non-equilibrium dynamics near critical regimes expected for
biological and social systems. However, we note that low-temperature spin
phases (e.g. the spin-glass phase in symmetric models) also impose limitations
on mean-field approximations [32, 41], which could be further explored with
methods like the ones presented here.

The generality of this framework allows us to picture other approximations
with atypical assumptions. For example, the Sessak-Monasson expansion [53] for
an equilibrium Ising model assumes a linear relation between $\alpha$ and spin
correlations. An equivalent equilibrium expansion could use an effective field
$h(\alpha)$ nonlinearly dependent on $\alpha$, satisfying linear
$\mathbf{C}_{t}(\alpha)=\alpha\mathbf{C}_{t}$ or
$\mathbf{D}_{t}(\alpha)=\alpha\mathbf{D}_{t}$ relations. As another extension,
Plefka2[$t$] could incorporate higher-order interactions. As Eq. 43, 53 are
each equivalent to two mean-field approximations with $s_{l,t-1}=\pm 1$
respectively, a generalized PlefkaM[$t$] would involve $2^{M-1}$ equations,
increasing accuracy but also computational costs. In general, reference models
$Q(\mathbf{s}_{t})$ set coupling parameters of the model to zero at some steps
of its dynamics. Other parameters (e.g. fields) are either free parameters
fitted as m-projection from $P(\mathbf{s}_{t})$, or preserved to their
original value (see Supplementary Note 7 for comparing free and fixed
parameters of each model). Augmenting accuracy by increasing parameters often
involves a computational cost. As a practical guideline for using each method,
Supplementary Note 7 compares their precision and computation time in the
forward and inverse problems (see also Supplementary Fig. 3,4).

Asides from its theoretical implications, our unified framework offers
analysis tools for diverse data-driven research fields. In neuroscience, it
has been popular to study the activity of ensembles of neurons by inferring an
equilibrium Ising model with homogeneous (fixed) parameters [23] or
inhomogeneous (time-dependent) parameters [54, 25] from empirical data.
Extended analyses based on the equilibrium model have reported that neurons
operate near a critical regime [6, 5]. However, studies of non-equilibrium
dynamics in neural spike trains are scarce [26, 7, 55], partly due to the lack
of systematic methods for analysing large-scale non-equilibrium data from
neurons exhibiting large fluctuations. The proposed pairwise model
Plefka2[$t$] is suitable for simulating such network activities, being more
accurate than previous methods in predicting the network evolution at
criticality (Fig. 3) and in testing if the system is near the maximally
fluctuating regime (Fig. 5). In particular, application of our methods for
computing entropy production in non-equilibrium systems could provide tools
for characterizing the non-equilibrium dynamics of neural systems [56].

In summary, a unified framework of mean-field theories offers a systematic way
to construct suitable mean-field methods in accordance with the statistical
properties of the systems researchers wish to uncover. This is expected to
foster a variety of tools to analyze large-scale non-equilibrium systems in
physical, biological, and social systems.

## Methods

### Boltzmann learning in the inverse Ising problem

Let $\mathbf{S}_{t}^{r}=\\{S_{1,t}^{r},S_{2,t}^{r},\ldots,S_{N,t}^{r}\\}$ for
$t=1,\dots,T$ be observed states of a process described by Eq. 1 at the $r$-th
trial ($r=1,\dots,R$). We also define $\mathbf{S}_{1:T}$ to represent the
processes from all trials. The inverse Ising problem consists in inferring the
external fields $\mathbf{H}$ and couplings $\mathbf{J}$ of the system. These
parameters can be estimated by maximizing the log-likelihood
$l(\mathbf{S}_{1:T})$ of the observed states under the model:

$\displaystyle l(\mathbf{S}_{1:T})=$
$\displaystyle\log\prod_{t=1}^{T}\prod_{r=1}^{R}P(\mathbf{S}_{t}^{r}|\mathbf{S}_{t-1}^{r})$
$\displaystyle=$
$\displaystyle\sum_{t}\sum_{r}\sum_{i}\left(S_{i,t}^{r}h_{i,t}^{r}-\log 2\cosh
h_{i,t}^{r}\right),$ (56)

with $h_{i,t}^{r}=H_{i}+\sum_{j}J_{ij}S_{j,t-1}^{r}$. The learning steps are
obtained as:

$\displaystyle\frac{\partial l(\mathbf{S}_{1:T})}{\partial H_{i}}=$
$\displaystyle RT(\langle S_{i,t}^{r}\rangle_{r,t}-\langle\tanh
h_{i,t}^{r}\rangle_{r,t}),$ (57) $\displaystyle\frac{\partial
l(\mathbf{S}_{1:T})}{\partial J_{il}}=$ $\displaystyle RT(\langle
S_{i,t}^{r}S_{l,t-1}^{r}\rangle_{r,t}-\langle\tanh
h_{i,t}^{r}S_{l,t-1}^{r,t}\rangle_{r,t}),$ (58)

where $\langle\cdot\rangle_{r}$ denotes average over trials. We solve the
inverse Ising problem by applying these equations as a gradient ascent rule
adjusting $\mathbf{H}$ and $\mathbf{J}$. The second terms of Eqs. 57, 58 need
to be computed at every iteration, thus the computational cost grows linearly
with $R\times T$. However, the use of mean-field approximations can
significantly reduce the cost when a large number of samples $R$ and time bins
$T$ are used to correctly estimate activation rates and correlations in large
networks. Here the second terms can be written as

$\displaystyle\langle\tanh h_{i,t}^{r}\rangle_{r,t}=$
$\displaystyle\sum_{\mathbf{s},\tilde{\mathbf{s}}}s_{i}P(\mathbf{s}|\tilde{\mathbf{s}})\overline{P}(\tilde{\mathbf{s}})=m_{i},$
(59) $\displaystyle\langle\tanh h_{i,t}^{r}S_{l,t-1}^{r}\rangle_{r,t}=$
$\displaystyle\sum_{\mathbf{s},\tilde{\mathbf{s}}}s_{i}\tilde{s}_{l}P(\mathbf{s}|\tilde{\mathbf{s}})\overline{P}(\tilde{\mathbf{s}})$
$\displaystyle=$ $\displaystyle D_{il}+m_{i}\tilde{m}_{l},$ (60)

where
$\overline{P}(\tilde{\mathbf{s}})=\frac{1}{RT}\sum_{r,t}\delta(\tilde{\mathbf{s}},\mathbf{S}_{t}^{r})$
is the empirical distribution averaged over trials and trajectories (with
$\delta$ being a Kronecker delta) and $\tilde{m}_{l}$ is the average
activation rate computed from the empirical distribution.
$P(\mathbf{s}|\tilde{\mathbf{s}})$ is defined as Eq. 1. We then approximate
$m_{i}$ and $D_{il}$ using the mean-field equations. Note that when we apply
the mean-field equations, we replaced all statistics related to the previous
step with those computed by the empirical distribution. By applying the mean-
field methods, we reduced the computation of $R$ trials of trajectories of
length $T$ into a single computation (instead of $RT$ calculations). In our
numerical tests, gradient ascent was executed using learning coefficients
$\eta_{H}=0.1/RT,\eta_{J}=1/(RT\sqrt{N})$, starting from
$\mathbf{H}=\mathbf{0},\mathbf{J}=\mathbf{0}$.

### Entropy production of the kinetic Ising model

The entropy production is defined as the KL divergence between the forward and
backward path, quantifying the irreversibility of the system [57, 17, 55]:

$\displaystyle\langle\sigma_{t}\rangle=$
$\displaystyle\frac{1}{2}\sum_{\mathbf{s}_{t},\mathbf{s}_{t-1}}\big{(}P(\mathbf{s}_{t-1})P(\mathbf{s}_{t}|\mathbf{s}_{t-1})-P(\mathbf{s}_{t})P_{\rm
B}(\mathbf{s}_{t-1}|\mathbf{s}_{t})\big{)}$
$\displaystyle\cdot\log\frac{P(\mathbf{s}_{t}|\mathbf{s}_{t-1})P(\mathbf{s}_{t-1})}{P_{\rm
B}(\mathbf{s}_{t-1}|\mathbf{s}_{t})P(\mathbf{s}_{t})}.$ (61)

where $P_{\rm B}(\mathbf{s}_{t-1}|\mathbf{s}_{t})$ is a probability of the
backward trajectory defined as in Eq. 1 but switching $\mathbf{s}_{t}$ and
$\mathbf{s}_{t-1}$. Assuming a non-equilibrium steady state, where
$P(\mathbf{s}_{t})=P(\mathbf{s}_{t-1})$, the entropy production of the kinetic
Ising system is computed as:

$\displaystyle\langle\sigma_{t}\rangle=$
$\displaystyle\sum_{\mathbf{s}_{t},\mathbf{s}_{t-1}}P(\mathbf{s}_{t},\mathbf{s}_{t-1})\Big{(}\sum_{i}H_{i}(s_{i,t}-s_{i,t-1})$
$\displaystyle+\sum_{ij}J_{ij}(s_{i,t}s_{j,t-1}-s_{i,t-1}s_{j,t})$
$\displaystyle+\log P(\mathbf{s}_{t-1})-\log P(\mathbf{s}_{t})$
$\displaystyle-\log(2\cosh(H_{i}+\sum_{j}J_{ij}s_{j,t-1}))$
$\displaystyle+\log(2\cosh(H_{i}+\sum_{j}J_{ij}s_{j,t}))\Big{)}$
$\displaystyle=$
$\displaystyle\sum_{ij}J_{ij}(D_{ij,t}-D_{ji,t})=\sum_{ij}(J_{ij}-J_{ji})D_{ij,t}.$
(62)

###### Acknowledgements.

The authors thank Yasser Roudi and Masanao Igarashi for valuable comments and
discussions on this manuscript. This work was supported in part by the
Cooperative Intelligence Joint Research between Kyoto University and Honda
Research Institute Japan, MEXT/JSPS KAKENHI Grant Number JP 20K11709, and the
grant of Joint Research by the National Institutes of Natural Sciences (NINS
Program No. 01112005). M.A. was funded by the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 892715 and the University of the Basque Country UPV/EHU post-
doctoral training program grant ESPDOC17/17, and supported in part by the
Basque Government project IT 1228-19 and project Outonomy
(PID2019-104576GB-I00) by the Spanish Ministry of Science and Innovation.

## Code Availability

The source code for implementing the results in this work is available under
GPL license at GitHub https://github.com/MiguelAguilera/kinetic-Plefka-
expansions [58] (DOI:10.5281/zenodo.4357634).

## Data Availability

The datasets generated and analysed during the current study are available
under CC BY license at Zenodo https://zenodo.org/record/4318983 [59]
(DOI:10.5281/zenodo.4318983).

## Competing Interests

The authors declare no competing interests.

## Author contributions

M.A., S.A.M. and H.S. designed and reviewed research; M.A. contributed
analytical and numerical results; M.A., S.A.M. and H.S. wrote the paper.

## Additional information

### Supplementary information

The online version contains supplementary material available at
https://doi.org/10.1038/s41467-021-20890-5.

## References

  * Nguyen  _et al._ [2016] J. P. Nguyen, F. B. Shipley, A. N. Linder, G. S. Plummer, M. Liu, S. U. Setru, J. W. Shaevitz, and A. M. Leifer, Whole-brain calcium imaging with cellular resolution in freely behaving Caenorhabditis elegans, Proceedings of the National Academy of Sciences 113, E1074 (2016).
  * Ahrens  _et al._ [2013] M. B. Ahrens, M. B. Orger, D. N. Robson, J. M. Li, and P. J. Keller, Whole-brain functional imaging at cellular resolution using light-sheet microscopy, Nature Methods 10, 413 (2013).
  * Stringer  _et al._ [2019] C. Stringer, M. Pachitariu, N. Steinmetz, M. Carandini, and K. D. Harris, High-dimensional geometry of population responses in visual cortex, Nature , 1 (2019).
  * Nicolis and Prigogine [1977] G. Nicolis and I. Prigogine,  _Self-Organization in Nonequilibrium Systems: From Dissipative Structures to Order through Fluctuations_ , 1st ed. (Wiley, New York, 1977).
  * Tkačik  _et al._ [2015] G. Tkačik, T. Mora, O. Marre, D. Amodei, S. E. Palmer, M. J. Berry, and W. Bialek, Thermodynamics and signatures of criticality in a network of neurons, Proceedings of the National Academy of Sciences 112, 11508 (2015).
  * Mora  _et al._ [2015] T. Mora, S. Deny, and O. Marre, Dynamical Criticality in the Collective Activity of a Population of Retinal Neurons, Physical Review Letters 114, 078105 (2015).
  * Hertz  _et al._ [2013] J. Hertz, Y. Roudi, and J. Tyrcha, Ising model for inferring network structure from spike data: In Principal of Neural Coding, in  _Principles of Neural Coding_ (CRC Press, 2013) pp. 527–546.
  * Roudi  _et al._ [2015] Y. Roudi, B. Dunn, and J. Hertz, Multi-neuronal activity and functional connectivity in cell assemblies, Current opinion in neurobiology 32, 38 (2015).
  * Bouchaud [2013] J.-P. Bouchaud, Crises and collective socio-economic phenomena: simple models and challenges, Journal of Statistical Physics 151, 567 (2013).
  * Evans  _et al._ [1993] D. J. Evans, E. G. D. Cohen, and G. P. Morriss, Probability of second law violations in shearing steady states, Physical review letters 71, 2401 (1993).
  * Jarzynski [1997] C. Jarzynski, Nonequilibrium equality for free energy differences, Physical Review Letters 78, 2690 (1997).
  * Crooks [1998] G. E. Crooks, Nonequilibrium measurements of free energy differences for microscopically reversible Markovian systems, Journal of Statistical Physics 90, 1481 (1998).
  * Crooks [1999] G. E. Crooks, Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences, Physical Review E 60, 2721 (1999).
  * Lebowitz and Spohn [1999] J. L. Lebowitz and H. Spohn, A Gallavotti–Cohen-type symmetry in the large deviation functional for stochastic dynamics, Journal of Statistical Physics 95, 333 (1999).
  * Ito  _et al._ [2020] S. Ito, M. Oizumi, and S.-i. Amari, Unified framework for the entropy production and the stochastic interaction based on information geometry, Physical Review Research 2, 033048 (2020).
  * Evans and Searles [2002] D. J. Evans and D. J. Searles, The fluctuation theorem, Advances in Physics 51, 1529 (2002).
  * Seifert [2012] U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Reports on progress in physics 75, 126001 (2012).
  * Gaspard [2007] P. Gaspard, Time Asymmetry in Nonequilibrium Statistical Mechanics, in  _Special Volume in Memory of Ilya Prigogine_ (John Wiley & Sons, Ltd, 2007) pp. 83–133.
  * Salinas [2001] S. R. A. Salinas, The Ising Model, in  _Introduction to Statistical Physics_, Graduate Texts in Contemporary Physics, edited by S. R. A. Salinas (Springer New York, New York, NY, 2001) pp. 257–276.
  * Ackley  _et al._ [1985] D. H. Ackley, G. E. Hinton, and T. J. Sejnowski, A learning algorithm for Boltzmann machines, Cognitive science 9, 147 (1985).
  * Witoelar and Roudi [2011] A. Witoelar and Y. Roudi, Neural network reconstruction using kinetic Ising models with memory, BMC Neuroscience 12, P274 (2011).
  * Donner and Opper [2017] C. Donner and M. Opper, Inverse Ising problem in continuous time: A latent variable approach, Physical Review E 96, 062104 (2017).
  * Schneidman  _et al._ [2006] E. Schneidman, M. J. Berry, R. Segev, and W. Bialek, Weak pairwise correlations imply strongly correlated network states in a neural population, Nature 440, 1007 (2006).
  * Cocco  _et al._ [2009] S. Cocco, S. Leibler, and R. Monasson, Neuronal couplings between retinal ganglion cells inferred by efficient inverse statistical physics methods, Proceedings of the National Academy of Sciences 106, 14058 (2009).
  * Shimazaki  _et al._ [2012] H. Shimazaki, S.-i. Amari, E. N. Brown, and S. Grün, State-Space Analysis of Time-Varying Higher-Order Spike Correlation for Multiple Neural Spike Train Data, PLOS Computational Biology 8, e1002385 (2012).
  * Tyrcha  _et al._ [2013] J. Tyrcha, Y. Roudi, M. Marsili, and J. Hertz, The effect of nonstationarity on models inferred from neural data, Journal of Statistical Mechanics: Theory and Experiment 2013, P03005 (2013).
  * Thouless  _et al._ [1977] D. J. Thouless, P. W. Anderson, and R. G. Palmer, Solution of ’Solvable model of a spin glass’, The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics 35, 593 (1977).
  * Kappen and Rodríguez [1998] H. J. Kappen and F. B. Rodríguez, Efficient Learning in Boltzmann Machines Using Linear Response Theory, Neural Computation 10, 1137 (1998).
  * Roudi  _et al._ [2009a] Y. Roudi, E. Aurell, and J. A. Hertz, Statistical Physics of Pairwise Probability Models, Frontiers in Computational Neuroscience 3, 10.3389/neuro.10.022.2009 (2009a).
  * Roudi  _et al._ [2009b] Y. Roudi, J. Tyrcha, and J. Hertz, Ising model for neural data: Model quality and approximate methods for extracting functional connectivity, Physical Review E 79, 051915 (2009b).
  * Donner  _et al._ [2017] C. Donner, K. Obermayer, and H. Shimazaki, Approximate Inference for Time-Varying Interactions and Macroscopic Dynamics of Neural Populations, PLOS Computational Biology 13, e1005309 (2017).
  * Plefka [1982] T. Plefka, Convergence condition of the TAP equation for the infinite-ranged Ising spin glass model, Journal of Physics A: Mathematical and general 15, 1971 (1982).
  * Tanaka [1998] T. Tanaka, Mean-field theory of Boltzmann machine learning, Physical Review E 58, 2302 (1998).
  * Tanaka [1999] T. Tanaka, A theory of mean field approximation, in  _Advances in Neural Information Processing Systems_ (1999) pp. 351–360.
  * Bhattacharyya and Keerthi [2000] C. Bhattacharyya and S. S. Keerthi, Information geometry and Plefka’s mean-field theory, Journal of Physics A: Mathematical and General 33, 1307 (2000).
  * Tanaka [2001] T. Tanaka, Information Geometry of Mean-Field Approximation, in  _Advanced Mean Field Methods: Theory and Practice_ (MIT Press, 2001).
  * Amari  _et al._ [2001] S. Amari, S. Ikeda, and H. Shimokawa, Information geometry of $\alpha$-projection in mean-field approximation, in  _Advanced Mean Field Methods: Theory and Practice_ (MIT Press, 2001).
  * Kappen and Spanjers [2000] H. J. Kappen and J. J. Spanjers, Mean field theory for asymmetric neural networks, Physical Review E 61, 5658 (2000).
  * Roudi and Hertz [2011a] Y. Roudi and J. Hertz, Dynamical TAP equations for non-equilibrium Ising spin glasses, Journal of Statistical Mechanics: Theory and Experiment 2011, P03031 (2011a).
  * Roudi and Hertz [2011b] Y. Roudi and J. Hertz, Mean Field Theory for Nonequilibrium Network Reconstruction, Physical Review Letters 106, 048702 (2011b).
  * Mézard and Sakellariou [2011] M. Mézard and J. Sakellariou, Exact mean-field inference in asymmetric kinetic Ising systems, Journal of Statistical Mechanics: Theory and Experiment 2011, L07001 (2011).
  * Mahmoudi and Saad [2014] H. Mahmoudi and D. Saad, Generalized mean field approximation for parallel dynamics of the Ising model, Journal of Statistical Mechanics: Theory and Experiment 2014, P07001 (2014).
  * Bachschmid-Romano  _et al._ [2016] L. Bachschmid-Romano, C. Battistin, M. Opper, and Y. Roudi, Variational perturbation and extended Plefka approaches to dynamics on random networks: the case of the kinetic Ising model, Journal of Physics A: Mathematical and Theoretical 49, 434003 (2016).
  * Pressé  _et al._ [2013] S. Pressé, K. Ghosh, J. Lee, and K. A. Dill, Principles of maximum entropy and maximum caliber in statistical physics, Reviews of Modern Physics 85, 1115 (2013).
  * Amari and Nagaoka [2007] S. Amari and H. Nagaoka,  _Methods of information geometry_ , Vol. 191 (American Mathematical Soc., 2007).
  * Amari [2016] S. Amari,  _Information geometry and its applications_ , Vol. 194 (Springer, 2016).
  * Amari  _et al._ [1992] S. Amari, K. Kurata, and H. Nagaoka, Information geometry of Boltzmann machines, IEEE Transactions on neural networks 3, 260 (1992).
  * Oizumi  _et al._ [2016] M. Oizumi, N. Tsuchiya, and S.-i. Amari, Unified framework for information integration based on information geometry, Proceedings of the National Academy of Sciences 113, 14817 (2016).
  * Saul and Jordan [1996] L. K. Saul and M. I. Jordan, Exploiting tractable substructures in intractable networks, in  _Advances in Neural Information Processing Systems_ (1996) pp. 486–492.
  * Zhang and Barato [2016] Y. Zhang and A. C. Barato, Critical behavior of entropy production and learning rate: Ising model with an oscillating field, Journal of Statistical Mechanics: Theory and Experiment 2016, 113207 (2016).
  * Crochik and Tomé [2005] L. Crochik and T. Tomé, Entropy production in the majority-vote model, Physical Review E 72, 057103 (2005).
  * Noa  _et al._ [2019] C. F. Noa, P. E. Harunari, M. de Oliveira, and C. Fiore, Entropy production as a tool for characterizing nonequilibrium phase transitions, Physical Review E 100, 012104 (2019).
  * Sessak and Monasson [2009] V. Sessak and R. Monasson, Small-correlation expansions for the inverse Ising problem, Journal of Physics A: Mathematical and Theoretical 42, 055001 (2009).
  * Granot-Atedgi  _et al._ [2013] E. Granot-Atedgi, G. Tkačik, R. Segev, and E. Schneidman, Stimulus-dependent Maximum Entropy Models of Neural Population Codes, PLOS Computational Biology 9, e1002922 (2013).
  * Cofré  _et al._ [2019] R. Cofré, L. Videla, and F. Rosas, An Introduction to the Non-Equilibrium Steady States of Maximum Entropy Spike Trains, Entropy 21, 884 (2019).
  * Lynn  _et al._ [2020] C. W. Lynn, E. J. Cornblath, L. Papadopoulos, M. A. Bertolero, and D. S. Bassett, Non-equilibrium dynamics and entropy production in the human brain, arXiv preprint arXiv:2005.02526 (2020).
  * Schnakenberg [1976] J. Schnakenberg, Network theory of microscopic and macroscopic behavior of master equation systems, Reviews of Modern Physics 48, 571 (1976).
  * Aguilera [2020a] M. Aguilera, A unifying framework for mean field theories of asymmetric kinetic Ising systems [Code], GitHub 10.5281/zenodo.4357634 (2020a).
  * Aguilera [2020b] M. Aguilera, A unifying framework for mean field theories of asymmetric kinetic Ising systems [Dataset], Zenodo 10.5281/zenodo.4318983 (2020b).
  * Nishimori [2001] H. Nishimori,  _Statistical physics of spin glasses and information processing: an introduction_ (Clarendon Press, 2001).
  * Brunetti  _et al._ [1992] R. Brunetti, G. Parisi, and F. Ritort, Asymmetric Little spin-glass model, Physical Review B 46, 5339 (1992).
  * Lecomte  _et al._ [2007] V. Lecomte, C. Appert-Rolland, and F. van Wijland, Thermodynamic Formalism for Systems with Markov Dynamics, Journal of Statistical Physics 127, 51 (2007).

## Supplementary Notes

### List of changes respect to the published version

This version of the Supplementary Information contains the following
corrections respect to the published version available at
https://doi.org/10.1038/s41467-021-20890-5.

  * •

A $g_{i,u}$ term was removed from Eq. SI.6.331.

  * •

Missing $s_{j,v-1}$ and $\mathrm{i}\mkern 1.0mu\hat{\vartheta}_{i,v}$ terms
were added to Eq. SI.6.358.

  * •

Missing $\sum_{i}$ were added and extra $N$ terms removed in the set of
equations starting at Eq. SI.6.369

  * •

$\mathbf{z}$ was substituted by $\mathbf{\tilde{z}}$ in Eq. SI.6.391

  * •

A missing $z$ term was added to Eq. SI.6.397

  * •

Missing limit operator has been added to Eq. SI.6.400 and SI.6.402

  * •

An intermediate step was added for deriving the critical exponents at Eq.
SI.6.409

##  Supplementary Note 1 General approach for Plefka expansions

Let $\mathbf{s}_{t}=\\{s_{i,t}\\}_{i},~{}i=1,\dots,N$ be the state of the
system at time $t$ and
$\mathbf{s}_{1:T}=\\{\mathbf{s}_{i,t}\\}_{i,t},~{}i=1,\dots,N,~{}t=1,\dots,T$
a trajectory of the system. Given an initial state $\mathbf{s}_{0}$, the
probability of a trajectory $\mathbf{s}_{1:T}$ of a kinetic Ising model is:

$\displaystyle P(\mathbf{s}_{1:T})=$
$\displaystyle\prod_{t}P(\mathbf{s}_{t}|\mathbf{s}_{t-1})=\prod_{t}\exp\left(\sum_{i}s_{i,t}h_{i,t}-\psi\right),$
(SI.1.63) $\displaystyle h_{i,t}=$ $\displaystyle
H_{i}+\sum_{j}J_{ij}s_{j,t-1},$ (SI.1.64) $\displaystyle\psi=$
$\displaystyle\sum_{t,i}\log 2\cosh h_{i,t}.$ (SI.1.65)

Note that $\psi$ depends on the specific trajectory $\mathbf{s}_{0:T}$. The
manifold $\mathcal{P}=\\{P(\mathbf{s}_{1:T}|\mathbf{H},\mathbf{J})\\}$ defines
the family of probability distributions of the trajectories of all kinetic
Ising models. Within this manifold, we consider a submanifold
$\mathcal{P}_{0}$ in which the probability distributions of a set of elements
of the system $\\{\\{s_{i,t}\\}_{i\in\mathcal{I}_{t}}\\}_{t}$ is tractable and
independent from the rest of the elements of the system. Here
$\mathcal{I}_{t}$ denotes the indices of the tractable elements of the system
at time $t$. Trajectories in the submanifold are defined as:

$\displaystyle P_{0}(\mathbf{s}_{1:T})=$
$\displaystyle\prod_{t}P_{0}(\mathbf{s}_{t}|\mathbf{s}_{1:t-1})=\prod_{t}\exp\left(\sum_{i\in\mathcal{I}_{t}}s_{i,t}\theta_{i,t}+\sum_{i\in\overline{\mathcal{I}}_{t}}s_{i,t}h_{i,t}-\psi_{0}\right),$
(SI.1.66) $\displaystyle\psi_{0}=$
$\displaystyle\sum_{t,i\in\mathcal{I}_{t}}\log
2\cosh\theta_{i,t}+\sum_{t,i\in\overline{\mathcal{I}}_{t}}\log 2\cosh
h_{i,t},$ (SI.1.67)

where $\overline{\mathcal{I}}_{t}$ is the complement set of $\mathcal{I}_{t}$
for the elements at time $t$. The mean-field definition of the new effective
field is $\theta_{i,t}=\Theta_{i,t}$, although other definitions are possible.
We restrict the function $\theta_{i,t}$ to maximum caliber models composed of
individual fields and pairwise couplings

$\theta_{i,t}=\sum_{I\in\mathcal{C}_{i,t}}\Theta_{i,t,I}s_{I},$ (SI.1.68)

where $\mathcal{C}_{i,t}$ is a set of couplings
$\\{(j,\tau)\\},\tau<t,j\in\mathcal{I}_{\tau}$ between $s_{i,t}$ and other
nodes $s_{j,\tau}$ from the set of nodes with tractable properties. We define
the individual field of a unit with
$\Theta_{i,t,\emptyset}\equiv\Theta_{i,t}$, by defining $s_{\emptyset}=1$. In
general, the effect from the past spiking activities at $\tau<t$ can be
modeled by this equation, which includes the generalized linear model for
conditionally independent Bernoulli processes. In the approximations explored
in this paper, however, we focus on the effect from the immediate past
$\tau=t-1$.

Different approximations are defined through different definitions of
$\theta_{i,t}$ using a model connecting $P$ and $P_{0}$. This model is defined
through a parameter $\alpha$:

$\displaystyle P_{\alpha}(\mathbf{s}_{1:T})=$
$\displaystyle\prod_{t}P_{\alpha}(\mathbf{s}_{t}|\mathbf{s}_{1:t-\tau-1})=\prod_{t}\exp\left(\sum_{i\in\mathcal{I}_{t}}s_{i,t}((1-\alpha)\theta_{i,t}+\alpha
h_{i,t})+\sum_{i\in\overline{\mathcal{I}}_{t}}s_{i,t}h_{i,t}-\psi_{\alpha}\right),$
(SI.1.69) $\displaystyle\psi_{\alpha}=$
$\displaystyle\sum_{t,i\in\mathcal{I}_{t}}\log
2\cosh((1-\alpha)\theta_{i,t}+\alpha
h_{i,t}))+\sum_{t,i\in\overline{\mathcal{I}}_{t}}\log 2\cosh h_{i,t},$
(SI.1.70)

such that $P_{1}(\mathbf{s}_{1:T})=P(\mathbf{s}_{1:T})$.

The model $P_{0}$ that better approximates $P$ is the one that minimizes the
Kullback Leibler divergence:

$D_{KL}(P||P_{0})=\sum_{\mathbf{s}_{1:T}}P(\mathbf{s}_{1:T})\log\frac{P(\mathbf{s}_{1:T})}{P_{0}(\mathbf{s}_{1:T})}.$
(SI.1.71)

Thus its parameters $\Theta_{i,t,I}$ meet

$\displaystyle\frac{\partial D_{KL}(P||P_{0})}{\partial\Theta_{i,t,I}}=$
$\displaystyle-\sum_{\mathbf{s}_{1:T}}P(\mathbf{s}_{1:T})s_{i,t}s_{I}+\sum_{\mathbf{s}_{1:T}}P(\mathbf{s}_{1:T})\sum_{\bm{\sigma}_{t}}P_{0}(\bm{\sigma}_{t}|\mathbf{s}_{1:t-1})\sigma_{i,t}s_{I}$
(SI.1.72) $\displaystyle=$ $\displaystyle\langle
s_{i,t}s_{I}\rangle_{0}^{t}-\langle s_{i,t}s_{I}\rangle_{1}^{t}=0,$ (SI.1.73)

where $\langle\dots\rangle_{\alpha}^{t}=\sum_{\mathbf{s}_{1:T}}\dots
P_{\alpha}(\mathbf{s}_{t}|\mathbf{s}_{1:t-1})P(\mathbf{s}_{1:t-1})$. Thus, the
closest approximation $P_{0}$ to $P$ is the one in which $\langle
s_{i,t}s_{I}\rangle_{0}^{t}=\langle s_{i,t}s_{I}\rangle_{1}^{t},\forall i,I$.
Furthermore, models $P_{0}$ are tractable for indices $\\{\mathcal{I}_{t}\\}$,
in the sense that knowing $\Theta_{i,t,I}$, it is easy to compute $\langle
s_{i,t}s_{I}\rangle_{0}^{t}$. The value of parameters $\Theta_{i,t,I}$ cannot
be computed directly in general, but they can be approximated by computing a
Plefka expansion. Approximating $\langle s_{i,t}s_{I}\rangle_{1}^{t}$ by the
Taylor expansion of $\langle s_{i,t}s_{I}\rangle_{\alpha}^{t}$ from
$\alpha=0$, we have

$\langle s_{i,t}s_{I}\rangle_{\alpha}^{t}=\langle
s_{i,t}s_{I}\rangle_{\alpha=0}^{t}+\sum_{k=1}^{n}\frac{\alpha^{k}}{k!}\frac{\partial^{k}\langle
s_{i,t}s_{I}\rangle_{\alpha=0}^{t}}{\partial\alpha^{k}}+\mathcal{O}(\alpha^{(n+1)}).$
(SI.1.74)

Evaluating it at $\alpha=1$ knowing that $\langle
s_{i,t}s_{I}\rangle_{0}^{t}=\langle s_{i,t}s_{I}\rangle_{1}^{t},\forall i,I$,
we have

$\sum_{k=1}^{n}\left[\frac{\alpha^{k}}{k!}\right]_{\alpha=1}\frac{\partial^{k}\langle
s_{i,t}s_{I}\rangle_{0}^{t}}{\partial\alpha^{k}}=0+\left[\mathcal{O}(\alpha^{(n+1)})\right]_{\alpha=1}.$
(SI.1.75)

Let us define $\Delta h_{i,t}=-\theta_{i,t}+h_{i,t}$, we have that

$\displaystyle\frac{\partial P_{\alpha}(\mathbf{s})}{\partial\alpha}=$
$\displaystyle\sum_{t}\sum_{i\in\mathcal{I}_{t}}\left(s_{i,t}\Delta
h_{i,t}-\langle s_{i,t}\Delta
h_{i,t}\rangle_{t,\alpha}\right)P_{\alpha}(\mathbf{s}),$ (SI.1.76)
$\displaystyle\frac{\partial^{2}P_{\alpha}(\mathbf{s})}{\partial\alpha^{2}}=$
$\displaystyle\left(\sum_{t}\sum_{i\in\mathcal{I}_{t}}s_{i,t}\Delta
h_{i,t}-\langle s_{i,t}\Delta
h_{i,t}\rangle_{t,\alpha}\right)^{2}P_{\alpha}(\mathbf{s})$ (SI.1.77)
$\displaystyle-\sum_{t}\left(\left\langle\left(\sum_{i\in\mathcal{I}_{t}}s_{i,t}\Delta
h_{i,t}\right)^{2}\right\rangle_{t,\alpha}-\left\langle\sum_{i\in\mathcal{I}_{t}}s_{i,t}\Delta
h_{i,t}\right\rangle_{t,\alpha}^{2}\right)P_{\alpha}(\mathbf{s})$ (SI.1.78)
$\displaystyle=$
$\displaystyle\sum_{t}\Bigg{(}\left(\left(\sum_{i\in\mathcal{I}_{t}}s_{i,t}\Delta
h_{i,t}\right)^{2}-\left\langle\left(\sum_{i\in\mathcal{I}_{t}}s_{i,t}\Delta
h_{i,t}\right)^{2}\right\rangle_{t,\alpha}\right)$ (SI.1.79)
$\displaystyle-2\left\langle\sum_{i\in\mathcal{I}_{t}}s_{i,t}\Delta
h_{i,t}\right\rangle_{t,\alpha}\left(\sum_{k\in\mathcal{I}_{t}}s_{k,t}\Delta
h_{k,t}-\left\langle\sum_{k\in\mathcal{I}_{t}}s_{k,t}\Delta
h_{k,t}\right\rangle_{t,\alpha}\right)\Bigg{)}P_{\alpha}(\mathbf{s}),$
(SI.1.80)

where $\langle\dots\rangle_{t,\alpha}=\sum_{\mathbf{s}_{t}}\dots
P_{\alpha}(\mathbf{s}_{t}|\mathbf{s}_{1:t-1})$. From these equations we derive
the first and second order approximations.

For the first order term, we have

$\displaystyle\frac{\partial\langle
s_{i,t}s_{I}\rangle_{\alpha=0}^{t}}{\partial\alpha}=$
$\displaystyle\sum_{\mathbf{s}_{1:T}}s_{i,t}s_{I}\frac{\partial
P_{0}(\mathbf{s})}{\partial\alpha}$ (SI.1.81) $\displaystyle=$
$\displaystyle\sum_{\mathbf{s}_{1:T}}\left(\sum_{t}\sum_{k\in\mathcal{I}_{t}}\left(s_{i,t}s_{I}(s_{k,t}\Delta
h_{k,t}-\langle s_{k,t}\Delta
h_{k,t}\rangle_{t,0})\right)\right)P_{0}(\mathbf{s}).$ (SI.1.82)

The second order term is:

$\displaystyle\frac{\partial^{2}\langle
s_{i,t}s_{I}\rangle_{\alpha=0}^{t}}{\partial\alpha^{2}}=$
$\displaystyle\sum_{\mathbf{s}_{1:T}}s_{i,t}s_{I}\frac{\partial^{2}P_{0}(\mathbf{s})}{\partial\alpha^{2}}$
(SI.1.83) $\displaystyle=$
$\displaystyle\sum_{\mathbf{s}_{1:T}}s_{i,t}s_{I}\Bigg{(}\Bigg{(}\sum_{k\in\mathcal{I}_{t}}s_{k,t}\Delta
h_{k,t}\Bigg{)}^{2}-\left\langle\Bigg{(}\sum_{k\in\mathcal{I}_{t}}s_{k,t}\Delta
h_{k,t}\Bigg{)}^{2}\right\rangle_{t,0}$ (SI.1.84)
$\displaystyle-2\left\langle\sum_{k\in\mathcal{I}_{t}}s_{k,t}\Delta
h_{k,t}\right\rangle_{t,0}\Bigg{(}\sum_{m\in\mathcal{I}_{t}}s_{m,t}\Delta
h_{m,t}-\left\langle\sum_{m\in\mathcal{I}_{t}}s_{m,t}\Delta
h_{m,t}\right\rangle_{t,0}\Bigg{)}\Bigg{)}P_{0}(\mathbf{s}).$ (SI.1.85)

Using the equations above to solve Supplementary Eq. SI.1.75 for different
orders and choices of $P_{0}$ will give us the different Plekfa
approximations.

##  Supplementary Note 2 Plefka[$t-1,t$]

This approximation uses the following approximated marginal probability
distribution:

$P_{\alpha}^{[t-1:t]}(\mathbf{s}_{t},\mathbf{s}_{t-1})=\sum_{\mathbf{s}_{t-2}}P_{\alpha}(\mathbf{s}_{t}|\mathbf{s}_{t-1})P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$
(SI.2.86)

where

$\displaystyle P_{\alpha}(s_{i,t}|\mathbf{s}_{t-1})=$
$\displaystyle\frac{\mathrm{e}^{s_{i,t}{h}_{i,t}(\alpha)}}{2\cosh{h}_{i,t}(\alpha)},$
(SI.2.87) $\displaystyle h_{i,t}(\alpha)=$
$\displaystyle(1-\alpha)\Theta_{i,t}+\alpha(H_{i}+\sum_{j}J_{ij}s_{j,t-1}).$
(SI.2.88)

Here, by increasing the value of $\alpha$ from $0$ to $1$, one can smoothly
connect the independent and coupled models. Further, $h_{i,t}(\alpha)$ can be
written as

$h_{i,t}(\alpha)=\Theta_{i,t}+\alpha\,\Delta h_{i,t},$ (SI.2.89)

where $\Delta h_{i,t}=-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}s_{j,t-1}$ represents
deviation from the independent model. We estimate $m_{i,t}$ by using its
$\alpha$-dependent approximation, defined as

$m_{i,t}(\alpha)=\sum_{\mathbf{s}_{t},\mathbf{s}_{t-1}}s_{i,t}P_{\alpha}^{[t-1:t]}(\mathbf{s}_{t},\mathbf{s}_{t-1})=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\tanh
h_{i,t}(\alpha)P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}).$
(SI.2.90)

Approximating its value by expanding around $\alpha=0$ yields

$m_{i,t}(\alpha)=m_{i,t}(\alpha=0)+\sum_{k=1}^{n}\frac{\alpha^{k}}{k!}\frac{\partial^{k}m_{i,t}(\alpha=0)}{\partial\alpha^{k}}+\mathcal{O}(\alpha^{(n+1)}),$
(SI.2.91)

By noting that $m_{i,t}(\alpha=0)=m_{i,t}(\alpha=1)$, we evaluate it at
$\alpha=1$. This results in

$\left[\sum_{k=1}^{n}\frac{\alpha^{k}}{k!}\frac{\partial^{k}m_{i,t}(\alpha=0)}{\partial\alpha^{k}}\right]_{\alpha=1}=0+\left[\mathcal{O}(\alpha^{(n+1)})\right]_{\alpha=1}.$
(SI.2.92)

The approximation yields the nMF equations when we ignore quadratic and higher
order terms in solving this equation, and the TAP equations when ignoring
third and higher order terms.

The first order derivative of $m_{i,t}(\alpha)$ is given as

$\displaystyle\frac{\partial m_{i,t}(\alpha)}{\partial\alpha}$
$\displaystyle=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\left[\frac{\partial\tanh
h_{i,t}(\alpha)}{\partial\alpha}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})+\tanh
h_{i,t}(\alpha)\frac{\partial
P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha}\right]P(\mathbf{s}_{t-2}).$
(SI.2.93)

Using the following equation,

$\displaystyle\frac{\partial\tanh h_{i,t}(\alpha)}{\partial\alpha}$
$\displaystyle=(1-\tanh^{2}h_{i,t}(\alpha))\Delta h_{i,t},$ (SI.2.94)

the first order derivative is given as

$\displaystyle\frac{\partial m_{i,t}(\alpha)}{\partial\alpha}$
$\displaystyle=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\Bigg{[}(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}+\tanh h_{i,t}(\alpha)\frac{\partial
P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha}\Bigg{]}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}).$
(SI.2.95)

Expectation of the first term at $\alpha=0$ is

$\displaystyle\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}(1-\tanh^{2}h_{i,t}(0))\Delta
h_{i,t}P_{0}(\mathbf{s}_{t-1})=(1-m_{i,t}^{2})(-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}).$
(SI.2.96)

The second term becomes zero at $\alpha=0$, since the derivative of
$P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})$ becomes independent of $\tanh
h_{i,t}(\alpha)$ for $\alpha=0$ and we know that
$\sum_{\mathbf{s}_{t-1}}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})=1$. Thus
we have

$\displaystyle\frac{\partial m_{i,t}(\alpha=0)}{\partial\alpha}$
$\displaystyle=(1-m_{i,t}^{2})(-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}).$
(SI.2.97)

From here, we obtain the nMF equations, yielding
$\left[\alpha(-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1})\right]_{\alpha=1}=0+\left[\mathcal{O}(\alpha^{2})\right]_{\alpha=1}$
and

$\displaystyle\Theta_{i,t}$
$\displaystyle=H_{i}+\sum_{j}J_{ij}m_{j,t-1}+\left[\mathcal{O}(\alpha^{1})\right]_{\alpha=1},$
(SI.2.98) $\displaystyle m_{i,t}$
$\displaystyle\approx\tanh[H_{i}+\sum_{j}J_{ij}m_{j,t-1}].$ (SI.2.99)

The second order derivative of $m_{i,t}(\alpha)$ is given as

$\displaystyle\frac{\partial^{2}m_{i,t}(\alpha)}{\partial\alpha^{2}}$
$\displaystyle=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\Bigg{[}\frac{\partial^{2}\tanh
h_{i,t}(\alpha)}{\partial\alpha^{2}}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})+2\frac{\partial\tanh
h_{i,t}(\alpha)}{\partial\alpha}\frac{\partial
P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha}$ (SI.2.100)
$\displaystyle\phantom{=}+\tanh
h_{i,t}(\alpha)\frac{\partial^{2}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha^{2}}\Bigg{]}P(\mathbf{s}_{t-2}).$
(SI.2.101)

Here we note that

$\displaystyle\frac{\partial^{2}\tanh h_{i,t}(\alpha)}{\partial\alpha^{2}}$
$\displaystyle=-2\tanh h_{i,t}(\alpha)(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}^{2}.$ (SI.2.102)

From these equations, the second derivative is computed as

$\displaystyle\frac{\partial^{2}m_{i,t}(\alpha)}{\partial\alpha^{2}}=$
$\displaystyle\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\Big{(}-2\tanh
h_{i,t}(\alpha)(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}^{2}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})$ (SI.2.103)
$\displaystyle\phantom{=}+2(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}\frac{\partial
P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha}$ (SI.2.104)
$\displaystyle\phantom{=}+\frac{\partial^{2}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha^{2}}\Big{)}P(\mathbf{s}_{t-2}).$
(SI.2.105)

We evaluate the second derivative at $\alpha=0$. Here at $\alpha=0$ the third
term is zero since
$\sum_{\mathbf{s}_{t-1}}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})=1$, so
its derivatives are equal to zero. Thus we have

$\displaystyle\frac{\partial^{2}m_{i,t}(\alpha=0)}{\partial\alpha^{2}}=$
$\displaystyle-2m_{i,t}(1-m_{i,t}^{2})((-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1})^{2}+\sum_{j}J_{ij}^{2}(1-m_{j,t-1}^{2}))$
(SI.2.106)
$\displaystyle\phantom{=}+2(1-m_{i,t}^{2})\sum_{j}J_{ij}\frac{\partial
m_{j,t-1}(\alpha)}{\partial\alpha},$ (SI.2.107)

where the last term comes from the $s_{j,t-1}$ terms in $h_{i,t}(\alpha)$
multiplied by $\frac{\partial
P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha}$.

Note that the second order term in Supplementary Eq. SI.2.107 contains the
expression
$\left[\alpha^{2}(-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1})^{2}\right]_{\alpha=1}$
and $\left[\alpha^{2}\frac{\partial
m_{k,t-1}(\alpha)}{\partial\alpha}\right]_{\alpha=1}$ that can be neglected as
terms with order higher than quadratic. This is due to the fact that the
second-order approximation is in the proximity of naive mean-field solution
which is in the first order of $\alpha$. Therefore, we know that
$\left[\alpha(1-m_{i}^{2})(-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1})\right]_{\alpha=1}=\left[\mathcal{O}(\alpha^{2})\right]_{\alpha=1}$
and thus
$\big{[}\alpha^{2}(-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1})^{2}\big{]}_{\alpha=1}=\left[\mathcal{O}(\alpha^{4})\right]_{\alpha=1}$
and $\left[\alpha^{2}\frac{\partial
m_{k,t-1}(\alpha)}{\partial\alpha}\right]_{\alpha=1}=\left[\mathcal{O}(\alpha^{3})\right]_{\alpha=1}$
which can be neglected for the second order approximation.

The combination of the first and second order derivatives of $m_{i,t}$
evaluated at $\alpha=0$ allows to solve Supplementary Eq. SI.2.92 for order
$n=2$, yielding the TAP equations:

$\displaystyle\Theta_{i,t}$
$\displaystyle=H_{i}+\sum_{j}J_{ij}m_{j,t-1}-m_{i,t}\sum_{jl}J_{ij}^{2}(1-m_{j,t-1}^{2})+\left[\mathcal{O}(\alpha^{2})\right]_{\alpha=1},$
(SI.2.108) $\displaystyle m_{i,t}$
$\displaystyle\approx\tanh[H_{i}+\sum_{j}J_{ij}m_{j,t-1}-m_{i,t}\sum_{jl}J_{ij}^{2}(1-m_{j,t-1}^{2})].$
(SI.2.109)

These results have a form similar to the TAP equations obtained for symmetric
and asymmetric networks [38, 39].

### Equal-time correlations

Here we approximate $C_{ik,t}$ by evaluating the Plefka expansion of
$C_{ik,t}(\alpha=1)$ around $\alpha=0$.

When $i=k$, we have $C_{ii,t}(\alpha)=1-m_{i,t}(\alpha)^{2}$. When $i\neq k$,
correlations in the system can be obtained by expanding the alpha-dependent
correlations

$C_{ik,t}(\alpha)=\sum_{\mathbf{s}_{t},\mathbf{s}_{t-1}}(s_{i,t}-m_{i,t}(\alpha))(s_{k,t}-m_{k,t}(\alpha))P_{\alpha}^{[t-1:t]}(\mathbf{s}_{t},\mathbf{s}_{t-1}),$
(SI.2.110)

around $\alpha=0$:

$C_{ik,t}(\alpha)=\sum_{n=0}^{m}\frac{\alpha^{n}}{n!}\frac{\partial^{n}C_{ik,t}(\alpha=0)}{\partial\alpha^{n}}+\mathcal{O}(\alpha^{(m+1)}).$
(SI.2.111)

The derivatives with different orders, and their values evaluated at
$\alpha=0$ are obtained as follows.

The zeroth order term is:

$\displaystyle C_{ik,t}(\alpha)=$
$\displaystyle\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}(\tanh
h_{i,t}(\alpha)-m_{i,t}(\alpha))(\tanh
h_{k,t}(\alpha)-m_{k,t}(\alpha))P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$
(SI.2.112) $\displaystyle C_{ik,t}(\alpha=0)=$
$\displaystyle(m_{i,t}-m_{i,t})(m_{k,t}-m_{k,t})=0.$ (SI.2.113)

The first order term is:

$\displaystyle\frac{\partial C_{ik,t}(\alpha)}{\partial\alpha}$
$\displaystyle=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\Big{(}\big{(}(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}-\frac{\partial m_{i,t}(\alpha)}{\partial\alpha}\big{)}\big{(}\tanh
h_{k,t}(\alpha)-m_{k,t}(\alpha)\big{)}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})$
(SI.2.114) $\displaystyle\phantom{=}+\big{(}(1-\tanh^{2}h_{k,t}(\alpha))\Delta
h_{k,t}-\frac{\partial m_{k,t}(\alpha)}{\partial\alpha}\big{)}\big{(}\tanh
h_{i,t}(\alpha)-m_{i,t}(\alpha)\big{)}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})$
(SI.2.115) $\displaystyle\phantom{=}+\sum_{m}(\tanh
h_{i,t}(\alpha)-m_{i,t}(\alpha))(\tanh
h_{k,t}(\alpha)-m_{k,t}(\alpha))\frac{\partial
P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha}\Big{)}P(\mathbf{s}_{t-2}).$
(SI.2.116)

Here, all terms are equal to zero for $\alpha=0$ because the terms $\tanh
h_{i,t}(\alpha=0)=m_{i,t}$ cancel out. Then we have

$\displaystyle\frac{\partial C_{ik,t}(\alpha=0)}{\partial\alpha}$
$\displaystyle=0.$ (SI.2.117)

Therefore the nMF equation is obtained as

$\displaystyle C_{ik,t}$ $\displaystyle\approx 0.$ (SI.2.118)

The second order term is:

$\displaystyle\frac{\partial^{2}C_{ik,t}(\alpha)}{\partial\alpha^{2}}=$
$\displaystyle\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\Big{(}2\big{(}(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}-\frac{\partial
m_{i,t}(\alpha)}{\partial\alpha}\big{)}\big{(}(1-\tanh^{2}h_{k,t}(\alpha))\Delta
h_{k,t}-\frac{\partial
m_{k,t}(\alpha)}{\partial\alpha}\big{)}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})$
(SI.2.119) $\displaystyle+\big{(}-2\tanh
h_{i,t}(\alpha)(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}^{2}-\frac{\partial^{2}m_{i,t}(\alpha)}{\partial\alpha^{2}}\big{)}\big{(}\tanh
h_{k,t}(\alpha)-m_{k,t}(\alpha)\big{)}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})$
(SI.2.120) $\displaystyle+\big{(}-2\tanh
h_{k,t}(\alpha)(1-\tanh^{2}h_{k,t}(\alpha))\Delta
h_{k,t}^{2}-\frac{\partial^{2}m_{k,t}(\alpha)}{\partial\alpha^{2}}\big{)}\big{(}\tanh
h_{i,t}(\alpha)-m_{i,t}(\alpha)\big{)}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})$
(SI.2.121) $\displaystyle+(\tanh h_{i,t}(\alpha)-m_{i,t}(\alpha))(\tanh
h_{k,t}(\alpha)-m_{k,t}(\alpha))\frac{\partial^{2}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha^{2}}$
(SI.2.122) $\displaystyle+\big{(}(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}-\frac{\partial m_{i,t}(\alpha)}{\partial\alpha}\big{)}(\tanh
h_{k,t}(\alpha)-m_{k,t}(\alpha))\frac{\partial
P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha}$ (SI.2.123)
$\displaystyle+(\tanh
h_{i,t}(\alpha)-m_{i,t}(\alpha))\big{(}(1-\tanh^{2}h_{k,t}(\alpha))\Delta
h_{k,t}-\frac{\partial m_{k,t}(\alpha)}{\partial\alpha}\big{)}\frac{\partial
P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha}\Big{)}P(\mathbf{s}_{t-2}).$
(SI.2.124)

Except for the term in the first line, for $\alpha=0$ all terms are equal to
zero because $\tanh h_{i,t}(\alpha=0)=m_{i,t}$ cancel out. This gives

$\displaystyle\frac{\partial^{2}C_{ik,t}(\alpha=0)}{\partial\alpha^{2}}=$
$\displaystyle
2(1-m_{i,t}^{2})(1-m_{k,t}^{2})\sum\limits_{j}J_{ij}J_{kj}(1-m_{j,t-1}^{2}).$
(SI.2.125)

Hence the correlations expanded up to the second order can be described as:

$\displaystyle C_{ik,t}(\alpha)$
$\displaystyle=C_{ik,t}(\alpha=0)+\alpha\frac{\partial
C_{ik,t}(\alpha=0)}{\partial\alpha}+\frac{\alpha^{2}}{2}\frac{\partial^{2}C_{ik,t}(\alpha=0)}{\partial\alpha^{2}}+\mathcal{O}(\alpha^{3})$
(SI.2.126)
$\displaystyle=\alpha^{2}(1-m_{i,t}^{2})(1-m_{k,t}^{2})\sum\limits_{j}J_{ij}J_{kj}(1-m_{j,t-1}^{2})+\mathcal{O}(\alpha^{3}).$
(SI.2.127)

Hence, the corresponding TAP approximation is obtained as

$\displaystyle C_{ik,t}$
$\displaystyle\approx(1-m_{i,t}^{2})(1-m_{k,t}^{2})\sum\limits_{j}J_{ij}J_{kj}(1-m_{j,t-1}^{2}).$
(SI.2.128)

### Time-delayed correlations

Similarly to the equal-time correlations, we approximate $D_{il,t}$ by
evaluating the Plefka expansion of $D_{il,t}(\alpha=1)$ around $\alpha=0$.

We describe time-delayed correlations of the system

$D_{il,t}(\alpha)=\sum_{\mathbf{s}_{t},\mathbf{s}_{t-1}}(s_{i,t}-m_{i,t}(\alpha))(s_{l,t-1}-m_{l,t-1}(\alpha))P_{\alpha}^{[t-1:t]}(\mathbf{s}_{t},\mathbf{s}_{t-1}),$
(SI.2.129)

using an expansion:

$D_{il,t}(\alpha)=\sum_{k=0}^{n}\frac{\alpha^{k}}{k!}\frac{\partial^{k}D_{il,t}(\alpha=0)}{\partial\alpha^{k}}+\mathcal{O}(\alpha^{(n+1)}).$
(SI.2.130)

Likewise, the zeroth order term yields:

$\displaystyle D_{il,t}(\alpha)$
$\displaystyle=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}(\tanh
h_{i,t}(\alpha)-m_{i,t}(\alpha))(s_{l,t-1}-m_{l,t-1}(\alpha))P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}).$
(SI.2.131) $\displaystyle D_{il,t}(\alpha=0)$ $\displaystyle=0.$ (SI.2.132)

The first order term is:

$\displaystyle\frac{\partial D_{il,t}(\alpha)}{\partial\alpha}$
$\displaystyle=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\Big{(}\big{(}(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}-\frac{\partial
m_{i,t}(\alpha)}{\partial\alpha}\big{)}(s_{l,t-1}-m_{l,t-1}(\alpha))P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})$
(SI.2.133) $\displaystyle\phantom{=}+(\tanh
h_{i,t}(\alpha)-m_{i,t}(\alpha))\big{(}(s_{l,t-1}-m_{l,t-1}(\alpha))\frac{\partial
P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha}-\frac{\partial
m_{l,t-1}(\alpha)}{\partial\alpha}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})\big{)}\Big{)}P(\mathbf{s}_{t-2}).$
(SI.2.134)

When evaluated at $\alpha=0$, the terms in the second line disappear because
the terms $\tanh h_{i,t}(\alpha=0)=m_{i,t}$ cancel out. Note that the term in
the first line is computed as

$\big{(}(1-\tanh^{2}h_{i,t}(0))\Delta h_{i,t}-\frac{\partial
m_{i,t}(\alpha)}{\partial\alpha}=(1-m_{i,t}^{2})\sum_{j}J_{ij}(s_{j,t-1}-m_{j,t-1}),$
(SI.2.135)

which is multiplied by $s_{l,t-1}-m_{l,t-1}(\alpha)$. Because we take
expectation over $\mathbf{s}_{t-1}$ with the independent model, only
expectation of the term $(s_{l,t-1}-m_{l,t-1})^{2}$ is preserved. Thus we have

$\displaystyle\frac{\partial D_{il,t}(\alpha=0)}{\partial\alpha}$
$\displaystyle=(1-m_{i,t}^{2})J_{il}(1-m_{l,t-1}^{2}).$ (SI.2.136)

Therefore the nMF equation is obtained as:

$D_{il,t}\approx J_{il}(1-m_{i,t}^{2})(1-m_{l,t-1}^{2}).$ (SI.2.137)

The second order term is:

$\displaystyle\frac{\partial^{2}D_{il,t}(\alpha)}{\partial\alpha^{2}}=$
$\displaystyle\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\Big{(}-2\tanh
h_{i,t}(\alpha)(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}^{2}-\frac{\partial^{2}m_{i,t}(\alpha)}{\partial\alpha^{2}}\big{)}(s_{l,t-1}-m_{l,t-1}(\alpha))P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})$
(SI.2.138) $\displaystyle+2\big{(}(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}-\frac{\partial
m_{i,t}(\alpha)}{\partial\alpha}\big{)}\big{(}(s_{l,t-1}-m_{l,t-1}(\alpha))\frac{\partial
P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha}$ (SI.2.139)
$\displaystyle-\frac{\partial
m_{l,t-1}(\alpha)}{\partial\alpha}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})\big{)}$
(SI.2.140) $\displaystyle+(\tanh
h_{i,t}(\alpha)-m_{i,t}(\alpha))\Big{(}\big{(}(s_{l,t-1}-m_{l,t-1}(\alpha))\frac{\partial^{2}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha^{2}}$
(SI.2.141) $\displaystyle-2\frac{\partial
m_{l,t-1}(\alpha)}{\partial\alpha}\frac{\partial
P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha}-\frac{\partial^{2}m_{l,t-1}(\alpha)}{\partial\alpha^{2}}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})\Big{)}P(\mathbf{s}_{t-2}).$
(SI.2.142)

Here, in the first line, only the $\Delta h_{i,t}^{2}(s_{l,t-1}-m_{l,t-1})$
survives as

$\displaystyle\Delta h_{i,t}^{2}(s_{l,t-1}-m_{l,t-1})$
$\displaystyle=\Big{(}-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}+\sum_{j}J_{ij}(s_{j,t-1}-m_{j,t-1})\Big{)}^{2}(s_{l,t-1}-m_{l,t-1})$
(SI.2.143)
$\displaystyle=\Big{(}-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}\Big{)}^{2}(s_{l,t-1}-m_{l,t-1})$
(SI.2.144)
$\displaystyle+2\Big{(}-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}\Big{)}\sum_{j}J_{ij}(s_{j,t-1}-m_{j,t-1})(s_{l,t-1}-m_{l,t-1})$
(SI.2.145)
$\displaystyle+\Big{(}\sum_{j}J_{ij}(s_{j,t-1}-m_{j,t-1})\Big{)}^{2}(s_{l,t-1}-m_{l,t-1}),$
(SI.2.146)

which results in
$2(-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1})(1-m_{l,t-1}^{2})-2J_{il}^{2}m_{l,t-1}(1-m_{l,t-1}^{2})$
when evaluated for the mean field model.

The term in the second and third lines of the previous equation can be
decomposed in two terms. The second part results in
$\sum_{j}J_{ij}(s_{j,t-1}-m_{j,t-1})\frac{\partial
m_{l,t-1}(\alpha)}{\partial\alpha}$, which results in zero. From the first
term, only survives the part containing
$\sum_{j}J_{ij}(s_{j,t-1}-m_{j,t-1})(s_{l,t-1}-m_{l,t-1})\frac{\partial
P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})}{\partial\alpha}$ which for
$\alpha=0$ results in

$\displaystyle\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\sum_{j}J_{ij}(s_{j,t-1}-m_{j,t-1})(s_{l,t-1}-m_{l,t-1})\dfrac{\partial{P_{0}(s_{t-1}|s_{t-2})}}{\partial{\alpha}}P(s_{t-2})$
(SI.2.147)
$\displaystyle=\sum_{\mathbf{s}_{t-2}}\sum_{jk}J_{ij}(s_{j,t-1}-m_{j,t-1})(s_{l,t-1}-m_{l,t-1})(s_{k,t-1}-m_{k,t-1})\Delta
h_{k,t-1}Q(s_{t-1})P(s_{t-2})$ (SI.2.148)
$\displaystyle=\sum_{\mathbf{s}_{t-1}}J_{il}(s_{l,t-1}-m_{l,t-1})^{3}Q(s_{t-1})\sum_{\mathbf{s}_{t-2}}\Delta
h_{k,t-1}P(s_{t-2})$ (SI.2.149)
$\displaystyle=-2J_{il}m_{l,t-1}(1-m_{l,t-1}^{2})(-\Theta_{i,t-1}+H_{i}+\sum_{j}J_{ij}m_{j,t-2}).$
(SI.2.150)

The term in the last two lines disappear because the terms $\tanh
h_{i,t}(\alpha=0)=m_{i,t}$ cancel out. Thus, evaluated at $\alpha=0$ we have

$\displaystyle\frac{\partial^{2}D_{il,t}(\alpha=0)}{\partial\alpha^{2}}=$
$\displaystyle 4m_{i,t}(1-m_{i,t}^{2})J_{il}^{2}m_{l,t-1}(1-m_{l,t-1}^{2})$
(SI.2.151)
$\displaystyle-4m_{i,t}(1-m_{i,t}^{2})J_{il}(1-m_{l,t-1}^{2})(-\Theta_{i,t}+H_{i}+\sum_{j}m_{j,t-1})$
(SI.2.152)
$\displaystyle-4(1-m_{i,t}^{2})J_{il}m_{l,t-1}(1-m_{l,t-1}^{2})(-\Theta_{i,t-1}+H_{i}+\sum_{j}J_{ij}m_{j,t-2}).$
(SI.2.153)

The second and third term above are equal to $0+\mathcal{O}(\alpha)$, which
makes them negligible (i.e. order larger than quadratic) when computing
$\alpha^{2}\frac{\partial^{2}D_{il,t}(\alpha=0)}{\partial\alpha^{2}}$.

Hence the second order expansion is

$\displaystyle D_{il,t}(\alpha)$
$\displaystyle=D_{il,t}(\alpha=0)+\alpha\frac{\partial
D_{il,t}(\alpha=0)}{\partial\alpha}+\frac{\alpha^{2}}{2}\frac{\partial^{2}D_{il,t}(\alpha=0)}{\partial\alpha^{2}}+\mathcal{O}(\alpha^{3})$
(SI.2.154)
$\displaystyle=\alpha(1-m_{i,t}^{2})J_{il}(1-m_{l,t-1}^{2})+\alpha^{2}2m_{i,t}(1-m_{i,t}^{2})J_{il}^{2}m_{l,t-1}(1-m_{l,t-1}^{2})+\mathcal{O}(\alpha^{3}).$
(SI.2.155)

Thus the TAP equation for the time-delayed correlations can be described as:

$D_{il,t}\approx
J_{il}(1-m_{i,t}^{2})(1-m_{l,t-1}^{2})(1+2J_{il}m_{i,t}m_{l,t-1}).$ (SI.2.156)

##  Supplementary Note 3 Plefka[$t$]

This approximation uses the following approximated marginal probability
distribution:

$P_{\alpha}^{[t]}(\mathbf{s}_{t})=\sum_{\mathbf{s}_{t-1}}P_{\alpha}(\mathbf{s}_{t}|\mathbf{s}_{t-1})P(\mathbf{s}_{t-1}),$
(SI.3.157)

where

$\displaystyle P_{\alpha}(s_{i,t}|\mathbf{s}_{t-1})=$
$\displaystyle\frac{\mathrm{e}^{s_{i,t}{h}_{i,t}(\alpha)}}{2\cosh{h}_{i,t}(\alpha)},$
(SI.3.158) $\displaystyle h_{i,t}(\alpha)=$
$\displaystyle(1-\alpha)\Theta_{i,t}+\alpha(H_{i}+\sum_{j}J_{ij}s_{j,t-1}).$
(SI.3.159)

We define $h_{i,t}(\alpha)=\Theta_{i,t}+\alpha\,\Delta h_{i,t}$, where $\Delta
h_{i,t}=-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}s_{j,t-1}$ represents deviation from
the independent model. We estimate $m_{i,t}$ by using its $\alpha$-dependent
approximation, whose element is defined as

$m_{i,t}(\alpha)=\sum_{\mathbf{s}_{t},\mathbf{s}_{t-1}}s_{i,t}P_{\alpha}^{[t]}(\mathbf{s}_{t},\mathbf{s}_{t-1})=\sum_{\mathbf{s}_{t-1}}\tanh
h_{i,t}(\alpha)P(\mathbf{s}_{t-1}).$ (SI.3.160)

Approximating its value by expanding around $\alpha=0$ yields

$m_{i,t}(\alpha)=m_{i,t}(\alpha=0)+\sum_{k=1}^{n}\frac{\alpha^{k}}{k!}\frac{\partial^{k}m_{i,t}(\alpha=0)}{\partial\alpha^{k}}+\mathcal{O}(\alpha^{(n+1)}).$
(SI.3.161)

By noting that $m_{i,t}(\alpha=0)=m_{i,t}(\alpha=1)$, we evaluate it at
$\alpha=1$. This results in

$\left[\sum_{k=1}^{n}\frac{\alpha^{k}}{k!}\frac{\partial^{k}m_{i,t}(\alpha=0)}{\partial\alpha^{k}}\right]_{\alpha=1}=0+\left[\mathcal{O}(\alpha^{(n+1)})\right]_{\alpha=1}.$
(SI.3.162)

The approximation yields the nMF equations when we ignore the quadratic term
and higher, and the TAP equations when ignoring the third and higher order
terms. The terms at each order and its evaluation at $\alpha=0$ are obtained
as follows.

The first order term is:

$\displaystyle\frac{\partial m_{i,t}(\alpha)}{\partial\alpha}$
$\displaystyle=\sum_{\mathbf{s}_{t-1}}\Big{(}1-\tanh^{2}{h_{i,t}(\alpha)}\Big{)}\Delta
h_{i,t}P(\mathbf{s}_{t-1}),$ (SI.3.163) $\displaystyle\frac{\partial
m_{i,t}(\alpha=0)}{\partial\alpha}$
$\displaystyle=(1-m_{i,t}^{2})\Big{(}-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}\Big{)}.$
(SI.3.164)

The first order or naive mean-field approximation is then

$\displaystyle\Theta_{i,t}$
$\displaystyle=H_{i}+\sum_{j}J_{ij}m_{j,t-1}+\left[\mathcal{O}(\alpha^{1})\right]_{\alpha=1},$
(SI.3.165) $\displaystyle m_{i}(t)$
$\displaystyle\approx\tanh{[H_{i}+\sum_{j}J_{ij}m_{j,t-1}]}.$ (SI.3.166)

For the second order approximation, we have

$\displaystyle\frac{\partial^{2}m_{i,t}(\alpha)}{\partial\alpha^{2}_{t}}$
$\displaystyle=-2\sum_{\mathbf{s}_{t-1}}\tanh{h_{i,t}(\alpha)}\Big{(}1-\tanh^{2}{h_{i,t}(\alpha)}\Big{)}\Delta
h_{i,t}^{2}P(\mathbf{s}_{t-1}),$ (SI.3.167)
$\displaystyle\frac{\partial^{2}m_{i,t}(\alpha=0)}{\partial\alpha^{2}_{t}}$
$\displaystyle=-2m_{i,t}\Big{(}1-m_{i,t}^{2}\Big{)}\bigg{[}\Big{(}-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}\Big{)}^{2}$
(SI.3.168)
$\displaystyle\phantom{=}+\sum_{j,k}J_{ij}J_{ik}C_{jk,t-1}\bigg{]}.$
(SI.3.169)

Then we solve the following second order equation:

$\frac{\partial
m_{i,t}(\alpha=0)}{\partial\alpha}+\frac{1}{2}\frac{\partial^{2}m_{i,t}(\alpha=0)}{\partial\alpha^{2}}=0+\mathcal{O}(\alpha^{3})\rvert_{\alpha=1}.$
(SI.3.170)

To compute the second order approximation, we can take advantage of the fact
that the second order term contains the expression
$\left[\alpha^{2}(-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j})^{2}\right]_{\alpha=1}$.
We know that
$\left[\alpha(1-m_{i}^{2})(-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j})\right]_{\alpha=1}=\left[\mathcal{O}(\alpha^{2})\right]_{\alpha=1}$
and thus
$\left[\alpha^{2}(-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j})^{2}\right]_{\alpha=1}=\left[\mathcal{O}(\alpha^{4})\right]_{\alpha=1}$.
This yields the TAP equations:

$\displaystyle\Theta_{i,t}$
$\displaystyle=H_{i}+\sum_{j}J_{ij}m_{j,t-1}-m_{i,t}\sum_{j,k}J_{ij}J_{ik}C_{jk,t-1}+\left[\mathcal{O}(\alpha^{2})\right]_{\alpha=1},$
(SI.3.171) $\displaystyle m_{i}(t)$
$\displaystyle\approx\tanh{[H_{i}+\sum_{j}J_{ij}m_{j,t-1}-m_{i,t}\sum_{j,k}J_{ij}J_{ik}C_{jk,t-1}]}.$
(SI.3.172)

This result is a reminiscence of the TAP equations obtained for symmetric
networks and asymmetric networks, which generally take the form
$\Theta_{i,t}\approx
H_{i}+\sum_{j}J_{ij}m_{j,t-1}-m_{i,t}\sum_{j}J_{ij}^{2}(1-m_{j}^{2})$ [38,
39]. The only difference is that previous results approximated either the
stationary state of the network, or the probability of trajectories over
several updates of the network dynamics. The consequence of these is that
previous results ignored correlations at previous states (since they were also
expanded from the independent model, thus $C_{jl}$ terms become zero when
$i\neq j$ and $1-m_{i}^{2}$ otherwise).

### Equal-time correlations

When $i=j$, $C_{ii,t}=1-m_{i,t}(\alpha)^{2}$. When $i\neq j$, correlations in
the system can be obtained by expanding

$C_{ik,t}(\alpha)=\sum_{\mathbf{s}_{t},\mathbf{s}_{t-1}}(s_{i,t}-m_{i,t}(\alpha))(s_{k,t}-m_{k,t}(\alpha))P_{\alpha}(\mathbf{s}_{t}|\mathbf{s}_{t-1})P(\mathbf{s}_{t-1}),$
(SI.3.173)

over $\alpha=0$:

$C_{ik,t}(\alpha)=\sum_{n=0}^{m}\frac{\alpha^{n}}{n!}\frac{\partial^{n}C_{ik,t}(\alpha=0)}{\partial\alpha^{n}}+\mathcal{O}(\alpha^{(m+1)}).$
(SI.3.174)

The expanded terms of each order and its evaluation at $\alpha=0$ are obtained
as follows.

The zeroth order term is:

$\displaystyle C_{ik,t}(\alpha)=$ $\displaystyle\sum_{\mathbf{s}_{t-1}}(\tanh
h_{i,t}(\alpha)-m_{i,t}(\alpha))(\tanh
h_{k,t}(\alpha)-m_{k,t}(\alpha))P(\mathbf{s}_{t-1}),$ (SI.3.175)
$\displaystyle C_{ik,t}(\alpha=0)=$
$\displaystyle(m_{i,t}-m_{i,t})(m_{k,t}-m_{k,t})=0.$ (SI.3.176)

The first order term yields

$\displaystyle\frac{\partial C_{ik,t}(\alpha)}{\partial\alpha}$
$\displaystyle=\sum_{\mathbf{s}_{t-1}}\Big{(}\big{(}(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}$ (SI.3.177) $\displaystyle\phantom{=}-\frac{\partial
m_{i,t}(\alpha)}{\partial\alpha}\big{)}\big{(}\tanh
h_{i,t}(\alpha)-m_{k,t}(\alpha)\big{)}$ (SI.3.178)
$\displaystyle\phantom{=}+\big{(}(1-\tanh^{2}h_{k,t}(\alpha))\Delta h_{k,t}$
(SI.3.179) $\displaystyle\phantom{=}-\frac{\partial
m_{k,t}(\alpha)}{\partial\alpha}\big{)}\big{(}\tanh
h_{i,t}(\alpha)-m_{i,t}(\alpha)\big{)}\Big{)}P(\mathbf{s}_{t-1}),$ (SI.3.180)
$\displaystyle\frac{\partial C_{ik,t}(\alpha=0)}{\partial\alpha}$
$\displaystyle=0.$ (SI.3.181)

Therefore the nMF equation is obtained as

$C_{ik,t}\approx 0.$ (SI.3.182)

The second order term is:

$\displaystyle\frac{\partial^{2}C_{ik,t}(\alpha)}{\partial\alpha^{2}}$
$\displaystyle=\sum_{\mathbf{s}_{t-1}}\Big{(}2\big{(}(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}$ (SI.3.183) $\displaystyle\phantom{=}-\frac{\partial
m_{i,t}(\alpha)}{\partial\alpha}\big{)}\big{(}(1-\tanh^{2}h_{k,t}(\alpha))\Delta
h_{k,t}-\frac{\partial m_{k,t}(\alpha)}{\partial\alpha}\big{)}$ (SI.3.184)
$\displaystyle\phantom{=}+\big{(}-2\tanh
h_{i,t}(\alpha)(1-\tanh^{2}h_{i,t}(\alpha))\Delta h_{i,t}^{2}$ (SI.3.185)
$\displaystyle\phantom{=}-\frac{\partial^{2}m_{i,t}(\alpha)}{\partial\alpha^{2}}\big{)}\big{(}\tanh
h_{k,t}(\alpha)-m_{k,t}(\alpha)\big{)}$ (SI.3.186)
$\displaystyle\phantom{=}+\big{(}-2\tanh
h_{k,t}(\alpha)(1-\tanh^{2}h_{k,t}(\alpha))\Delta h_{k,t}^{2}$ (SI.3.187)
$\displaystyle\phantom{=}-\frac{\partial^{2}m_{k,t}(\alpha)}{\partial\alpha^{2}}\big{)}\big{(}\tanh
h_{i,t}(\alpha)-m_{i,t}(\alpha)\big{)}\Big{)}P(\mathbf{s}_{t-1}),$ (SI.3.188)
$\displaystyle\frac{\partial^{2}C_{ik,t}(\alpha=0)}{\partial\alpha^{2}}$
$\displaystyle=2(1-m_{i,t}^{2})(1-m_{k,t}^{2})(\sum\limits_{ln}J_{il}J_{jn}C_{ln,t-1}).$
(SI.3.189)

So the TAP equation for the correlations can be described as:

$\displaystyle C_{ik,t}$
$\displaystyle\approx(1-m_{i,t}^{2})(1-m_{k,t-1}^{2})\sum\limits_{ln}J_{il}J_{jn}C_{ln,t-1}.$
(SI.3.190)

To obtain this result, correlations of order $3$ (and more)
$C_{ijk,t}=\sum_{\mathbf{s}_{t}}(s_{i,t}-m_{i,t})(s_{j,t}-m_{j,t})(s_{k,t}-m_{k,t})P(\mathbf{s}_{t})=0+\left[\mathcal{O}(\alpha^{3})\right]_{\alpha=1}$
are ignored in the TAP equations if $i\neq j\neq k$. Otherwise we used
$C_{iik,t}=-2m_{i,t}C_{ik,t}$ and $C_{iii,t}=-2m_{i,t}(1-m_{i,t}^{2})$.

The obtained expression has a form that is similar to the equations in [40]
but presents some differences since is computed from the same expansion as the
TAP equations, instead of performing a new expansion from the approximation of
$\Theta_{i,t}$.

### Time-delayed correlations

We introduce the $\alpha$-dependent time-delayed correlations as

$D_{il,t}(\alpha)=\sum_{\mathbf{s}_{t},\mathbf{s}_{t-1}}(s_{i,t}-m_{i,t}(\alpha))(s_{l,t-1}-m_{l,t-1})P_{\alpha}(\mathbf{s}_{t}|\mathbf{s}_{t-1})P(\mathbf{s}_{t-1}).$
(SI.3.191)

We approximate $D_{il,t}$ by expanding this equation around $\alpha=0$.

The zeroth order term yields:

$\displaystyle D_{il,t}(\alpha)=$ $\displaystyle\sum_{\mathbf{s}_{t-1}}(\tanh
h_{i,t}(\alpha)-m_{i,t}(\alpha))(s_{l,t-1}-m_{l,t-1})P(\mathbf{s}_{t-1}),$
(SI.3.192) $\displaystyle D_{il,t}(\alpha=0)=$ $\displaystyle 0.$ (SI.3.193)

The first order term is:

$\displaystyle\frac{\partial D_{il,t}(\alpha)}{\partial\alpha}$
$\displaystyle=\sum_{\mathbf{s}_{t-1}}\sum_{\mathbf{s}_{t-1}}\Big{(}\big{(}(1-\tanh^{2}h_{i,t}(\alpha))\Delta
h_{i,t}$ (SI.3.194) $\displaystyle\phantom{=}-\frac{\partial
m_{i,t}(\alpha)}{\partial\alpha}\big{)}(s_{l,t-1}-m_{l,t-1})\big{)}P(\mathbf{s}_{t-1}),$
(SI.3.195) $\displaystyle\frac{\partial D_{il,t}(\alpha=0)}{\partial\alpha}$
$\displaystyle=(1-m_{i,t}^{2})\sum\limits_{j}J_{ij}C_{jl,t-1},$ (SI.3.196)

therefore obtaining that for the nMF equation:

$D_{il,t}\approx(1-m_{i,t}^{2})\sum\limits_{j}J_{ij}C_{jl,t-1}.$ (SI.3.197)

The second order term is:

$\displaystyle\frac{\partial^{2}D_{il,t}(\alpha)}{\partial\alpha^{2}}$
$\displaystyle=\sum_{\mathbf{s}_{t-1}}\Big{(}\big{(}-2\tanh
h_{i,t}(\alpha)(1-\tanh^{2}h_{i,t}(\alpha))\Delta h_{i,t}^{2}$ (SI.3.198)
$\displaystyle\phantom{=}-\frac{\partial^{2}m_{i,t}(\alpha)}{\partial\alpha^{2}}\big{)}(s_{l,t-1}-m_{l,t-1})\Big{)}P(\mathbf{s}_{t-1})$
(SI.3.199) $\displaystyle=-2\tanh
h_{i,t}(\alpha)(1-\tanh^{2}h_{i,t}(\alpha))\sum_{\mathbf{s}_{t-1}}\Delta
h_{i,t}^{2}(s_{l,t-1}-m_{l,t-1})P(\mathbf{s}_{t-1}).$ (SI.3.200)

Here we note that

$\displaystyle\sum_{\mathbf{s}_{t-1}}\Delta
h_{i,t}^{2}(s_{l,t-1}-m_{l,t-1})P(\mathbf{s}_{t-1})$ (SI.3.201)
$\displaystyle=\sum_{\mathbf{s}_{t-1}}\Big{(}-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}+\sum_{j}J_{ij}(s_{j,t-1}-m_{j,t-1})\Big{)}^{2}(s_{l,t-1}-m_{l,t-1})P(\mathbf{s}_{t-1})$
(SI.3.202)
$\displaystyle=\Big{(}-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}\Big{)}^{2}\sum_{\mathbf{s}_{t-1}}(s_{l,t-1}-m_{l,t-1})P(\mathbf{s}_{t-1})$
(SI.3.203)
$\displaystyle+2\Big{(}-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}\Big{)}\sum_{\mathbf{s}_{t-1}}\sum_{j}J_{ij}(s_{j,t-1}-m_{j,t-1})(s_{l,t-1}-m_{l,t-1})P(\mathbf{s}_{t-1})$
(SI.3.204)
$\displaystyle+\sum_{\mathbf{s}_{t-1}}\Big{(}\sum_{j}J_{ij}(s_{j,t-1}-m_{j,t-1})\Big{)}^{2}(s_{l,t-1}-m_{l,t-1})P(\mathbf{s}_{t-1})$
(SI.3.205)
$\displaystyle=2\Big{(}-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}\Big{)}\sum_{j}J_{ij}C_{jl,t-1}+\sum_{j,k}J_{ij}J_{ik}C_{jkl,t-1}.$
(SI.3.206)

where
$C_{jkl,t-1}\equiv\sum_{\mathbf{s}_{t-1}}(s_{j,t-1}-m_{j,t-1})(s_{k,t-1}-m_{k,t-1})(s_{l,t-1}-m_{l,t-1})P(\mathbf{s}_{t-1})$.
We then obtain

$\displaystyle\frac{\partial^{2}D_{il,t}(\alpha=0)}{\partial\alpha^{2}}$
$\displaystyle=-2m_{i,t}(1-m_{i,t}^{2})\sum\limits_{j,k}J_{ij}J_{ik}C_{jkl,t-1}-4m_{i,t}(1-m_{i,t}^{2})(-\Theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j})\sum\limits_{j}J_{ij}C_{jl,t-1}.$
(SI.3.207)

Here, the second term above is equal to $0+\mathcal{O}(\alpha)$, which makes
it negligible (i.e. order larger than quadratic) when computing
$\alpha^{2}\frac{\partial^{2}D_{il,t}(\alpha=0)}{\partial\alpha^{2}}$. Thus
the second order expansion for delayed correlations can be described as:

$D_{il,t}=(1-m_{i,t}^{2})\big{(}\sum\limits_{j}J_{ij}C_{jl,t-1}-m_{i,t}\sum\limits_{j,k}J_{ij}J_{ik}C_{jkl,t-1}\big{)}+\left[\mathcal{O}(\alpha^{3})\right]_{\alpha=1}.$
(SI.3.208)

Moreover, if we consider that $C_{jkl,t-1}$ terms are going to be equal to
$0+\left[\mathcal{O}(\alpha^{3})\right]_{\alpha=1}$ when $j\neq l\neq n$, we
have the TAP approximation

$D_{il,t}\approx(1-m_{i,t}^{2})\sum\limits_{j}J_{ij}C_{jl,t-1}(1+2J_{il}m_{i,t}m_{l,t-1})+\left[\mathcal{O}(\alpha^{3})\right]_{\alpha=1}.$
(SI.3.209)

The nMF Equation for time-delayed correlations is similar to the first order
approximation obtained by [40]. The second order approximation differs since
they do a new expansion over one time step of the obtained TAP expression
obtained for the whole trajectory. In our case, since we apply the TAP
expansion for one-step updates in all cases, we can derive an expression from
the same expansion that obtains the TAP equation for updating the mean fields
of the system.

##  Supplementary Note 4 Plefka[$t-1$]

This approximation uses the following approximated marginal probability
distribution.

$P_{\alpha}^{[t-1]}(\mathbf{s}_{t},\mathbf{s}_{t-1})=\sum_{\mathbf{s}_{t-2}}P(\mathbf{s}_{t}|\mathbf{s}_{t-1})P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$
(SI.4.210)

where activity at time $t$ is normally defined by

$\displaystyle P(s_{i,t}|\mathbf{s}_{t-1})=$
$\displaystyle\frac{\mathrm{e}^{s_{i,t}{h}_{i,t}}}{2\cosh{h}_{i,t}},$
(SI.4.211) $\displaystyle h_{i,t}=$ $\displaystyle
H_{i}+\sum_{j}J_{ij}s_{j,t-1},$ (SI.4.212)

whereas activity at time $t-1$ is mediated by the parameter $\alpha$

$\displaystyle P_{\alpha}(s_{i,t-1}|\mathbf{s}_{t-2})=$
$\displaystyle\frac{\mathrm{e}^{s_{i,t-1}{h}_{i,t-1}(\alpha)}}{2\cosh{h}_{i,t-1}(\alpha)},$
(SI.4.213) $\displaystyle h_{i,t-1}(\alpha)=$
$\displaystyle(1-\alpha)\Theta_{i,t-1}+\alpha(H_{i}+\sum_{j}J_{ij}s_{j,t-2}).$
(SI.4.214)

Computing the values of $\mathbf{\Theta}_{t-1}$ is equivalent to the
calculations in 3. However, in the calculations below we will see that we will
need only $m_{i,t-1}$ to compute the statistics at time $t$. Namely,
$m_{i,t-1}$ is the only value we need to know to make the approximation at
$t-1$.

Now, we calculate $m_{i,t}$ by using its $\alpha$-dependent approximation,
defined as

$m_{i,t}(\alpha)=\sum_{\mathbf{s}_{t},\mathbf{s}_{t-1}}s_{i,t}P_{\alpha}^{[t-1]}(\mathbf{s}_{t},\mathbf{s}_{t-1})=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\tanh
h_{i,t}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}).$
(SI.4.215)

Approximating its value by expanding around $\alpha=0$ yields

$m_{i,t}(\alpha)=m_{i,t}(\alpha=0)+\sum_{k=1}^{n}\frac{\alpha^{k}}{k!}\frac{\partial^{k}m_{i,t}(\alpha=0)}{\partial\alpha^{k}}+\mathcal{O}(\alpha^{(n+1)}).$
(SI.4.216)

We solve this equation at $\alpha=1$. The approximation yields the nMF
equations when we ignore quadratic terms and higher, and the TAP equations
when ignoring third and higher order terms. However, in this case we will only
compute the terms in the nMF equation, since the second order yields marginals
that are complicated to evaluate. In the case of the first order term, we will
show how to estimate the corresponding marginals at the end of this
Supplementary Note.

The derivatives of this distribution with respect to $\alpha$ is given by

$\displaystyle\frac{\partial
P_{\alpha}(s_{i,t-1}|\mathbf{s}_{t-2})}{\partial\alpha}=$
$\displaystyle(s_{i,t-1}-\tanh{h}_{i,t-1}(\alpha))\Delta
h_{i,t-1}P_{\alpha}(s_{i,t-1}|\mathbf{s}_{t-2}).$ (SI.4.217)

Again, we define $h_{i,t-1}(\alpha)=\Theta_{i,t-1}+\alpha\,\Delta h_{i,t-1}$,
where $\Delta h_{i,t-1}=-\Theta_{i,t-1}+H_{i}+\sum_{j}J_{ij}s_{j,t-2}$
represents deviation from the independent model. The derivatives of each order
are derived as follows.

The zeroth order term is:

$\displaystyle m_{i,t}(\alpha)$
$\displaystyle=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\tanh
h_{i,t}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$
(SI.4.218) $\displaystyle m_{i,t}(\alpha=0)$
$\displaystyle=\sum_{\mathbf{s}_{t-1}}\tanh h_{i,t}Q(\mathbf{s}_{t-1}).$
(SI.4.219)

The first order terms are:

$\displaystyle\frac{\partial m_{i,t}(\alpha)}{\partial\alpha}$
$\displaystyle=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\tanh
h_{i,t}\sum_{k}(s_{k,t-1}-\tanh h_{k,t-1}(\alpha))\Delta
h_{k,t-1}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$
(SI.4.220) $\displaystyle\frac{\partial m_{i,t}(\alpha=0)}{\partial\alpha}$
$\displaystyle=\sum_{k}(-\Theta_{k,t-1}+H_{k}+\sum_{l}J_{kl}m_{l,t-2})\sum_{\mathbf{s}_{t-1}}\tanh
h_{i,t}(s_{k,t-1}-m_{k,t-1})Q(\mathbf{s}_{t-1}).$ (SI.4.221)

Since
$\left[\alpha(H_{k}-\Theta_{k,t-1}+\sum_{l}J_{kl}m_{l,t-2})\right]_{\alpha=1}=0+\left[\mathcal{O}(\alpha^{2})\right]_{\alpha=1}$
(3) we can ignore the first order terms, leading to the first order expansion

$\displaystyle m_{i,t}$
$\displaystyle\approx\Gamma_{i,t}+\left[\mathcal{O}(\alpha^{2})\right]_{\alpha=1},$
(SI.4.222)

where $\Gamma_{i,t}=\sum_{\mathbf{s}_{t-1}}\tanh h_{i,t}Q(\mathbf{s}_{t-1})$.
Moreover, since $Q(\mathbf{s}_{t-1})$ is an independent distribution, for a
large number of units, applying the central limit theorem, it can be
approximated by a Gaussian distribution (see last section of this
Supplementary Note).

### Equal-time correlations

When $i\neq k$, correlations in the system are calculated as

$C_{ik,t}(\alpha)=\sum_{\mathbf{s}_{t},\mathbf{s}_{t-1}}(s_{i,t}-m_{i,t}(\alpha))(s_{k,t}-m_{k,t}(\alpha))P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})\mathbf{s}_{t},\mathbf{s}_{t-1}(\mathbf{s}_{t},\mathbf{s}_{t-1}).$
(SI.4.223)

Again, we compute this using a Plefka expansion:

$C_{ik,t}(\alpha)=\sum_{n=0}^{m}\frac{\alpha^{n}}{n!}\frac{\partial^{n}C_{ik,t}(\alpha=0)}{\partial\alpha^{n}}+\mathcal{O}(\alpha^{(m+1)}).$
(SI.4.224)

Otherwise, when $i=j$, we have $C_{ii,t}=1-m_{i,t}^{2}$. The derivatives of
each order are derived as follows.

The zeroth order term:

$\displaystyle C_{ik,t}(\alpha)$
$\displaystyle=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}(\tanh
h_{i,t}-m_{i,t}(\alpha))(\tanh
h_{k,t}-m_{k,t}(\alpha))P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$
(SI.4.225) $\displaystyle C_{ik,t}(\alpha=0)$
$\displaystyle=\sum_{\mathbf{s}_{t-1}}\tanh h_{i,t}\tanh
h_{k,t}Q(\mathbf{s}_{t-1})-m_{i,t}m_{k,t}\equiv\Gamma_{ik,t}.$ (SI.4.226)

The first order term yields:

$\displaystyle\frac{\partial C_{ik,t}(\alpha)}{\partial\alpha}=$
$\displaystyle\sum_{s(t-1)}\Big{(}(\tanh h_{i,t}-m_{i,t}(\alpha))(\tanh
h_{k,t}-m_{k,t}(\alpha))\sum_{m}(s_{m,t-1}-\tanh h_{m,t-1}(\alpha))\Delta
h_{m,t-1}$ (SI.4.227) $\displaystyle-\frac{\partial
m_{i,t}(\alpha)}{\partial\alpha}(\tanh h_{k,t}-m_{k,t}(\alpha))-(\tanh
h_{i,t}-m_{i,t}(\alpha))\frac{\partial
m_{k,t}(\alpha)}{\partial\alpha}\Big{)}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$
(SI.4.228) $\displaystyle\frac{\partial C_{ik,t}(\alpha=0)}{\partial\alpha}=$
$\displaystyle-\sum_{m}(H_{m}-\Theta_{m,t-1}+\sum_{n}J_{mn}m_{n,t-2})(\sum_{\mathbf{s}_{t-1}}(\tanh
h_{k,t}-m_{k,t}$ (SI.4.229) $\displaystyle+\tanh
h_{i,t}-m_{i,t})(s_{m,t-1}-m_{m,t-1})Q(\mathbf{s}_{t-1}).$ (SI.4.230)

Since
$\left[\alpha(H_{m}-\Theta_{m,t-1}+\sum_{n}J_{mn}m_{n,t-2})\right]_{\alpha=1}=0+\left[\mathcal{O}(\alpha^{2})\right]_{\alpha=1}$
(3), the nMF equation is:

$C_{ik,t}\approx\Gamma_{ik,t}.$ (SI.4.231)

### Time-delayed correlations

Time-delayed correlations in the system are calculated as

$D_{il,t}(\alpha)=\sum_{\mathbf{s}_{t},\mathbf{s}_{t-1}}(s_{i,t}-m_{i,t}(\alpha))(s_{l,t-1}-m_{l,t-1}(\alpha))P_{\alpha}^{[t-1]}(\mathbf{s}_{t},\mathbf{s}_{t-1}).$
(SI.4.232)

We compute this using a Plefka expansion:

$D_{il,t}(\alpha)=\sum_{n=0}^{m}\frac{\alpha^{n}}{n!}\frac{\partial^{n}D_{il,t}(\alpha=0)}{\partial\alpha^{n}}+\mathcal{O}(\alpha^{(m+1)}).$
(SI.4.233)

The zeroth order terms are:

$\displaystyle D_{il,t}(\alpha)$
$\displaystyle=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}(\tanh
h_{i,t}-m_{i,t}(\alpha))(s_{l,t-1}-m_{l,t-1}(\alpha))P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$
(SI.4.234) $\displaystyle D_{il,t}(\alpha=0)$
$\displaystyle=\sum_{\mathbf{s}_{t-1}}(\tanh
h_{i,t}-m_{i,t})(s_{l,t-1}-m_{l,t-1})Q(\mathbf{s}_{t-1})\equiv\Gamma_{i,t}^{(l)}.$
(SI.4.235)

The first order terms are:

$\displaystyle\frac{\partial D_{il,t}(\alpha)}{\partial\alpha}$
$\displaystyle=\sum_{\mathbf{s}_{t-1},\mathbf{s}_{t-2}}\Big{(}(\tanh
h_{i,t}-m_{i,t}(\alpha))(s_{l,t-1}-m_{l,t-1}(\alpha))\sum_{k}(s_{k,t-1}-\tanh
h_{k,t-1}(\alpha))\Delta h_{k,t-1}$ (SI.4.236)
$\displaystyle\phantom{=}-\frac{\partial
m_{i,t}(\alpha))}{\partial\alpha}(s_{l,t-1}-m_{l,t-1}(\alpha))-(\tanh
h_{i,t}-m_{i,t}(\alpha))\frac{\partial