m_{l,t-1}(\alpha)}{\partial\alpha}\Big{)}P_{\alpha}(\mathbf{s}_{t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$
(SI.4.237) $\displaystyle\frac{\partial D_{il,t}(\alpha=0)}{\partial\alpha}$
$\displaystyle=\sum_{k}\sum_{\mathbf{s}_{t-1}}(\tanh
h_{i,t}-m_{i,t})(s_{l,t-1}-m_{l,t-1}(\alpha))(s_{k,t-1}$ (SI.4.238)
$\displaystyle\phantom{=}-m_{k,t-1}(\alpha))(-\Theta_{k,t-1}+H_{k}+\sum_{n}J_{kn}m_{n,t-2})Q(\mathbf{s}_{t-1}).$
(SI.4.239)

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

$D_{il,t}\approx\Gamma_{i,t}^{(l)}.$ (SI.4.240)

### Gaussian approximations

The integrals for computing the first order means and correlations can be
directly obtained applying the central limit theorem to approximate a set of
independent binary signals to a Gaussian distribution.

Thus we obtain

$\Gamma_{i,t}=\sum_{\mathbf{s}_{t-1}}\tanh
h_{i,t}Q(\mathbf{s}_{t-1})\approx\int\mathrm{D}_{x}\tanh[g_{i,t}+x\sqrt{\Delta_{i,t}}],$
(SI.4.241)

where $\mathrm{D}_{x}$ denotes an integral using a Gaussian distribution:

$\mathrm{D}_{x}=\frac{\mathrm{d}x}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}x^{2}\right).$
(SI.4.242)

The other parameters are $g_{i,t}=H_{i}+\sum_{j}J_{ij}m_{j,t-1}$ and
$\Delta_{i,t}=\mathrm{Var}[h_{i,t}]_{\alpha=0}=\sum_{j}J_{ij}^{2}(1-m_{j,t-1}^{2})$.

Similarly, we have

$\displaystyle\Gamma_{ik,t}$ $\displaystyle=\sum_{\mathbf{s}_{t-1}}\tanh
h_{i,t}\tanh h_{k,t}Q(\mathbf{s}_{t-1})$ (SI.4.243)
$\displaystyle\approx\int\mathrm{D}_{xy}^{\rho_{ik}}\tanh[g_{i,t}+x\sqrt{\Delta_{i,t}}]\tanh[g_{k,t}+y\sqrt{\Delta_{k,t}}]$
(SI.4.244)
$\displaystyle=\int\mathrm{D}_{x}\mathrm{D}_{y}\tanh[g_{i,t}+(x\sqrt{\frac{1+\rho_{ik}}{2}}+y\sqrt{\frac{1-\rho_{ik}}{2}})\sqrt{\Delta_{i,t}}]\tanh[g_{k,t}+(x\sqrt{\frac{1+\rho_{ik}}{2}}-y\sqrt{\frac{1-\rho_{ik}}{2}})\sqrt{\Delta_{k,t}}],$
(SI.4.245)

where $\mathrm{D}_{xy}^{\rho_{ik}}$ denotes an integral using a bivariate
Gaussian distribution:

$\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).$
(SI.4.246)

The other parameters are
$\rho_{ij}=\frac{\Delta_{ij,t}}{\sqrt{\Delta_{ii,t}\Delta_{kk,t}}}$, with
$\Delta_{ik,t}=\mathrm{Cov}(h_{i,t},h_{k,t})_{\alpha=0}=\sum_{j}J_{ij}J_{kj}(1-m_{j,t-1}^{2})$.
In the following step, we applied an orthogonal transformation to the previous
bivariate distribution, involving a change of variables
$x^{\prime}=x\sqrt{\frac{1+\rho_{ik}}{2}}+y\sqrt{\frac{1-\rho_{ik}}{2}}$ and
$y^{\prime}=x\sqrt{\frac{1+\rho_{ik}}{2}}-y\sqrt{\frac{1-\rho_{ik}}{2}}$ that
removes the coupling terms in the bivariate normal distribution.

The integrals of the Gaussian approximations of $\Gamma_{i,t}$ and
$\Gamma_{ik,t}$ are easy to compute. The problem arises when we deal with the
more complex terms $\Gamma_{i,t}^{(l)}$ in Supplementary Eq. SI.4.240, which
can be computed by similar approximations but multiplies the number of
integrals to be solved.

For large system sizes, terms obtained from the first and second order
expressions like $\Gamma_{i,t}^{(l)}$ (Supplementary Eq. SI.4.232) can be
obtained by assuming that the values of individual weights are small (e.g.
$J_{ij}=\mathcal{O}(1/N)$). We can compute

$\displaystyle
D_{il,t}\approx\Gamma_{i,t}^{(l)}=\sum_{\mathbf{s}_{t-1}}(s_{l,t}-m_{l,t})\tanh
h_{i,t}Q(\mathbf{s}_{t-1}),$ (SI.4.247)

where the constant term $m_{i,t}$ was removed from Supplementary Eq. SI.4.232
as it yields zero when averaging its product with $(s_{l,t}-m_{l,t})$.

With independent spins at $t-1$, we can approximate this quantity with the aid
of a bivariate Gaussin distribution of fields $h_{k,t},h_{i,t}$. Knowing that
$h_{k,t}-\langle h_{k,t}\rangle=\sum_{l}J_{kl}(s_{l,t-1}-m_{l,t-1})$, we can
approximate the following quantity

$\displaystyle\sum_{l}J_{kl}D_{il,t}\approx$
$\displaystyle\sum_{\mathbf{s}_{t-1}}\sum_{l}J_{kl}(s_{l,t-1}-m_{l,t-1})\tanh
h_{i,t}Q(\mathbf{s}_{t-1})$ (SI.4.248) $\displaystyle=$
$\displaystyle\sum_{\mathbf{s}_{t-1}}(h_{k,t}-\langle h_{k,t}\rangle)\tanh
h_{i,t}Q(\mathbf{s}_{t-1})$ (SI.4.249) $\displaystyle\approx$
$\displaystyle\int\mathrm{D}_{xy}^{\rho_{ik}}y\sqrt{\Delta_{k,t}}\tanh[g_{i,t}+x\sqrt{\Delta_{i,t}}].$
(SI.4.250)

where $\langle h_{k,t}\rangle=H_{k}+\sum_{l}m_{l,t-1}$.

As in [41], if we assume that $\rho_{ik}$ is
small,$\mathrm{D}_{xy}^{\rho_{ik}}$ can be approximated as:

$\mathrm{D}_{xy}^{\rho_{ik}}\approx\mathrm{D}_{x}\mathrm{D}_{y}(1+\rho_{ik}xy).$
(SI.4.251)

Then

$\displaystyle\sum_{l}J_{kl}D_{il,t}\approx$
$\displaystyle\int\mathrm{D}_{xy}^{\rho_{ik}}y\sqrt{\Delta_{j,t}}\tanh[g_{i,t}+x\sqrt{\Delta_{i,t}}]$
(SI.4.252) $\displaystyle\approx$
$\displaystyle\int\mathrm{D}_{x}\mathrm{D}_{y}(y+\rho_{ik}xy^{2})\sqrt{\Delta_{k,t}}\tanh[g_{i,t}+x\sqrt{\Delta_{i,t}}]$
(SI.4.253) $\displaystyle=$ $\displaystyle\rho_{ik}\sqrt{\Delta_{k,t}}\int
D_{x}x\tanh[g_{i,t}+x\sqrt{\Delta_{i,t}}]$ (SI.4.254) $\displaystyle=$
$\displaystyle\rho_{ik}\sqrt{\Delta_{i,t}\Delta_{k,t}}\int
D_{x}\left(1-\tanh^{2}[g_{i,t}+x\sqrt{\Delta_{i,t}}]\right),$ (SI.4.255)

where the last step was obtained by partial integration.

As $\rho_{ik}\sqrt{\Delta_{i,t}\Delta_{k,t}}=\sum_{jl}J_{ij}J_{kl}C_{jl,t-1}$
we have that

$\sum_{l}J_{kl}D_{il,t}\approx a_{i,t}\sum_{jl}J_{ij}J_{kl}C_{jl,t-1},$
(SI.4.256)

where

$a_{i,t}=\int D_{x}\left(1-\tanh^{2}[g_{i,t}+x\sqrt{\Delta_{i,t}}]\right).$
(SI.4.257)

Therefore, we have the following approximation for the delayed correlations,

$D_{il,t}\approx a_{i,t}\sum_{j}J_{ij}C_{jl,t-1}.$ (SI.4.258)

##  Supplementary Note 5 Pairwise Plefka expansions, Plefka2[$t$]

### Pairwise delayed spin correlations

Instead of a manifold of independent distributions, in this approximation we
consider a manifold $\mathcal{Q}$ with a pairwise probability distribution:

$Q(s_{i,t},s_{l,t-1})=Q(s_{i,t}|s_{l,t-1})Q(s_{l,t-1})=\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}},$
(SI.5.259)

where $\theta_{i,t}(s_{l,t-1})=\Theta_{i,t}+\Delta_{il,t}s_{l,t-1}$. Here
$Q(s_{l,t-1})$ is the independent probability distribution for $s_{l,t-1}$
computed as in 3, and $Q(s_{i,t}|s_{l,t-1})$ is a conditional probability
distribution we use to construct the pairwise probability distribution
$Q(s_{i,t},s_{l,t-1})$ using the chain rule.

We use this pairwise model to approximate the distribution

$P(s_{i,t},s_{l,t-1})=\sum_{\mathbf{s}_{\setminus
l,t-1}}P(s_{i,t},\mathbf{s}_{t-1})=\sum_{\mathbf{s}_{\setminus
l,t-1}}P(s_{i,t}|\mathbf{s}_{t-1})P(\mathbf{s}_{t-1}).$ (SI.5.260)

with $s_{\setminus l,t-1}$ containing all elements of $\mathbf{s}_{t-1}$
except $s_{l,t-1}$.

As in previous cases, we want to find an approximation of the probability
distribution at time $t$ that minimizes the relative entropy

$\displaystyle D(P(s_{i,t},s_{l,t-1})||Q(s_{i,t},s_{l,t-1}))$
$\displaystyle=\sum_{\begin{subarray}{c}s_{i,t}\\\
s_{l,t-1}\end{subarray}}P(s_{i,t},s_{l,t-1})\log\frac{P(s_{i,t},s_{l,t-1})}{Q(s_{i,t},s_{l,t-1})}.$
(SI.5.261)

Specifically, the mean-field approximation that minimizes the relative entropy
is the one that satisfies

