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⁴⁴). The system consists of N interacting binary variables (down or up of Ising spins or inactive or active of neural units) s_i,t ∈ { − 1, + 1}, i = 1, 2, . . , N, evolving in discrete-time steps t with parallel dynamics. Given the configuration of spins at t − 1, s_t−1 = {s_1,t−1, s_2,t−1, …, s_N,t−1}, spins s_t at time t are conditionally independent random variables, updated as a discrete-time Markov chain, following

The parameters H = {H_i} and 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 ≠ 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(s_t−1), the distribution of the current state is:

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 m_t, correlations between pairs of units (covariance function) C_t, and delayed correlations D_t given by

Note that m_t and 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 C_t as they are commonly used to describe neural systems, and are investigated by some of the mean-field approximations in the literature⁴⁰. 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^15,47,48. 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 submanifolds $$ {\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(s_t∣H, J) in the manifold $$ {\mathcal{P}}_t$$ and its statistical properties m_t, C_t, D_t.

Fig. 1

A geometric view of the approximations based on Plefka expansions.

The point P(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(s_t) and Q^*(s_t) on $$ {\mathcal{Q}}_t$$, where for independent $$ {\mathcal{Q}}_t$$ all points on $$ \mathcal{A}$$ share the same mean values m_t. Geometrically, $$ \mathcal{A}$$ constitutes the m-projection from P(s_t) to $$ {\mathcal{Q}}_t$$, defining Q^*(s_t) as the closest point in the submanifold $$ {\mathcal{Q}}_t$$ to the point P(s_t)⁴⁷. The Plefka expansion is defined by expanding an α-dependent distribution P_α(s_t) that satisfies P_α=0(s_t) = Q^*(s_t) and P_α=1(s_t) = P(s_t).

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

Our first goal is to find the average activation rates of the target distribution P(s_t∣H, J). It turns out that they can be obtained from the independent model Q(s_t∣Θ_t) that minimizes the following Kullback-Leibler (KL) divergence from P(s_t):

The independent model $$ Q\left({\mathbf{s}}_t\mid {\mathbf{\Theta }}_t^{*}\right)\left(\equiv Q^{*}\left({\mathbf{s}}_t\right)\right)$$ that minimizes the KL divergence has activation rates m_t identical to those of the target distribution P(s_t)³⁸ because the minimizing points $$ {\mathrm{\Theta }}_{i,t}^{*}$$ satisfy (for i = 1, …, N)

From an information geometric point of view, given m_t (or $$ {\mathbf{\Theta }}_t^{*}$$), we may consider a family of points defined as a linear mixture of P(s_t) and $$ Q\left({\mathbf{s}}_t\mid {\mathbf{\Theta }}_t^{*}\right)$$ for which 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^37,47. Thus, the previous search of $$ Q\left({\mathbf{s}}_t\mid {\mathbf{\Theta }}_t^{*}\right)$$ given P(s_t∣H, J) is equivalent to finding the orthogonal projection point from P(s_t∣H, J) to the manifold $$ {\mathcal{Q}}_t$$ of independent models^36,37.

The Plefka expansion

Although the m-projection provides the exact and unique average activation rates, its calculation in practice requires the complete distribution P(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(s_t∣H, J) and $$ Q\left({\mathbf{s}}_t\mid {\mathbf{\Theta }}_t^{*}\right)$$ (the solid line in Fig. 1). This distribution is defined using a new conditional distribution introducing a parameter α that connects a distribution on the manifold $$ {\mathcal{Q}}_t$$ with the original distribution P(s_t):

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

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))^36,49. 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(s_t∣H, J) to the submanifold $$ {\mathcal{Q}}_t$$, which does not result in $$ Q\left({\mathbf{s}}_t\mid {\mathbf{\Theta }}_t^{*}\right)$$. Namely, minimizing D(Q∣∣P), as well as minimization of other α-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 m_t using Eq. (12) up to an arbitrary order.

In the subsequent sections we show that different approximations of the marginal distribution P(s_t) in Eq. (3) can be constructed by replacing P(s_i,τ∣s_τ−1) with P_α(s_i,τ∣s_τ−1) for different pairs i, τ (here we will explore the cases of τ = t and τ = 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(s_t−k+1, …, 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(s_t−k+1, …, 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 m_t, C_t, and 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). Here we should note that assuming an approximation where previous distributions (e.g., t − 2, t − 3, … ) 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 }^{\left[t-1:t\right]}$$. Using Eqs. (3) and (10), we write

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

We apply the same expansion to approximate the correlations, expanding C_ik,t(α) and D_il,t(α) around α = 0 up to the first order using $$ {\mathrm{\Theta }}_{i,t}={\mathrm{\Theta }}_{i,t}^{*}$$. Then we obtain

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

Having $$ {\mathrm{\Theta }}_{i,t}^{*}$$, we can incorporate TAP approximations of the correlations by expanding C_ik,t(α) and D_il,t(α) (see Supplementary Note 2 for details) as:

In these approximations, Eqs. (16) and (20) of activation rates 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), (21), and (22)) are novel contributions from applying the Plefka expansion to correlations.