$\displaystyle\frac{\partial
D(P(s_{i,t},s_{l,t-1})||Q(s_{i,t},s_{l,t-1}))}{\partial\Theta_{i,t}}$
$\displaystyle=-\sum_{\begin{subarray}{c}s_{i,t}\\\
s_{l,t-1}\end{subarray}}\Big{(}s_{i,t}-\tanh[\theta_{i,t}(s_{l,t-1})]\Big{)}P(s_{i,t},s_{l,t-1})$
(SI.5.262) $\displaystyle=\langle s_{i,t}\rangle_{Q}-\langle
s_{i,t}\rangle_{P}=0,$ (SI.5.263) $\displaystyle\frac{\partial
D(P(s_{i,t},s_{l,t-1})||Q(s_{i,t},s_{l,t-1}))}{\partial\Delta_{il,t}}$
$\displaystyle=-\sum_{\begin{subarray}{c}s_{i,t}\\\
s_{l,t-1}\end{subarray}}\Big{(}s_{i,t}s_{l,t-1}-\tanh[\theta_{i,t}(s_{l,t-1})]s_{l,t-1}\Big{)}P(s_{i,t},s_{l,t-1})$
(SI.5.264)
$\displaystyle=\langle\tanh[\theta_{i,t}(s_{l,t-1})]s_{l,t-1}\rangle_{P}-\langle
s_{i,t}s_{l,t-1}\rangle_{P}$ (SI.5.265) $\displaystyle=\langle
s_{i,t}s_{l,t-1}\rangle_{(Q\cdot P)}-\langle s_{i,t}s_{l,t-1}\rangle_{P}=0,$
(SI.5.266)

where $\langle s_{i,t}s_{l,t-1}\rangle_{(Q\cdot
P)}=\sum_{s_{i,t},s_{l,t-1}}s_{i,t}s_{l,t-1}Q(s_{i,t}|s_{l,t-1})P(s_{l,t-1})$,
and we used the equivalence
$\tanh[\theta_{i,t}(s_{l,t-1})]=\sum_{s_{i,t}}s_{i,t}\frac{\mathrm{e}^{s_{i,t}\theta_{i,t}(s_{l,t-1})}}{2\cosh\theta_{i,t}(s_{l,t-1})}$.

This equation states that the closest factorized model has its first and
second order moments equal to the first moments of the target distribution
$P$. That is, $\langle s_{i,t}\rangle_{Q}=\langle s_{i,t}\rangle_{P}=m_{i,t}$
and $\langle s_{i,t}s_{l,t-1}\rangle_{(Q\cdot P)}=\langle
s_{i,t}s_{l,t-1}\rangle_{P}=D_{il,t}+m_{i,t}m_{l,t-1}$. This is equivalent to
having the marginalized distribution for spins $i$ and $l$ equal to the model
$Q$, i.e., $P(s_{i,t},s_{l,t-1})=Q(s_{i,t},s_{l,t-1})$. If we assume that the
distribution $P$ is close to $Q$, we can compute $P(s_{i,t},s_{l,t-1})$ as an
expansion with respect to $\alpha$ of the probability distribution:

$P_{\alpha}(s_{i,t},s_{l,t-1})=\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})P(\mathbf{s}_{\setminus
l,t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$ (SI.5.267)

with

$\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.5.268) $\displaystyle h_{i,t}(\alpha)$
$\displaystyle=(1-\alpha)\theta_{i,t}(s_{l,t-1})+\alpha(H_{i}+\sum_{j}J_{ij}s_{j,t-1}),$
(SI.5.269)

and

$\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)},$
(SI.5.270) $\displaystyle h_{l,t-1}(\alpha)$
$\displaystyle=(1-\alpha)\Theta_{l,t-1}+\alpha(H_{l}+\sum_{n}J_{ln}s_{n,t-2}).$
(SI.5.271)

When $\alpha$ is set to zero,
$P_{\alpha=0}(s_{i,t},s_{l,t-1})=Q(s_{i,t},s_{l,t-1})$, whereas when
$\alpha=1$, $P_{\alpha=1}(s_{i,t},s_{l,t-1})=P(s_{i,t},s_{l,t-1})$. We
approximate the values of $\Theta_{i,t},\Theta_{l,t-1},\Delta_{il,t}$ as
follows:

$\displaystyle\frac{\partial P_{\alpha}(s_{i,t},s_{l,t-1})}{\partial\alpha}=$
$\displaystyle\sum_{\begin{subarray}{c}\mathbf{s}_{\setminus l,t-1}\\\
\mathbf{s}_{t-1}\end{subarray}}\Big{(}(s_{i,t}-\tanh
h_{i,t}(\alpha))(-\theta_{i,t}(s_{l,t-1})+H_{i}+\sum_{j}J_{ij}s_{j,t-1})$
(SI.5.272) $\displaystyle+(s_{l,t-1}-\tanh
h_{l,t-1}(\alpha))(-\Theta_{l,t-1}+H_{l}+\sum_{n}J_{ln}s_{n,t-2})\Big{)}$
(SI.5.273) $\displaystyle\cdot
P_{\alpha}(s_{i,t}|\mathbf{s}_{t-1})P_{\alpha}(s_{l,t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{\setminus
l,t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$ (SI.5.274)
$\displaystyle\frac{\partial
P_{\alpha}(s_{i,t},s_{l,t-1})}{\partial\alpha}\Big{\rvert}_{\alpha=0}=$
$\displaystyle\Big{(}(s_{i,t}-\tanh\theta_{i,t}(s_{l,t-1}))(-\theta_{i,t}(s_{l,t-1})+J_{il}(s_{l,t-1}-m_{l,t-1})+H_{i}+\sum_{j}J_{ij}m_{j,t-1})$
(SI.5.275)
$\displaystyle+(s_{l,t-1}-m_{l,t-1})(-\Theta_{l,t-1}+H_{l}+\sum_{n}J_{ln}m_{n,t-2})\Big{)}Q(s_{i,t},s_{l,t-1}).$
(SI.5.276)

From here, we obtain the nMF equations:

$\displaystyle\theta_{i,t}(s_{l,t-1})$ $\displaystyle\approx
H_{i}+\sum_{j}J_{ij}m_{j,t-1}+J_{il}(s_{l,t-1}-m_{l,t-1}),$ (SI.5.277)
$\displaystyle\Theta_{l,t-1}$ $\displaystyle\approx
H_{l}+\sum_{n}J_{ln}m_{n,t-1}.$ (SI.5.278)

The second order expressions are obtained by expanding:

$\displaystyle\frac{\partial^{2}P_{\alpha}(s_{i,t},s_{l,t-1})}{\partial\alpha^{2}}=$
$\displaystyle\sum_{\begin{subarray}{c}\mathbf{s}_{\setminus l,t-1}\\\
\mathbf{s}_{t-1}\end{subarray}}\Big{(}-(1-\tanh^{2}h_{i,t}(\alpha))(-\theta_{i,t}(s_{l,t-1})+H_{i}+\sum_{j}J_{ij}s_{j,t-1})^{2}$
(SI.5.279)
$\displaystyle-(1-\tanh^{2}h_{l,t-1}(\alpha))(-\Theta_{l,t-1}+H_{l}+\sum_{n}J_{ln}s_{n,t-2})^{2}$
(SI.5.280) $\displaystyle+\big{(}(s_{i,t}-\tanh
h_{i,t}(\alpha))(-\theta_{i,t}(s_{l,t-1})+H_{i}+\sum_{j}J_{ij}s_{j,t-1})$
(SI.5.281) $\displaystyle+(s_{l,t-1}-\tanh
h_{l,t-1}(\alpha))(-\Theta_{l,t-1}+H_{l}+\sum_{n}J_{ln}s_{n,t-2})\big{)}^{2}\Big{)}$
(SI.5.282) $\displaystyle\cdot
P_{\alpha}(s_{i,t}|\mathbf{s}_{t-1})P_{\alpha}(s_{l,t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{\setminus
l,t-1}|\mathbf{s}_{t-2})P(\mathbf{s}_{t-2}),$ (SI.5.283)
$\displaystyle\frac{\partial^{2}P_{\alpha}(s_{i,t},s_{l,t-1})}{\partial\alpha^{2}}\Big{\rvert}_{\alpha=0}=$
$\displaystyle\Big{(}-2\tanh\theta_{i,t}(s_{l,t-1})(s_{i,t}-\tanh\theta_{i,t}(s_{l,t-1}))\big{(}W_{i,t}^{2}+V_{ii,t}\big{)}$
(SI.5.284)
$\displaystyle-2m_{l,t-1}(s_{l,t-1}-m_{l,t-1})\big{(}W_{l,t-1}^{2}+V_{ll,t-1}\big{)}$
(SI.5.285)
$\displaystyle+2(s_{i,t}-\tanh\theta_{i,t}(s_{l,t-1}))(s_{l,t-1}-m_{l,t-1})\big{(}W_{i,t}W_{l,t-1}$
(SI.5.286) $\displaystyle+\sum_{j\neq
l,n}J_{ij}J_{ln}D_{jn,t-1}\big{)}\Big{)}Q(s_{i,t},s_{l,t-1}),$ (SI.5.287)

where $V_{ii,t}=\sum_{j\neq l,n\neq l}J_{ij}J_{in}C_{jn,t-1}$,
$V_{ll,t-1}=\sum_{mn}J_{lm}J_{ln}C_{mn,t-2}$,
$W_{i,t}=-\theta_{i,t}(s_{l,t-1})+H_{i}+\sum_{j}J_{ij}m_{j,t-1}+J_{il}(s_{l,t-1}-m_{l,t-1})$
and $W_{l,t-1}=-\Theta_{l,t-1}+H_{l}+\sum_{n}J_{ln}m_{n,t-1}$. Note that the
first term in the second equation comes from the combination of the term in
the first line with the squared term in the third line of the first equation.
Similarly, the second term in the second equation derives from the combination
of the second line and the squared term in the fourth line of the first
equation.

We know that
$\big{[}\alpha^{2}W_{i,t}^{2}\big{]}_{\alpha=1}=0+\big{[}\mathcal{O}(\alpha^{4})\big{]}_{\alpha=1}$,
that
$\big{[}\alpha^{2}W_{l,t-1}^{2}\big{]}_{\alpha=1}=0+\big{[}\mathcal{O}(\alpha^{4})\big{]}_{\alpha=1}$
and
$\big{[}\alpha^{2}W_{i,t}W_{l,t-1}\big{]}_{\alpha=1}=0+\big{[}\mathcal{O}(\alpha^{4})\big{]}_{\alpha=1}$.
This allows us to dismiss terms of the equations above. Then, grouping
together the terms that multiply $(s_{i,t}-\tanh\theta_{i,t}(s_{l,t-1}))$ and
those that do not, we obtain the TAP equations:

$\displaystyle\theta_{i,t}(s_{l,t-1})=$ $\displaystyle
H_{i}+\sum_{j}J_{ij}m_{j,t-1}+\left(J_{il}+\sum_{j\neq
l,n}J_{ij}J_{ln}D_{jn,t-1}\right)(s_{l,t-1}-m_{l,t-1})$ (SI.5.288)
$\displaystyle-\tanh\left[\theta_{i,t}(s_{l,t-1})\right]\sum_{jn\neq
l}J_{ij}J_{in}C_{jn,t-1}+\big{[}\mathcal{O}(\alpha^{2})\big{]}_{\alpha=1},$
(SI.5.289) $\displaystyle\Theta_{l,t-1}=$ $\displaystyle
H_{l}+\sum_{n}J_{ln}m_{n,t-2}-m_{l,t-1}\sum_{mn}J_{lm}J_{ln}C_{mn,t-2}+\big{[}\mathcal{O}(\alpha^{2})\big{]}_{\alpha=1}.$
(SI.5.290)

With these parameters we can compute the TAP approximations of the system’s
sufficient statistics as follows:

$\displaystyle m_{i,t}\approx$
$\displaystyle\sum_{s_{l,t-1}}\tanh[\theta_{i,t}(s_{l,t-1})]Q(s_{l,t-1}),$
(SI.5.291) $\displaystyle m_{l,t-1}\approx$
$\displaystyle\tanh\Theta_{l,t-1},$ (SI.5.292) $\displaystyle D_{il,t}\approx$
$\displaystyle\sum_{s_{l,t-1}}(\tanh[\theta_{i,t}(s_{l,t-1})]-m_{i,t})(s_{l,t-1}-m_{l,t-1})Q(s_{l,t-1}),$
(SI.5.293)

where $Q(s_{l,t-1})=P(s_{l,t-1})=\frac{1+s_{l,t-1}m_{l,t-1}}{2}$ is the
factorised distribution for $s_{l,t-1}$.

### Pairwise equal-time spin correlations

We can use a similar approximation in the case above to compute equal-time
spin correlations. We do this by substituting $s_{l,t-1}$ in the approximation
above by $s_{k,t}$. As $s_{i,t}$, $s_{k,t}$ are conditionally independent, we
can assume that $s_{k,t}$ is computed first and that then $s_{i,t}$ is
computed conditioned on $s_{k,t}$, (i.e., consider
$P(s_{i,t}|s_{k,t},\mathbf{s}_{t-1})=P(s_{i,t}|\mathbf{s}_{t-1})$). In this
approximation we consider a manifold $\mathcal{Q}$ with a pairwise probability
distribution:

$Q(s_{i,t},s_{k,t})=Q(s_{i,t}|s_{k,t})Q(s_{k,t})=\frac{\mathrm{e}^{s_{i,t}\theta_{i,t}(s_{k,t})}}{2\cosh\theta_{i,t}(s_{k,t})}\frac{\mathrm{e}^{s_{k,t}\Theta_{k,t}}}{2\cosh\Theta_{k,t}},$
(SI.5.294)

where $\theta_{i,t}(s_{k,t})=\Theta_{i,t}+\Delta_{ik,t}s_{k,t}$. Here
$Q(s_{k,t})$ is the independent probability distribution for $s_{k,t}$
computed as in 3, and $Q(s_{i,t}|s_{k,t})$ is a conditional probability
distribution we use to construct the pairwise probability distribution
$Q(s_{i,t},s_{k,t})$ using the chain rule.

We use this pairwise model to approximate the distribution

$P(s_{i,t},s_{k,t})=\sum_{\mathbf{s}_{t-1}}P(s_{i,t},s_{k,t},\mathbf{s}_{t-1})=\sum_{\mathbf{s}_{t-1}}P(s_{i,t}|\mathbf{s}_{t-1})P(s_{k,t},\mathbf{s}_{t-1}).$
(SI.5.295)

As in previous cases, we want to find an approximation of the probability
distribution at time $t$ that minimizes the relative entropy

$\displaystyle D(P(s_{i,t},s_{k,t})||Q(s_{i,t},s_{k,t}))$
$\displaystyle=\sum_{\begin{subarray}{c}s_{i,t}\\\
s_{k,t}\end{subarray}}P(s_{i,t},s_{k,t})\log\frac{P(s_{i,t},s_{k,t})}{Q(s_{i,t},s_{k,t})}.$
(SI.5.296)

Specifically, the mean-field approximation that minimizes the relative entropy
is the one that satisfies

$\displaystyle\frac{\partial
D(P(s_{i,t},s_{k,t})||Q(s_{i,t},s_{k,t}))}{\partial\Theta_{i,t}}$
$\displaystyle=-\sum_{\begin{subarray}{c}s_{i,t}\\\
s_{k,t}\end{subarray}}\Big{(}s_{i,t}-\tanh[\theta_{i,t}(s_{k,t})])\Big{)}P(s_{i,t},s_{k,t})$
(SI.5.297) $\displaystyle=\langle s_{i,t}\rangle_{Q}-\langle
s_{i,t}\rangle_{P}=0,$ (SI.5.298) $\displaystyle\frac{\partial
D(P(s_{i,t},s_{k,t})||Q(s_{i,t},s_{k,t}))}{\partial\Delta_{ik,t}}$
$\displaystyle=-\sum_{\begin{subarray}{c}s_{i,t}\\\
s_{k,t}\end{subarray}}\Big{(}s_{i,t}s_{k,t}-\tanh[\theta_{i,t}(s_{k,t})])s_{k,t}\Big{)}P(s_{i,t},s_{k,t})$
(SI.5.299) $\displaystyle=\langle s_{i,t}s_{k,t}\rangle_{(Q\cdot P)}-\langle
s_{i,t}s_{k,t}\rangle_{P}=0,$ (SI.5.300)

where $\langle s_{i,t}s_{k,t}\rangle_{(Q\cdot
P)}=\sum_{s_{i,t},s_{k,t}}Q(s_{i,t}|s_{k,t})P(s_{k,t})$.

This equation states that the closest factorized model has its first moments
equal to the first and second order moments of the target distribution $P$.
That is, $\langle s_{i,t}\rangle_{Q}=\langle s_{i,t}\rangle_{P}=m_{i,t}$ and
$\langle s_{i,t}s_{k,t}\rangle_{(Q\cdot P)}=\langle
s_{i,t}s_{k,t}\rangle_{P}=C_{ik,t}+m_{i,t}m_{k,t}$. This is equivalent to
having the marginalized distribution for spins $i$ and $l$ equal to the model
$Q$, i.e. $P(s_{i,t},s_{k,t})=Q(s_{i,t},s_{k,t})$. If we assume that the
distribution $P$ is close to $Q$, we can compute $P(s_{i,t},s_{k,t})$ as an
expansion with respect to $\alpha$ of the probability distribution:

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

with

$\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)}$ (SI.5.302) $\displaystyle h_{i,t}(\alpha)$
$\displaystyle=(1-\alpha)\theta_{i,t}(s_{k,t})+\alpha(H_{i}+\sum_{j}J_{ij}s_{j,t-1}),$
(SI.5.303)

and

$\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)},$
(SI.5.304) $\displaystyle h_{k,t}(\alpha)$
$\displaystyle=(1-\alpha)\Theta_{k,t}+\alpha(H_{k}+\sum_{l}J_{kl}s_{l,t-1}).$
(SI.5.305)

When $\alpha$ is set to zero,
$P_{\alpha=0}(s_{i,t},s_{k,t})=Q(s_{i,t},s_{k,t})$, whereas when $\alpha=1$,
$P_{\alpha=1}(s_{i,t},s_{k,t})=P(s_{i,t},s_{k,t})$. We approximate the values
of $\Theta_{i,t},\Theta_{k,t},\Delta_{ik,t}$ as follows:

$\displaystyle\frac{\partial P_{\alpha}(s_{i,t},s_{k,t})}{\partial\alpha}=$
$\displaystyle\sum_{\mathbf{s}_{t-1}}\Big{(}(s_{i,t}-\tanh
h_{i,t}(\alpha))(-\theta_{i,t}(s_{k,t})+H_{i}+\sum_{j}J_{ij}s_{j,t-1})$
(SI.5.306) $\displaystyle+(s_{k,t}-\tanh
h_{k,t}(\alpha))(-\Theta_{k,t}+H_{k}+\sum_{l}J_{kl}s_{l,t-1})\Big{)}$
(SI.5.307) $\displaystyle\cdot
P_{\alpha}(s_{i,t}|s_{k,t},\mathbf{s}_{t-1})P_{\alpha}(s_{k,t}|\mathbf{s}_{t-1})P(\mathbf{s}_{t-1}),$
(SI.5.308) $\displaystyle\frac{\partial
P_{\alpha}(s_{i,t},s_{k,t})}{\partial\alpha}\Big{\rvert}_{\alpha=0}=$
$\displaystyle\Big{(}(s_{i,t}-\tanh\theta_{i,t}(s_{k,t}))(-\theta_{i,t}(s_{k,t})+H_{i}+\sum_{j}J_{ij}m_{j,t-1})$
(SI.5.309)
$\displaystyle+(s_{k,t}-m_{k,t})(-\Theta_{k,t}+H_{k}+\sum_{l}J_{kl}m_{l,t-1})\Big{)}Q(s_{i,t},s_{k,t}).$
(SI.5.310)

From here, we obtain the nMF expansions:

$\displaystyle\theta_{i,t}(s_{k,t})$
$\displaystyle=H_{i}+\sum_{j}J_{ij}m_{j,t-1}+\big{[}\mathcal{O}(\alpha^{1})\big{]}_{\alpha=1},$
(SI.5.311) $\displaystyle\Theta_{k,t}$
$\displaystyle=H_{k}+\sum_{l}J_{kl}m_{l,t-1}+\big{[}\mathcal{O}(\alpha^{1})\big{]}_{\alpha=1}.$
(SI.5.312)

The second order expressions are obtained by expanding:

$\displaystyle\frac{\partial^{2}P_{\alpha}(s_{i,t},s_{k,t})}{\partial\alpha^{2}}=$
$\displaystyle\sum_{\mathbf{s}_{t-1}}\Big{(}-(1-\tanh^{2}h_{i,t}(\alpha))(-\theta_{i,t}(s_{k,t})+H_{i}+\sum_{j}J_{ij}s_{j,t-1})^{2}$
(SI.5.313)
$\displaystyle-(1-\tanh^{2}h_{k,t}(\alpha))(-\Theta_{k,t}+H_{k}+\sum_{l}J_{kl}s_{l,t-1})^{2}$
(SI.5.314) $\displaystyle+\big{(}(s_{i,t}-\tanh
h_{i,t}(\alpha))(-\theta_{i,t}(s_{k,t})+H_{i}+\sum_{j}J_{ij}s_{j,t-1})$
(SI.5.315) $\displaystyle+(s_{k,t}-\tanh
h_{k,t}(\alpha))(-\Theta_{k,t}+H_{k}+\sum_{l}J_{kl}s_{n,t-1})\big{)}^{2}\Big{)}$
(SI.5.316) $\displaystyle\cdot
P_{\alpha}(s_{i,t}|s_{k,t},\mathbf{s}_{t-1})P_{\alpha}(s_{k,t}|\mathbf{s}_{t-1})P(\mathbf{s}_{t-1}),$
(SI.5.317)
$\displaystyle\frac{\partial^{2}P_{\alpha}(s_{i,t},s_{k,t})}{\partial\alpha^{2}}\Big{\rvert}_{\alpha=0}=$
$\displaystyle\Big{(}-2\tanh\theta_{i,t}(s_{k,t})(s_{i,t}-\tanh\theta_{i,t}(s_{k,t}))\big{(}W_{i,t}^{2}+V_{ii,t}\big{)}$
(SI.5.318)
$\displaystyle-2m_{k,t}(s_{k,t}-m_{k,t})\big{(}W_{k,t}^{2}+V_{kk,t}\big{)}$
(SI.5.319)
$\displaystyle+2(s_{i,t}-\tanh\theta_{i,t}(s_{k,t}))(s_{k,t}-m_{k,t})\big{(}W_{i,t}W_{k,t}+V_{ik,t}\big{)}\Big{)}Q(s_{i,t},s_{k,t}),$
(SI.5.320)

where $V_{ik,t}=\sum_{jl}J_{ij}J_{kl}C_{jl,t-1}$ and
$W_{i,t}=-\theta_{i,t}+H_{i}+\sum_{j}J_{ij}m_{j,t-1}$. Again, note that the
first term in the second equation comes from the combination of the term in
the first line with the squared term in the third line of the first equation.
Similarly, the second term in the second equation derives from the combination
of the second line and the squared term in the fourth line of the first
equation.

We know that
$\big{[}\alpha^{2}W_{i,t}^{2}\big{]}_{\alpha=1}=0+\big{[}\mathcal{O}(\alpha^{4})\big{]}_{\alpha=1}$,
that
$\big{[}\alpha^{2}W_{k,t}^{2}\big{]}_{\alpha=1}=0+\big{[}\mathcal{O}(\alpha^{4})\big{]}_{\alpha=1}$
and
$\big{[}\alpha^{2}W_{i,t}W_{k,t}\big{]}_{\alpha=1}=0+\big{[}\mathcal{O}(\alpha^{4})\big{]}_{\alpha=1}$.
This allows us to dismiss terms of the equations above. Again, grouping terms
that contain $(s_{i,t}-\tanh\theta_{i,t}(s_{l,t-1}))$ and those that not we
obtain the follwoing Plefka expansions for the TAP approximations:

$\displaystyle\theta_{i,t}(s_{k,t})=$ $\displaystyle
H_{i}+\sum_{j}J_{ij}m_{j,t-1}+(s_{k,t}-m_{k,t})\sum_{jl}J_{ij}J_{kl}C_{jl,t-1}$
(SI.5.321)
$\displaystyle-\tanh\left[\theta_{i,t}(s_{k,t})\right]\sum_{jl}J_{ij}J_{il}C_{jl,t-1}+\big{[}\mathcal{O}(\alpha^{2})\big{]}_{\alpha=1},$
(SI.5.322) $\displaystyle\Theta_{k,t}=$ $\displaystyle
H_{k}+\sum_{l}J_{kl}m_{l,t-1}-m_{k,t}\sum_{jl}J_{kj}J_{kl}C_{jl,t-1}+\big{[}\mathcal{O}(\alpha^{2})\big{]}_{\alpha=1}.$
(SI.5.323)

With these parameters we can the TAP approximations of the system’s statistics
as follows:

$\displaystyle m_{i,t}\approx$
$\displaystyle\sum_{s_{k,t}}\tanh[\theta_{i,t}(s_{k,t})]Q(s_{k,t}),$
(SI.5.324) $\displaystyle m_{k,t}\approx$ $\displaystyle\tanh\Theta_{k,t},$
(SI.5.325) $\displaystyle C_{ik,t}\approx$
$\displaystyle\sum_{s_{k,t}}(\tanh[\theta_{i,t}(s_{k,t})]-m_{i,t})(s_{k,t}-m_{k,t})Q(s_{k,t}),$
(SI.5.326)

where $Q(s_{k,t})=P(s_{k,t})=\frac{1+s_{k,t}m_{k,t}}{2}$ is the factorised
probability of unit $s_{k,t}$.

##  Supplementary Note 6 Solution of the asymmetric kinetic Sherrington-
Kirkpatrick model

The infinite kinetic Ising model with Gaussian couplings used in this article
is generally referred in its symmetric version as the Sherrington-Kirkpatrick
(SK) model. The SK model behaviour is well studied in statistical mechanics,
and its solution can be obtained using the replica trick [60]. In the
symmetric case, dynamics can be represented as a bipartite network that can
also be solved using the replica trick [61]. Here, we extend the solution to
the kinetic, asymmetric version of the model. As this model does not have an
equilibrium distribution or a free energy defined in classical terms as the SK
model, we need to recur to an dynamical equivalent in the form of a generating
functional.

### Generating functional

We start with a kinetic Ising model

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

where $\beta$ is the inverse temperature, $H_{i,t}$ are time-varying fields
and $J_{ij}$ couplings. The results apply similarly to fixed fields.

The probability of a specific trajectory $\mathbf{s}_{0:t}$, is defined as

$\displaystyle P(\mathbf{s}_{0:t})=$
$\displaystyle\prod_{u=1}^{t}P(\mathbf{s}_{u}|\mathbf{s}_{u-1})P(\mathbf{s}_{0})$
(SI.6.329) $\displaystyle=$
$\displaystyle\exp\Bigg{(}\sum_{u}\left(\sum_{i}s_{i,u}\left(\beta
H_{i,u}+\beta\sum_{j}J_{ij}s_{j,u-1}\right)\right)$ (SI.6.330)
$\displaystyle-\sum_{i,u}\log
2\cosh\left(\beta\left(H_{i,u}+\sum_{j}J_{ij}s_{j,u-1}\right)\right)\Bigg{)}P(\mathbf{s}_{0}).$
(SI.6.331)

Instead of a partition function, the distribution of the trajectories of the
asymmetric SK model can be described by defining a generating functional or a
dynamical partition function:

$\displaystyle Z_{t}(\mathbf{g})=$
$\displaystyle\sum_{\mathbf{s}_{0:t}}\exp{\left(\sum_{iu}g_{i,u}s_{i,u}\right)}P(\mathbf{s}_{0:t})$
(SI.6.332) $\displaystyle=$
$\displaystyle\sum_{\mathbf{s}_{0:t}}\exp\Bigg{(}\sum_{u}\left(\sum_{i}s_{i,u}\left(\beta
H_{i,u}+g_{i,u}+\beta\sum_{j}J_{ij}s_{j,u-1}\right)\right)$ (SI.6.333)
$\displaystyle-\sum_{i,u}\log
2\cosh\left(\beta\left(H_{i,u}+\sum_{j}J_{ij}s_{j,u-1}\right)\right)\Bigg{)}P(\mathbf{s}_{0}).$
(SI.6.334)

Note that $Z_{t}(\mathbf{0})=1$. For simplicity we will assume an initial
distribution with a single possible state $\mathbf{s}_{0}$, i.e.
$P(\mathbf{s}_{0}^{\prime})=\delta_{\mathbf{s}_{0}^{\prime},\mathbf{s}_{0}}$,
which allows us to drop the $P(\mathbf{s}_{0})$ term in the equation above. In
the $t\to\infty$ limit, the logarithmic dynamical partition function converges
to the large deviation function

$\lim_{t\to\infty}\frac{1}{t}\log Z_{t}(\mathbf{g})=\varphi(\mathbf{g}),$
(SI.6.335)

which plays the role of a free-energy function for trajectories [62].

The generating functional can be used to compute moments of the system as:

$\displaystyle m_{i,u}=$
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial
Z_{t}(\mathbf{g})}{\partial
g_{i,u}}=\lim_{\mathbf{g}\to\mathbf{0}}\left\langle
s_{i,u}\right\rangle_{\mathbf{g}}=\left\langle s_{i,u}\right\rangle,$
(SI.6.336) $\displaystyle R_{ij,uv}=$
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial^{2}Z_{t}(\mathbf{g})}{\partial
g_{i,u}\partial g_{j,v}}=\lim_{\mathbf{g}\to\mathbf{0}}\left\langle
s_{i,u}s_{j,v}\right\rangle_{\mathbf{g}}=\left\langle
s_{i,u}s_{j,v}\right\rangle,$ (SI.6.337)

where the brackets are defined as

$\displaystyle\left\langle f(\mathbf{s})\right\rangle_{\mathbf{g}}=$
$\displaystyle\sum_{\mathbf{s}_{0:t}}f(\mathbf{s})\exp{\left(\sum_{iu}g_{i,u}s_{i,u}\right)}P(\mathbf{s}_{0:t}),$
(SI.6.338) $\displaystyle\left\langle f(\mathbf{s})\right\rangle=$
$\displaystyle\sum_{\mathbf{s}_{0:t}}f(\mathbf{s})P(\mathbf{s}_{0:t}).$
(SI.6.339)

It also allows us to derive identities that will be helpful in eliminating
spurious solutions

$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial
Z_{t}(\mathbf{g})}{\partial H_{i,u}}=\beta\left(m_{i,u}-m_{i,u}\right)=$
$\displaystyle 0,$ (SI.6.340)
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial^{2}Z_{t}(\mathbf{g})}{\partial
H_{i,u}\partial H_{j,v}}=\beta\left(\frac{\partial m_{i,u}}{\partial
H_{j,v}}-\frac{\partial m_{i,u}}{\partial H_{j,v}}\right)=$ $\displaystyle 0.$
(SI.6.341)

### Path integral

In the asymmetric SK model, the couplings $J_{ij}$ are quenched variables with
a Gaussian distribution function

$P(J_{ij})=\frac{1}{\sqrt{2\pi
J_{\sigma}^{2}/N}}\exp{\left(-\frac{1}{2J_{\sigma}^{2}/N}\left(J_{ij}-\frac{J_{0}}{N}\right)^{2}\right)},$
(SI.6.342)

where the mean and the variance are proportional to $1/N$.

In the thermodynamic limit, one can study the system by computing the
configurational average

$\left[Z_{t}(\mathbf{g})\right]=\int\prod_{ij}\mathrm{d}J_{ij}P(J_{ij})Z_{T}(\mathbf{g}).$
(SI.6.343)

The quenched average can be solved using path integral methods, by inserting
an appropriate delta integral for the effective fields of each unit

$1=\int\mathrm{d}\bm{\theta}\prod_{i,u}\delta(\theta_{i,u}-\beta
H_{i,u}-\beta\sum_{j}J_{ij}s_{j,u-1})=\frac{1}{\left(2\pi\right)^{Nt}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\exp\left(\sum_{iu}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,t}(\theta_{i,t}-\beta
H_{i,u}-\beta\sum_{j}J_{ij}s_{j.t-1})\right).$ (SI.6.344)