Using the equations above, we can compute the approximate statistical properties of the system at t (m_t, C_t, D_t) from m_t−1. Therefore, the system evolution is described by recursively computing m_t from an initial state m₀ (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 expansion that uses a model containing independent units at time t − 1 and t to construct a marginal probability distribution $$ P_{\alpha }^{\left[t-1:t\right]}\left({\mathbf{s}}_t\right)$$. This is, however, not the only possible choice of approximation. As we mentioned above, other approximations have been introduced in the literature. In ref. ⁴⁰, expressions are provided for the nonstationary delayed correlations D_t as a function of C_t−1. In ref. ⁴¹, an approximation is derived by assuming that units at state s_t−1 are independent while correlations of s_t are preserved.

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 submanifold $$ {\mathcal{Q}}_t$$ at each time step. Fig. 2 shows the corresponding submanifolds $$ {\mathcal{Q}}_{t-1:t}$$ of possible approximations, where gray lines represent interactions that are affected by α 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 α. 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.

Fig. 2

Unified mean-field framework.

Original model (A) and family of generalized Plefka expansions (B–E). Gray lines represent connections that are proportional to α 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 naive and TAP 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 ref. ⁴⁰. Plefka[t − 1] (D) assumes independent activity at t-1, and in its first order approximation reproduces the results in ref. ⁴¹. Plefka2[t] (E) represents a novel pairwise approximation which performs better in approximating correlations.

Plefka[t]: expansion around an independent model at time t

For the Plefka[t − 1, t] approximation, explained above, the system becomes independent for α = 0 at t as well as t − 1. This leads to approximations of m_t, C_t, D_t being specified by m_t−1, while being independent of C_t−1 and D_t−1. In ref. ⁴⁰, the authors describe a mean-field approximation by performing new expansion over the classical nMF and TAP equations that takes into account previous correlations 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(s_t−1) (Fig. 2C). Therefore, we denote this approximation as $$ P_{\alpha }^{\left[t\right]}$$ and consider

Note that Eqs. (24) and (25) are the same as the nMF Plefka[t − 1, t] equations. Equation (26) includes C_t−1, being exactly the same result obtained in ref. ⁴⁰, Eq. (4). The second-order approximations leads to:

All update rules include the effect of 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 the Plefka[t − 1, t] approximation in the previous section. In contrast with ref. ⁴⁰, 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 ref. ⁴¹, a mean-field method is proposed by approximating the effective field 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(s_t−1) to the submanifold $$ {\mathcal{Q}}_{t-1}$$, using a model where only s_t−1 are independent (Fig. 2D). In this case (see Supplementary Note 4), the effective field h_t at the submanifold is a sum of independent terms, which for large N yields a Gaussian distribution.

By defining

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:

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:

These results are related to previous work⁴³ 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 C_t−1 and 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 s_t−1, and then s_i,t is computed conditioned on s_k,t, s_t−1 (Fig. 2E, right):

Here θ_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):

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

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}\left(-\beta H_0,\beta H_0\right)$$, H₀ = 0.5, whereas coupling terms J_ij are sampled from independent Gaussian distributions $$ \mathcal{N}\left(\beta \frac{J_0}{N},{\beta }^2\frac{J_{\sigma }^2}{N}\right)$$, J₀ = 1, J_σ = 0.1, where β 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^5,6. 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 β_c ≈ 1.1108 in thermodynamic limit (see Supplementary Note 6, Supplementary Fig. 1). The uniformly distributed bias terms H shift the phase transition point from β = 1 obtained at H = 0. By simulation of the finite size systems, we confirmed that the maximum fluctuations in the model are found near the theoretical β_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¹⁹. 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⁵⁰ 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 β_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 H and J. To study the nonstationary transient dynamics of the model, we start from s₀ = 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⁶ trials for 21 values of β in the range [0.7β_c, 1.3β_c] (except for the reconstruction of the phase transition where we used R = 10⁵ and 201 values of β in the same range). Using the R samples, we computed the statistical moments and cumulants of the system, m_t, C_t, and D_t at each time step. We then computed their averages over the system units, i.e., $$ {⟨m_{i,t}⟩}_i$$, $$ {⟨C_{ik,t}⟩}_{ik}$$ and $$ {⟨D_{il,t}⟩}_{il}$$, where the angle bracket denotes average over indices of its subscript.