where $\bm{\theta}$ represent the effective fields of units and
$\bm{\hat{\theta}}$ are conjugates of these fields.

By inserting the above equation, the configurational average is written as

$\displaystyle\left[Z_{t}(\mathbf{g})\right]=$
$\displaystyle\frac{1}{\left(2\pi\right)^{Nt}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\prod_{ij}\mathrm{d}J_{ij}P(J_{ij})\sum_{\mathbf{s}_{1:t}}\exp\Bigg{(}\sum_{iu}s_{i,u}(g_{i,u}+\theta_{i,u})-\sum_{iu}\log
2\cosh\left(\theta_{i,u}\right)$ (SI.6.345)
$\displaystyle+\sum_{iu}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}(\theta_{i,u}-\beta
H_{i,u}-\beta\sum_{j}J_{ij}s_{j,u-1})\Bigg{)}$ (SI.6.346) $\displaystyle=$
$\displaystyle\frac{1}{\left(2\pi\right)^{Nt}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\sum_{\mathbf{s}_{1:t}}\exp\Bigg{(}\sum_{iu}s_{i,u}(g_{i,u}+\theta_{i,u})-\sum_{iu}\log
2\cosh\left(\theta_{i,u}\right)+\sum_{iu}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}(\theta_{i,u}-\beta H_{i,u})$ (SI.6.347)
$\displaystyle-\sum_{iu}\mathrm{i}\mkern 1.0mu\frac{\beta
J_{0}}{N}\hat{\theta}_{i,u}\sum_{j}s_{j,u-1}+\sum_{ij}\frac{\beta^{2}J_{\sigma}^{2}}{2N}\left(\mathrm{i}\mkern
1.0mu\sum_{u}\hat{\theta}_{i,u}s_{j,u-1}\right)^{2}\Bigg{)}$ (SI.6.348)
$\displaystyle=$
$\displaystyle\frac{1}{\left(2\pi\right)^{Nt}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\sum_{\mathbf{s}_{1:t}}\exp\Bigg{(}\sum_{iu}s_{i,u}(g_{i,u}+\theta_{i,u})-\sum_{iu}\log
2\cosh\left(\theta_{i,u}\right)+\sum_{iu}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}(\theta_{i,u}-\beta H_{i,u})$ (SI.6.349)
$\displaystyle-N\beta J_{0}\sum_{u}\frac{1}{N}\sum_{i}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\frac{1}{N}\sum_{j}s_{j,u-1}+N\frac{\beta^{2}J_{\sigma}^{2}}{2}\sum_{uv}\frac{1}{N}\sum_{i}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,v}\frac{1}{N}\sum_{j}s_{j,u-1}s_{j,v-1}\Bigg{)}.$
(SI.6.350)

#### Gaussian integral

We want to simplify the expression above by introducing new variables that
will become order parameters. We do so by introducing a double Gaussian
integral with the form:

$\displaystyle\exp\left(Cxy\right)$
$\displaystyle=\exp\left(\frac{C}{2}\left(\frac{1}{2}(x+y)^{2}+\frac{1}{2}(\mathrm{i}\mkern
1.0mu(x-y))^{2}\right)\right)$ (SI.6.351)
$\displaystyle=\frac{C}{\pi}\int\mathrm{d}z_{R}\mathrm{d}z_{I}\exp\left(\frac{C}{2}\left(-\frac{1}{2}z_{R}^{2}-\frac{1}{2}z_{I}^{2}+(x+y)z_{R}+\mathrm{i}\mkern
1.0mu(x-y)z_{I}\right)\right)$ (SI.6.352)
$\displaystyle=\frac{C}{\pi}\int\mathrm{d}z_{R}\mathrm{d}z_{I}\exp\left(\frac{C}{2}\left(-\frac{1}{2}z_{R}^{2}-\frac{1}{2}z_{I}^{2}+x(z_{R}+\mathrm{i}\mkern
1.0muz_{I})+y(z_{R}-\mathrm{i}\mkern 1.0muz_{I})\right)\right)$ (SI.6.353)
$\displaystyle=\frac{C}{\pi}\int\mathrm{d}z_{1}\mathrm{d}z_{2}\exp\left(C\left(-z_{1}z_{2}+xz_{1}+yz_{2}\right)\right),$
(SI.6.354)

where in the last step we applied a change of variables
$z_{1}=\frac{1}{2}(z_{R}+\mathrm{i}\mkern
1.0muz_{I}),z_{2}=\frac{1}{2}(z_{R}-\mathrm{i}\mkern 1.0muz_{I})$.

We apply the Gaussian integral above:

$\displaystyle\exp\left(N\beta
J_{0}\frac{1}{N}\sum_{j}s_{j,u-1}\frac{1}{N}\sum_{i}(-\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u})\right)=$ $\displaystyle\frac{N\beta
J_{0}}{\pi}\int\mathrm{d}\mu_{u}\mathrm{d}m_{u-1}\exp\Bigg{(}-N\beta
J_{0}\mu_{u}m_{u-1}$ (SI.6.355) $\displaystyle+\beta
J_{0}\mu_{u}\sum_{j}s_{j,u-1}-\beta J_{0}m_{u-1}\sum_{i}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\Bigg{)},$ (SI.6.356)

by using $z_{1}=\mu_{u}$, $z_{2}=m_{u-1}$, $x=\frac{1}{N}\sum_{j}s_{j,u-1}$,
$y=\frac{1}{N}\sum_{i}(-\mathrm{i}\mkern 1.0mu\hat{\theta}_{i,u})$, and
$C=N\beta J_{0}$. Similarly, we have

$\displaystyle\exp\left(N\beta^{2}J_{\sigma}^{2}\frac{1}{N}\sum_{j}s_{j,u-1}s_{j,v-1}\frac{1}{N}\sum_{i}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\mathrm{i}\mkern 1.0mu\hat{\theta}_{i,v}\right)=$
$\displaystyle\frac{N\beta^{2}J_{\sigma}^{2}}{\pi}\int\mathrm{d}\rho_{u,v}\mathrm{d}q_{u-1,v-1}\exp\Bigg{(}-N\beta^{2}J_{\sigma}^{2}\rho_{u,v}q_{u-1,v-1}$
(SI.6.357)
$\displaystyle+\beta^{2}J_{\sigma}^{2}\mu_{u}\sum_{j}s_{j,u-1}s_{j,v-1}+\beta^{2}J_{\sigma}^{2}q_{u-1,v-1}\sum_{i}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\mathrm{i}\mkern 1.0mu\hat{\theta}_{i,v}\Bigg{)}.$
(SI.6.358)

With these transformations, we can rewrite

$\displaystyle\left[Z_{t}(\mathbf{g})\right]=$ $\displaystyle\frac{(N\beta
J_{0})^{t}(N\beta^{2}J_{\sigma}^{2})^{t}}{\pi^{2t}}\int\mathrm{d}\mathbf{m}\mathrm{d}\bm{\mu}\mathrm{d}\mathbf{q}\mathrm{d}\bm{\rho}\exp\Bigg{(}-N\beta
J_{0}\sum_{u}\mu_{u}m_{u-1}-N\beta^{2}J_{\sigma}^{2}\sum_{u>v}\rho_{u,v}q_{u-1,v-1}$
(SI.6.359)
$\displaystyle+\log\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}\Bigg{)},$
(SI.6.360)

where the rest of the terms from the Gaussian integral are contained in the
terms

$\displaystyle\Phi(\mathbf{s},\bm{\theta})=$
$\displaystyle\sum_{iu}\left(g_{i,u}+\theta_{i,u}\right)s_{i,u}+\sum_{iu}\beta
J_{0}\mu_{u}s_{i,u-1}+\sum_{i,u>v}\beta^{2}J_{\sigma}^{2}\rho_{u,v}s_{i,u-1}s_{i,v-1}-\sum_{iu}\log
2\cosh\left(\theta_{i,u}\right),$ (SI.6.361)
$\displaystyle\Omega(\bm{\hat{\theta}},\bm{\theta})=$
$\displaystyle\sum_{iu}(\theta_{i,u}-\beta H_{i,u}-\beta
J_{0}m_{u-1})\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}+\frac{\beta^{2}J_{\sigma}^{2}}{2}\sum_{i,u}\left(\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\right)^{2}+\beta^{2}J_{\sigma}^{2}\sum_{i,u>v}q_{u-1,v-1}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\mathrm{i}\mkern 1.0mu\hat{\theta}_{i,v}-Nt\log 2\pi.$
(SI.6.362)

We will see now that the values of $\mathbf{m},\bm{\mu},\mathbf{q},\bm{\rho}$
will act as order parameters of the system.

#### Order parameters

The exponent of the above integrand is proportional to $N$, being it possible
to evaluate the integral by steepest descent, giving the saddle-point solution
as

$\displaystyle\left[Z_{t}(\mathbf{g})\right]=$
$\displaystyle\exp\Bigg{\\{}-N\beta
J_{0}\sum_{u}\mu_{u}m_{u-1}-N\beta^{2}J_{\sigma}^{2}\sum_{u>v}\rho_{u,v}q_{u-1,v-1}+\log\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}\Bigg{\\}},$
(SI.6.363)

where the values of $\mathbf{m},\bm{\mu},\mathbf{q},\bm{\rho}$ are chosen to
extremize (maximize or minimize) the quantity between the braces $\\{\\}$.
Notice that integration over disordered connections has removed coupling
between units and replaced it with same-unit temporal couplings
$\mathbf{\rho}$ and varying effective fields, which are also independent
between units, resulting in a mean-field solution where the activity of
different spins is independent.

From here, knowing that
$\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial\left[Z_{t}(\mathbf{g})\right]}{\partial
g_{i,u}}=\left[\left\langle s_{i,u}\right\rangle\right]$,
$\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial^{2}\left[Z_{t}(\mathbf{g})\right]}{\partial
g_{i,u}\partial g_{j,v}}=\left[\left\langle
s_{i,u}s_{j,v}\right\rangle\right]$, and taking into account that
$\lim_{\mathbf{g}\to\mathbf{0}}\left[Z_{t}(\mathbf{g})\right]=1$ we can
compute the order parameters of the system as

$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial\left[Z_{t}(\mathbf{g})\right]}{\partial
g_{i,u}}=$
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}s_{i,u}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}\left[Z_{t}(\mathbf{g})\right]=\left\langle
s_{i,u}\right\rangle_{*}=\left[\left\langle s_{i,u}\right\rangle\right],$
(SI.6.364)
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial^{2}\left[Z_{t}(\mathbf{g})\right]}{\partial
g_{i,u}\partial g_{j,v}}=$
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}s_{i,u}s_{j,v}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}\left[Z_{t}(\mathbf{g})\right]=\left\langle
s_{i,u}s_{j,v}\right\rangle_{*}=\left[\left\langle
s_{i,u}s_{j,v}\right\rangle\right],$ (SI.6.365)
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial\left[Z_{t}(\mathbf{g})\right]}{\partial
H_{i,u}}=$
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}-\beta\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}\left[Z_{t}(\mathbf{g})\right]=-\beta\left\langle\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\right\rangle_{*}=0,$ (SI.6.366)
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial^{2}\left[Z_{t}(\mathbf{g})\right]}{\partial
H_{i,u}\partial H_{j,v}}=$
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\beta^{2}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\mathrm{i}\mkern
1.0mu\hat{\theta}_{j,v}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}\left[Z_{t}(\mathbf{g})\right]=\beta^{2}\left\langle\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\mathrm{i}\mkern
1.0mu\hat{\theta}_{j,v}\right\rangle_{*}=0,$ (SI.6.367)

where

$\displaystyle\left\langle
f(\mathbf{s},\bm{\hat{\theta}})\right\rangle_{*}=\frac{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}f(\mathbf{s},\bm{\hat{\theta}})\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}.$
(SI.6.368)

Here we should note that, as there is no coupling between units, for $i\neq j$
we have a factorised solution $\left[\left\langle
s_{i,u}s_{j,v}\right\rangle\right]=\left\langle
s_{i,u}s_{j,v}\right\rangle_{*}=\left\langle s_{i,u}\right\rangle\left\langle
s_{j,v}\right\rangle_{*}$.

To obtain the values of the order parameters, we extremize the contents of the
brackets, finding

$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial\log\left[Z_{t}(\mathbf{g})\right]}{\partial\mu_{u+1}}=$
$\displaystyle-N\beta
J_{0}m_{u}+\lim_{\mathbf{g}\to\mathbf{0}}\frac{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\beta
J_{0}\sum_{i}s_{i,u}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}$
(SI.6.369) $\displaystyle=$ $\displaystyle\beta
J_{0}\left(\sum_{i}\left\langle
s_{i,u}\right\rangle_{*}-Nm_{u}\right)=0;\qquad
m_{u}=\frac{1}{N}\sum_{i}\left[\left\langle s_{i,u}\right\rangle\right],$
(SI.6.370)
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial\log\left[Z_{t}(\mathbf{g})\right]}{\partial
m_{u-1}}=$ $\displaystyle-N\beta
J_{0}\mu_{u}+\lim_{\mathbf{g}\to\mathbf{0}}\frac{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}-\beta
J_{0}\sum_{i}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}$
(SI.6.371) $\displaystyle=$ $\displaystyle-\beta
J_{0}\left(\sum_{i}\left\langle\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\right\rangle_{*}+N\mu_{u}\right)=0;\qquad\mu_{u}=0,$
(SI.6.372)
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial\log\left[Z_{t}(\mathbf{g})\right]}{\partial\rho_{u+1,v+1}}=$
$\displaystyle-N\beta^{2}J_{\sigma}^{2}q_{u,v}+\lim_{\mathbf{g}\to\mathbf{0}}\frac{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\beta^{2}J_{\sigma}^{2}\sum_{i}s_{i,u}s_{i,v}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}$
(SI.6.373) $\displaystyle=$
$\displaystyle\beta^{2}J_{\sigma}^{2}\left(\sum_{i}\left\langle
s_{i,u}s_{i,v}\right\rangle_{*}-Nq_{uv}\right)=0;\qquad
q_{uv}=\frac{1}{N}\sum_{i}\left[\left\langle
s_{i,u}s_{i,v}\right\rangle\right],$ (SI.6.374)
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial\log\left[Z_{t}(\mathbf{g})\right]}{\partial
q_{u-1,v-1}}=$
$\displaystyle-N\beta^{2}J_{\sigma}^{2}\rho_{u,v}+\lim_{\mathbf{g}\to\mathbf{0}}\frac{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\beta^{2}J_{\sigma}^{2}\sum_{i}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\mathrm{i}\mkern
1.0mu\hat{\theta}_{j,v}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}{\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}}$
(SI.6.375) $\displaystyle=$
$\displaystyle\beta^{2}J_{\sigma}^{2}\left(\sum_{i}\left\langle\mathrm{i}\mkern
1.0mu\theta_{i,u}\mathrm{i}\mkern
1.0mu\theta_{i,v}\right\rangle_{*}-N\rho_{uv}\right)=0;\qquad\rho_{uv}=0.$
(SI.6.376)

### Saddle-point solution

After solving the saddle-point integral, we have the following dynamical
partition function

$\displaystyle\left[Z_{t}(\mathbf{g})\right]$
$\displaystyle=\sum_{\mathbf{s}_{1:t}}\int\mathrm{d}\bm{\theta}\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})+\Omega(\bm{\hat{\theta}},\bm{\theta})}.$
(SI.6.377)

At this point, we want to remove the effective conjugate fields
$\bm{\hat{\theta}}$ by recovering a delta function. We first rewrite

$\displaystyle\mathrm{e}^{\Omega(\bm{\hat{\theta}},\bm{\theta})}=$
$\displaystyle\frac{1}{(2\pi)^{Nt}}\exp\left(\sum_{iu}(\theta_{i,u}-\beta
H_{i,u}-\beta J_{0}m_{u-1})\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}+\frac{\beta^{2}J_{\sigma}^{2}}{2}\sum_{i,uv}q_{u-1,v-1}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}\mathrm{i}\mkern 1.0mu\hat{\theta}_{i,v}\right),$
(SI.6.378)

where we defined $q_{u-1,u-1}=1$ and $q_{u-1,v-1}=q_{v-1,u-1}$. We can remove
the quadratic terms of $\bm{\hat{\theta}}$ by applying a multivariate Gaussian
integral of the form

$\displaystyle\mathrm{e}^{\frac{1}{2}\sum_{uv}K_{uv}x_{u}x_{v}}=\frac{1}{\sqrt{(2\pi)^{t}|K^{-1}|}}\int\mathrm{d}\mathbf{z}\mathrm{e}^{-\frac{1}{2}\sum_{uv}K_{uv}z_{u}z_{v}+\sum_{uv}K_{uv}x_{u}z_{v}}.$
(SI.6.379)

For $x_{u}=-\beta J_{\sigma}\mathrm{i}\mkern 1.0mu\hat{\theta}_{i,u}$ and
$K_{uv}=q_{u-1,v-1}$, we get

$\displaystyle\int\mathrm{d}\bm{\hat{\theta}}\mathrm{e}^{\Omega(\bm{\hat{\theta}},\bm{\theta})}=$
$\displaystyle\frac{1}{(2\pi)^{Nt}}\int\mathrm{d}\bm{\hat{\theta}}\mathrm{d}\mathbf{z}P(\mathbf{z})\exp\Bigg{(}\sum_{iu}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}(\theta_{i,u}-\beta H_{i,u}-\beta J_{0}m_{u-1})-\beta
J_{\sigma}\sum_{i,uv}q_{u-1,v-1}\mathrm{i}\mkern
1.0mu\hat{\theta}_{i,u}z_{i,v}\Bigg{)}$ (SI.6.380) $\displaystyle=$
$\displaystyle\int\mathrm{d}\mathbf{z}P(\mathbf{z})\prod_{iu}\delta\Bigg{(}\theta_{i,u}-\hat{H}_{i,u}(\mathbf{z})\Bigg{)},$
(SI.6.381)

where $\hat{H}_{i,u}(\mathbf{z})$ is

$\displaystyle\hat{H}_{i,u}(\mathbf{z})=$ $\displaystyle\beta H_{i,u}+\beta
J_{0}m_{u-1}+\beta J_{\sigma}\sum_{v}z_{i,v}q_{u-1,v-1}.$ (SI.6.382)

Here the distribution $P(\mathbf{z})$ is a multivariate Gaussian
$\mathcal{N}(0,\mathbf{\Sigma})$, with
$(\mathbf{\Sigma}^{-1})_{uv}=q_{u-1,v-1}=\frac{1}{N}\sum_{i}\left[\left\langle
s_{i,u-1}s_{i,v-1}\right\rangle\right]$.

To simplify calculations, we can perform a change of variables

$\displaystyle\tilde{z}_{i,u}=\sum_{v}z_{i,v}q_{u-1,v-1}.$ (SI.6.383)

Here, we can find that the covariance matrix of $\bm{\tilde{z}}_{i}$ is the
inverse of the covariance matrix of $\bm{z}_{i}$. We can show this which we
can describe using the characteristic function of the multivariate Gaussian
$P(\mathbf{z})$

$\displaystyle\phi_{\mathbf{z}_{i}}(\bm{\tau})=\left\langle\exp\left(\mathrm{i}\mkern
1.0mu\sum_{u}\tau_{u}z_{i,u}\right)\right\rangle_{\mathbf{z}_{i}}=\exp\left(-\frac{1}{2}\sum_{uv}\tau_{u}\tau_{v}\Sigma_{uv}\right),$
(SI.6.384)

and computing then the characteristic function of $P(\bm{\tilde{z}}_{i})$

$\displaystyle\phi_{\bm{\tilde{z}}_{i}}(\bm{\tau})=$
$\displaystyle\left\langle\exp\left(\mathrm{i}\mkern
1.0mu\sum_{u}\tau_{u}\tilde{z}_{i,u}\right)\right\rangle_{\bm{\tilde{z}}_{i}}=\left\langle\exp\left(\mathrm{i}\mkern
1.0mu\sum_{u}\tau_{u}\sum_{v}z_{i,v}q_{u-1,v-1}\right)\right\rangle_{\mathbf{z}_{i}}$
(SI.6.385) $\displaystyle=$
$\displaystyle\exp\left(-\frac{1}{2}\sum_{uu^{\prime}}\tau_{u}\tau_{u^{\prime}}\sum_{vv^{\prime}}q_{u-1,v-1}\Sigma_{vv^{\prime}}q_{u^{\prime}-1,v^{\prime}-1}\right)$
(SI.6.386) $\displaystyle=$
$\displaystyle\exp\left(-\frac{1}{2}\sum_{uv}\tau_{u}\tau_{v}q_{u-1,v-1}\right),$
(SI.6.387)

where in the last step we took into account that
$\bm{\tau}^{T}\mathbf{q}^{T}\mathbf{\Sigma}\mathbf{q}\bm{\tau}=\bm{\tau}^{T}\mathbf{q}\bm{\tau}$
as $\mathbf{\Sigma}^{-1}=\mathbf{q}$. This result indicates that
$P(\bm{\tilde{z}}_{i})$ is a multivariate Gaussian
$\mathcal{N}(0,\mathbf{q})$.

After solving the saddle-point equation we have

$\displaystyle\Phi(\mathbf{s},\bm{\theta})=$
$\displaystyle\sum_{iu}s_{i,u}\left(g_{i,u}+\theta_{i,u}\right)-\sum_{iu}\log
2\cosh\left(\theta_{i,u}\right),$ (SI.6.388)

which leads us to

$\displaystyle\left[Z_{t}(\mathbf{g})\right]=$
$\displaystyle\int\mathrm{d}\bm{\theta}\exp\left(\log\sum_{\mathbf{s}_{1:t}}\mathrm{e}^{\Phi(\mathbf{s},\bm{\theta})}\right)\int\mathrm{d}\mathbf{\tilde{z}}P(\mathbf{\tilde{z}})\prod_{iu}\delta\Bigg{(}\theta_{i,u}-\hat{H}_{i,u}(\mathbf{\tilde{z}})\Bigg{)}$
(SI.6.389) $\displaystyle=$
$\displaystyle\int\mathrm{d}\mathbf{\tilde{z}}P(\mathbf{\tilde{z}})\exp\Bigg{(}\log\Bigg{(}\sum_{\mathbf{s}_{1:t}}\exp\left(\sum_{iu}s_{i,u}\left(g_{i,u}+\hat{H}_{i,u}(\mathbf{\tilde{z}})\right)\right)-\sum_{iu}\log
2\cosh\left(\hat{H}_{i,u}(\mathbf{\tilde{z}})\right)\Bigg{)}$ (SI.6.390)
$\displaystyle=$
$\displaystyle\int\mathrm{d}\mathbf{\tilde{z}}P(\mathbf{\tilde{z}})\prod_{iu}\frac{\cosh\left(g_{i,u}+\hat{H}_{i,u}(\mathbf{\tilde{z}})\right)}{\cosh\left(\hat{H}_{i,u}(\mathbf{\tilde{z}})\right)}.$
(SI.6.391)

And, as the diagonal of the covariance matrix of $\mathbf{\tilde{z}}$ is equal
to 1, we can derive

$\displaystyle m_{i,u}=$
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial
Z_{t}(\mathbf{g})}{\partial g_{i,u}}=\int\mathrm{D}z\tanh(\hat{H}_{i,u}(z)),$
(SI.6.392) $\displaystyle R_{ii,uv}=$
$\displaystyle\lim_{\mathbf{g}\to\mathbf{0}}\frac{\partial^{2}Z_{t}(\mathbf{g})}{\partial
g_{i,u}\partial
g_{i,v}}=\int\mathrm{D}xy(q_{u-1,v-1})\tanh(\hat{H}_{i,u}(x))\tanh(\hat{H}_{i,v}(y)),$
(SI.6.393)

where

$\displaystyle\mathrm{D}z=$
$\displaystyle\frac{1}{\sqrt{2\pi}}\mathrm{e}^{-\frac{1}{2}z^{2}},$ (SI.6.394)
$\displaystyle\mathrm{D}xy(q)=$
$\displaystyle\frac{1}{2\pi\sqrt{1-q^{2}}}\mathrm{e}^{-\frac{x^{2}+y^{2}-2qxy}{2(1-q^{2})}},$
(SI.6.395) $\displaystyle\hat{H}_{i,u}(z)=$ $\displaystyle\beta H_{i,u}+\beta
J_{0}m_{u-1}+\beta J_{\sigma}z.$ (SI.6.396)

The resulting equations are similar to the symmetric SK model, but the same-
spin coupling parameter vanishes for the computation of $m_{u}$.

$\displaystyle m_{u}=$
$\displaystyle\frac{1}{N}\sum_{i}\int\mathrm{D}z\tanh\left(\beta H_{i,u}+\beta
J_{0}m_{u-1}+\beta J_{\sigma}z\right),$ (SI.6.397) $\displaystyle q_{uv}=$
$\displaystyle\frac{1}{N}\sum_{i}\int\mathrm{D}xy(q_{u-1,v-1})\tanh\left(\beta
H_{i,u}+\beta J_{0}m_{u-1}+\beta J_{\sigma}x\right)\tanh\left(H_{i,v}+\beta
J_{0}m_{v-1}+\beta J_{\sigma}y\right).$ (SI.6.398)

This is consistent with findings of the asymmetric SK model lacking a spin-
glass phase as $m_{u}$ is independent of $q_{u,v}$. Note also that all
$q_{u,v}$ only depends on the previous $q_{u-1,v-1}$, this in the $t\to\infty$
limit the value of $q_{t,t-d}$ tends to the same value for any finite $d$.

### Ferromagnetic critical phase transition in the infinite kinetic Ising
model with Gaussian couplings and uniform weights

We define an Ising network of infinite size, with randomly defined bias
$H_{i,u}=H_{i}$, where $H_{i}$ has a distribution $\mathcal{U}(-H_{0},H_{0})$
and couplings $J_{ij}$ with a Gaussian distribution
$\mathcal{N}(\frac{1}{N},\frac{J_{\sigma}^{2}}{N})$. We choose $H_{0}=0.5$,
$J_{\sigma}=0.1$.

As we have found that the asymmetric SK model with arbitrary fields follows a
mean-field solution, calculating the effects of disorder in the fields becomes
much easier, as we can approximate the update equations of the order
parameters in the thermodynamic limit $N\to\infty$ as:

$\displaystyle m_{t}=\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}m_{i,t}=$
$\displaystyle\frac{1}{2H_{0}}\int_{-H_{0}}^{H_{0}}\mathrm{d}h\int\mathrm{D}z\tanh\left(\beta\left(h+J_{0}m_{t-1}+J_{\sigma}z\right)\right)$
(SI.6.399) $\displaystyle=$ $\displaystyle\frac{1}{2\beta
H_{0}}\int\mathrm{D}z\log\left(\frac{\cosh\left(\beta\left(H_{0}+J_{0}m_{t-1}+J_{\sigma}z\right)\right)}{\cosh\left(\beta\left(-H_{0}+J_{0}m_{t-1}+J_{\sigma}z\right)\right)}\right).$
(SI.6.400)

Similarly, the delayed self-correlation parameter is given as

$\displaystyle
q_{t,t^{\prime}}=\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}R_{ii,t,t^{\prime}}=$
$\displaystyle\frac{1}{2H_{0}}\int_{-H_{0}}^{H_{0}}\mathrm{d}h\int\mathrm{D}xy(q_{t-1,t^{\prime}-1})\tanh\left(\beta\left(h+J_{0}m_{t-1}+J_{\sigma}x\right)\right)\tanh\left(\beta\left(h+J_{0}m_{t^{\prime}-1}+J_{\sigma}y\right)\right)$
(SI.6.401) $\displaystyle=$ $\displaystyle\frac{1}{2\beta
H_{0}}\int\mathrm{D}xy(q_{t-1,t^{\prime}-1})\Bigg{(}\frac{1}{2}\left(\log\left(\frac{1+\tanh\left(\beta\left(H_{0}+H_{t,x}\right)\right)}{1-\tanh\left(\beta\left(H_{0}+H_{t,x}\right)\right)}\frac{1-\tanh\left(\beta\left(-H_{0}+H_{t,x}\right)\right)}{1+\tanh\left(\beta\left(-H_{0}+H_{t,x}\right)\right)}\right)\right)$
(SI.6.402)
$\displaystyle-\frac{1}{\tanh\left(\beta\left(H_{t^{\prime},y}-H_{t,x}\right)\right)}\log\left(\frac{1+\tanh\left(\beta\left(H_{t^{\prime},y}-H_{t,x}\right)\right)\tanh\left(\beta\left(H_{0}+H_{t,x}\right)\right)}{1+\tanh\left(\beta\left(H_{t^{\prime},y}-H_{t,x}\right)\right)\tanh\left(\beta\left(-H_{0}+H_{t,x}\right)\right)}\right)\Bigg{)}.$
(SI.6.403)

with $H_{t,x}=J_{0}m_{t-1}+J_{\sigma}x$ and
$H_{t^{\prime},y}=J_{0}m_{t^{\prime}-1}+J_{\sigma}y$.

By recursively updating Supplementary Eqs. SI.6.400 and SI.6.402-SI.6.403, we
can obtain the order parameters of the system. In Supplementary Fig. 6, we
show the order parameters for $J_{0}=1$, $J_{\sigma}=0.1$ and $H_{0}=0.5$,
which resulted in a ferromagnetic transition around $\beta_{c}\approx 1.1108$.

A | B
---|---
|
C | D
|
Figure 6: Analytical results of the behaviour of the infinite kinetic Ising
model. A critical point is found for approximately $\beta_{c}=1.1108$ with
parameters $J_{0}=1$, $J_{\sigma}=0.1$ and $H_{0}=0.5$.

Similarly, assuming a nonequilibrium steady state (NESS) in which
$m_{t}=m_{t-1}=m$ we can obtain the critical point of the system by computing
the non-zero solutions of the first order Taylor expansion around $m=0$ of the
right hand part of Supplementary Eq. SI.6.400,

$\displaystyle m\approx$ $\displaystyle\frac{1}{2\beta
H_{0}}\int\mathrm{D}z\log\left(\frac{\cosh\left(\beta\left(H_{0}+J_{\sigma}z\right)\right)}{\cosh\left(\beta\left(-H_{0}+J_{\sigma}z\right)\right)}\right)$
(SI.6.404) $\displaystyle+$ $\displaystyle\frac{1}{2\beta
H_{0}}\int\mathrm{D}z\left(\tanh\left(\beta\left(H_{0}+J_{\sigma}z\right)\right)-\tanh\left(\beta\left(-H_{0}+J_{\sigma}z\right)\right)\right)\beta
J_{0}m$ (SI.6.405) $\displaystyle=$
$\displaystyle\frac{1}{H_{0}}\int\mathrm{D}z\tanh\left(\beta\left(H_{0}+J_{\sigma}z\right)\right)J_{0}m,$
(SI.6.406)

that yields the same solution for $\beta_{c}$ solving the equation

$\displaystyle\frac{1}{H_{0}}\int\mathrm{D}z\tanh\left(\beta_{c}\left(H_{0}+J_{\sigma}z\right)\right)J_{0}=1.$
(SI.6.407)

We can characterize some critical exponents of the system using the reduced
temperature $\tau=-\frac{\beta-\beta_{c}}{\beta_{c}}$ , and finding that near
the critical point the first order Taylor expansion around $\beta_{c}$

$\displaystyle\frac{1}{H_{0}}\int\mathrm{D}z\tanh\left(\beta\left(H_{0}+J_{\sigma}z\right)\right)J_{0}\approx$
$\displaystyle
1+\frac{1}{H_{0}}\int\mathrm{D}z\left(1-\tanh^{2}\left(\beta\left(H_{0}+J_{\sigma}z\right)\right)\right)\left(H_{0}+J_{\sigma}z\right)J_{0}\left(\beta-\beta_{c}\right)$
(SI.6.408) $\displaystyle=$ $\displaystyle 1-K^{\prime}\tau$ (SI.6.409)

and computing the value of $m$ around $\beta=\beta_{c}$ with the third order
Taylor expansion around $m=0$

$\displaystyle m\approx$
$\displaystyle\frac{1}{H_{0}}\int\mathrm{D}z\tanh\left(\beta\left(H_{0}+J_{\sigma}z\right)\right)J_{0}m$
(SI.6.410) $\displaystyle-$ $\displaystyle\frac{1}{3\beta
H_{0}}\int\mathrm{D}z\tanh\left(\beta\left(H_{0}+J_{\sigma}z\right)\right)(1-\tanh^{2}\left(\beta\left(H_{0}+J_{\sigma}z\right)\right))\left(\beta
J_{0}m\right)^{3}$ (SI.6.411) $\displaystyle=$
$\displaystyle(1-K^{\prime}\tau)m-K^{\prime\prime}m^{3},$ (SI.6.412)
$\displaystyle m\propto$ $\displaystyle\left(-\tau\right)^{-\frac{1}{2}}.$
(SI.6.413)

This yields a critical exponent $\beta^{\prime}=\frac{1}{2}$ (note that this
is different from the inverse temperature $\beta$), consistent with the mean-
field universality class.

Similarly, we can compute the susceptibility to a uniform magnetic field $B$
added on top of $H_{i}$, having that the derivative of $m$ with respect to $B$
evaluated at $B=0$

$\displaystyle\frac{\partial m}{\partial B}=\frac{1}{2\beta
H_{0}}\int\mathrm{D}z\left(\tanh\left(\beta\left(H_{0}+J_{0}m+J_{\sigma}z\right)\right)-\tanh\left(\beta\left(-H_{0}+J_{0}m+J_{\sigma}z\right)\right)\right)\left(\beta+\beta
J_{0}\frac{\partial m}{\partial B}\right).$ (SI.6.414)

The susceptibility evaluated at $m=0$ yields, for the limit $\tau\to 0$,

$\displaystyle\frac{\partial m}{\partial B}=$
$\displaystyle\left(1-K^{\prime}\tau\right)\left(\frac{1}{J_{0}}+\frac{\partial
m}{\partial B}\right),$ (SI.6.415) $\displaystyle\frac{\partial m}{\partial
B}\propto$ $\displaystyle\frac{1-K^{\prime}\tau}{\tau}\approx\tau^{-1}.$
(SI.6.416)

The result retrieves the $\gamma=1$ exponent consistent with the mean-field
universality class.

#### Finite-size fluctuations

Finally, to evaluate the behaviour around $\beta_{c}$ in the finite network
used in the article, we numerically simulated the behaviour of a network of
size $N=512$, with parameters $H_{0}=0.5$, $J_{0}=1,J_{\sigma}=0.1$ as in the
numerical results from the main text. In Supplementary Fig. 7 we show the
values of $\mathbf{m}_{t}$. $\mathbf{C}_{t}$ and $\mathbf{D}_{t}$, showing
that correlations peak around the value of $\beta_{c}$ computed theoretically.

A | B | C
---|---|---
|  |
Figure 7: Finite-size simulation of the asymmetric Ising model. Experimental
results of the behaviour around the critical point ($\beta_{c}=1.1108$) of an
asymmetric Ising model with $N=512$ and parameters $H_{0}=0.5$, $J_{0}=1$ and
$J_{\sigma}=0.1$.

##  Supplementary Note 7 Complexity of the models and their computational
costs

Here we clarify the structure of each model, and quantitatively compare the
relation between model complexity and their computational cost. We start by
noting that the reference model in Plefka[$t-1,t$] is a submodel of the
reference models in the other three approximation models used in Plefka[$t$],
Plefka[$t-1$], and Plefka2[$t$]. The latter approximation models (Plefka[$t$],
Plefka[$t-1$], Plefka2[$t$]) are not mutually inclusive, rather each model
manifests a distinct assumption.

To clarify the relation between the performance of each approximation and its
model complexity, we provide the number of free and fixed parameters of each
model. We then examine the relation between accuracy of the models in
forward/inverse Ising problems and their computational time required to solve
these problems.

### Number of free and fixed parameters in the models

To compare the structure of each approximated model, let us first revisit the
definition of each model, and introduce the number of their free or fixed
parameters. Each approximation of $P$ (a marginal distribution of spins
covering some time steps $t$ of the model) defines a reference manifold using
$Q$ to perform a Plefka expansion. $Q$ is defined at each time point, and
obtained by fixing some parameters (couplings) of $P$ at zero. The rest of the
parameters are either free parameters that are fitted as an m-projection from
$P$, or preserved as the original values of $P$. In the Supplementary Table 1,
we listed the number of free and fixed (either to zero or to their original
value) parameters of reference models. For example, in Plefka[$t$], we
preserve one field parameter ($H_{i}$) and $N$ couplings ($J_{ij}$) for each
$i$-th spin at time $t-1$ with their original values, while fixing $N$
couplings of the $i$-th spin at time $t$ at zero. A field parameter of the
$i$-th neuron at time $t$ ($\Theta_{i,t}$) is a free parameter obtained as an
m-projection from $P$ (i.e., it is fitted so that $Q$ and $P$ have the same
expectation $\mathbf{m}$). We need to perform this projection for all neurons
(i.e., $N$ times).

|  Number of times the approximation is computed |  Free parameters to be fitted (at each approximation) |  Parameters fixed at zero (at each approximation) |  Parameters fixed at their original value (for each approximation)
---|---|---|---|---
Plefka[$t-1,t$] |  $2N$ ($N$ neurons at time $t$ and $N$ neurons at time $t-1$) |  $1$ field |  $N$ couplings |  $0$
Plefka[$t$] |  $N$ ($N$ neurons at time $t$) |  $1$ field |  $N$ couplings |  $1$ field + $N$ couplings (for $N$ neurons at $t-1$)
Plefka[$t-1$] |  $N$ ($N$ neurons at time $t-1$) |  $1$ field |  $N$ couplings |  $1$ field + $N$ couplings (for $N$ neurons at $t$)
Plefka2[$t$] |  $2N^{2}$ (once for every pair of neurons at time $t$ once for every pair at $t,t-1$) |  $2$ fields (time $t$ and $t-1$) and $1$ coupling |  $N-1$ couplings at $t$, $N$ couplings at $t-1$ |  $1$ field and $N$ couplings (for $N-1$ neurons except for $l$ at $t-1$)
Table 1: Comparison of methods. The number of free and fixed parameters of
each approximation model.

### Accuracy versus computational costs

The structure of the reference model $Q$ influences both the accuracy of the
approximation methods and their computational complexity. For example,
Plefka[$t$] is computationally more demanding than Plefka[$t-1,t$] since the
former requires to compute multiplications over the correlation matrices
$\mathbf{C}_{t-1}$, while the latter only makes use of the means
$\mathbf{m}_{t-1}$. Plefka[$t-1$] involves performing a Gaussian integral,
which becomes computationally demanding, especially for the equal-time
correlations, where one needs to perform a 2-dimensional integral.
Plefka2[$t$] is computationally more expensive than Plefka[$t$] because, while
the latter needs to be computed $N$ times (once per spin), the former needs to
be computed $2N^{2}$ times (once per each pair of delayed and same-time
spins).

In Supplementary Fig. 8, we show the trade-off between computational time and
accuracy for the different methods. These results are obtained by simulating
the model for $T=128$ steps at $\beta=\beta_{c}$ (in principle the most
challenging point of our model). We computed the average squared error of its
statistical moments and cumulants
$\epsilon_{\mathbf{m}}=\langle\langle(m^{o}_{i,t}-m^{p}_{i,t})^{2}\rangle_{i}\rangle_{t}$,
$\epsilon_{\mathbf{C}}=\langle\langle(C^{o}_{ik,t}-C^{p}_{ik,t})^{2}\rangle_{ik}\rangle_{t}$
and
$\epsilon_{\mathbf{D}}=\langle\langle(D^{o}_{ik,t}-D^{p}_{ik,t})^{2}\rangle_{il}\rangle_{t}$.

As expected, Plefka[$t-1,t$] is the fastest method, although it yields the
most inaccurate approximation. Plefka[$t$] significantly improves accuracy
with a slight increase in computational time. Plefka2[$t$] increases further
both computational cost and accuracy. Finally, Plefka[$t-1$] is relatively
slow compared to the other methods, offering a poor performance in the forward
Ising problem. In sum, we recommend to use Plefka2[$t$] if users pursue the
accuracy, and Plefka[$t$] if users pursue the speed in the forward Ising
problem.

Figure 8: Accuracy vs computational time in the forward Ising problem. Average
squared error in (A) the Ising model average magnetizations, (B) equal-time
correlations and (C) delayed correlations at $t=1,\dots,128$ versus
computation time in minutes for different mean-field models for
$\beta=\beta_{c}$. Results were obtained using a 2-core Intel(R) Core(TM)
i7-5500U CPU @ 2.40GHz processor.

In Supplementary Fig. 9, we show the results for the inverse Ising problem at
$\beta=\beta_{c}$. The (in)accuracy is evaluated by the average squared error
between the inferred parameters and true parameters:
$\epsilon_{\mathbf{H}}=\langle\langle(H^{o}_{i,t}-H^{p}_{i,t})^{2}\rangle_{i}\rangle_{t}$
and
$\epsilon_{\mathbf{J}}=\langle\langle(J^{o}_{ik,t}-J^{p}_{ik,t})^{2}\rangle_{il}\rangle_{t}$.
We found that the results are different from the forward Ising model. In the
inference problem, Plefka[$t-1,t$] shows a poor performance in accuracy as in
the forward problem, but it was no longer the fastest method. Plefka[$t$],
Plefka[$t-1$] and Plefka2[$t$] offer similar performances.

There are several reasons why the results in the inverse Ising problem differ
from the previous forward problem. First, the computational time in the
inverse problem depends on how fast the learning dynamics reach convergence.
These dynamics differ across the methods (e.g., due to the slowest speed of
convergence, Plefka[$t-1,t$] is slightly slower than Plefka[$t$]). In
addition, we need to compute only $\mathbf{m}_{t}$ and $\mathbf{D}_{t}$, and
do not need to compute $\mathbf{C}_{t}$ in the inverse problem. Since
Plefka[$t-1$] is very slow in computing $\mathbf{C}_{t}$, by freeing from its
computational load, Plefka[$t-1$] becomes much faster in the inverse problem.

Figure 9: Accuracy vs computational time in the inverse Ising problem. Average
squared error in the inverse Ising inferred (A) external fields and (B)
couplings versus computation time in minutes for different mean-field models
for $\beta=\beta_{c}$. Results were obtained using a 2-core Intel(R) Core(TM)
i7-5500U CPU @ 2.40GHz processor.