The black solid lines in Fig. 3A–C display nonstationary dynamics of these averaged statistics from t = 0, …, 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 m₀ = 1, C₀ = 0 and D₀ = 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 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 C_t, 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.

Fig. 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 β = β_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 β = β_c and t = 128. Black lines represent the identity line. Bottom: Average squared error of the magnetizations $$ {ϵ}_{\mathbf{m}}={⟨{⟨{\left(m_{i,t}^o-m_{i,t}^p\right)}^2⟩}_i⟩}_t$$ (G), equal-time correlations $$ {ϵ}_{\mathbf{C}}={⟨{⟨{\left(C_{ik,t}^o-C_{ik,t}^p\right)}^2⟩}_{ik}⟩}_t$$ (H), and delayed correlations $$ {ϵ}_{\mathbf{D}}={⟨{⟨{\left(D_{ik,t}^o-D_{ik,t}^p\right)}^2⟩}_{il}⟩}_t$$ (I) for 21 values of β in the range [0.7β_c, 1.3β_c].

Performance of the methods in predicting individual activation rates and correlations are displayed in Fig. 3D–F by comparing vectors m_t, C_t and 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 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 C_t, D_t are best predicted by Plefka2[t] (Fig. 3E, F).

The above results are obtained at the critical β = β_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 β (the same random parameters are applied for all β). Fig. 3G, H, I, respectively, show the averaged squared errors (averaged over time and units) of the activation rates ϵ_m, equal-time correlations ϵ_C and delayed correlations ϵ_D between the original model and approximations, averaged over units and time for 21 values of β in the range [0.7β_c, 1.3β_c]. Fig. 3G–I shows that Plefka2[t] outperforms the other methods in computing m_t, C_t, D_t (with the exception of a certain region of β > β_c in which Plefka[t] is slightly better), yielding consistently a low error bound for all values of β. Errors of these approximations are smaller when the system is away from β_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⁶ trials to infer the parameters of the model, H and J. The model parameters are estimated by the Boltzmann learning method under the maximum likelihood principle: H and J are updated to minimize the differences between the average rates m_t or delayed correlations 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 β = β_c (Fig. 4A, B), we observe that the classical Plefka[t − 1, t] approximation adds significant offset values to the fields H and couplings J. In contrast, Plefka[t], Plefka[t − 1] and Plefka2[t] are all precise in estimating the values of H and J.

Fig. 4

Inverse Ising problem.

Top: Inferred external fields (A) and couplings (B) found by different mean-field models plotted versus the real ones for β = β_c. Black lines represent the identity line. Bottom: Average squared error of inferred external fields $$ {ϵ}_{\mathbf{H}}={⟨{\left(H_i^o-H_i^p\right)}^2⟩}_i$$ (C) and couplings $$ {ϵ}_{\mathbf{J}}={⟨{\left(J_{ij}^o-J_{ij}^p\right)}^2⟩}_{ij}$$ (D) for 21 values of β in the range [0.7β_c, 1.3β_c].

Fig. 4C, D shows the mean squared error ϵ_H, ϵ_J for bias terms and couplings between the original model and the inferred values for different β. In this case, errors are large in the estimation of 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 (β ≈ β_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. 5A–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 β = β_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 C_t and 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 C_t and 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 β_c and Plefka2[t] is more precise than the other methods.

Fig. 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 β = β_c. Bottom (D–F): The same as above using the reconstructed network H, J by solving the inverse Ising problem at β = β_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.

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 take the values of H, J from solving the inverse problem from the data sampled at β = β_c, and next we solve again the forward problem with those estimated parameters rescaled by a new inverse temperature $$ \tilde{\beta }$$. The results for the correlations (Fig. 5D, 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. 5F), 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.

Boltzmann learning in the inverse Ising problem

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

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^17,55,57: