Introduction

While isolated systems tend toward thermodynamic equilibrium, many physical, chemical, and biological processes operate far from equilibrium. Such nonequilibrium systems – from molecules to organisms and machines – persist by exchanging matter and energy with their surroundings, being causally driven by time-varying external stimuli or by their past states (e.g., the adaptive action of sensor and effector interfaces1). Nonequilibrium processes inherently break time reversal symmetry, describing spatial and temporal patterns with a definite past-future order, and being thus strikingly different from the reversible dynamics found at thermodynamic equilibrium. Understanding these dissipative processes – from chemical reactions to neural dynamics or flocks of birds – brings critical insights into the self-organization of open systems2. Although these ideas have attracted the interest of disparate fields, from evolutionary dynamics3 to neuroscience4,5,6,7, little is known about the thermodynamic description of nonequilibrium systems comprising many interacting particles. While stochastic thermodynamics has been greatly influential in the study of small systems with appreciable fluctuations8, the thermodynamics of large-scale nonequilibrium systems and their phase transitions has attracted attention only very recently9,10,11.

When the elements of a system are numerous, characterizing its nonequilibrium states is challenging due to the expansion of its state space. Inspired by the success of the equilibrium Ising model in investigating disordered systems in the thermodynamic limit, we study the nonequilibrium thermodynamics of a stochastic, kinetic Ising model with both synchronous and asynchronous updates. The Ising model is a cornerstone of statistical mechanics, originally conceived as a model describing phase transitions in magnetic materials12. A natural extension of the model introducing Markovian dynamics either in discrete or continuous time is the kinetic Ising model, a prototypical model of both equilibrium and nonequilibrium systems such as recurrent neural networks13 or genetic regulatory networks14. With time-independent parameters and symmetric couplings (under synchronous or asynchronous updates in the absence of lagged self-couplings), the kinetic Ising model results in an equilibrium process exhibiting a variety of complex phenomena, including ordered (ferromagnetic), disordered (paramagnetic), and quenched disordered states (known as spin glasses). The celebrated Sherrington-Kirkpatrick (SK) model, characterized by quenched random couplings resulting in a spin-glass phase15, can be solved using the replica mean-field method16. A kinetic version of this symmetric-coupling model has been represented as a bipartite network, also solved using the replica trick17.

The kinetics of equilibrium Ising systems are indistinguishable when observed in a forward or backward direction in time, i.e., they are invariant under the reversal of the arrow of time. This time-symmetry breaks down under time-varying external fields or asymmetric couplings comprising history-dependent, non-conservative forces. Such time-asymmetric processes violate detailed balance, leading to nonequilibrium dynamics yielding a positive entropy production8,18,19,20. In the latter case of asymmetric couplings with constant fields, the system may relax towards a steady state known as a nonequilibrium steady state after some time. Time-asymmetric trajectories in steady state are linked with entropy change of the heat baths under ‘local detailed balance’ for a system coupled to equilibrium reservoirs or heat baths21,22, suggesting that steady-state entropy production is critical for unveiling the interaction of out-of-equilibrium systems with their environments. Yet, unlike its equilibrium counterpart, the properties of irreversible Ising dynamics remain unclear due to the lack of theoretical description of its entropy production.

Here, we study the kinetics of the SK model with asymmetric connections under synchronous and asynchronous updates as a prototypical model of nonlinear and nonequilibrium processes. As the model does not have a free energy defined in classical terms, we resort to a dynamical equivalent in the form of a generating functional. We apply a path integral approach to obtain exact solutions on its statistical moments and nonequilibrium thermodynamic properties. Unlike the replica method, the generating functional for fully asymmetric couplings has exact solutions in the thermodynamic limit without additional assumptions like analytic continuation and replica symmetry breaking23. Previous studies using this method24,25 have shown that the asymmetric kinetic Ising model with asynchronous updates does not have a spin glass phase. In this manuscript, we will extend the generating functional path integral method to confirm this result in both cases of synchronous and asynchronous updates and further retrieve an exact solution of the entropy production of the system.

One of the open questions in empirical studies is whether an increase in entropy production observed in specific nonequilibrium systems under investigation is linked with the critical properties of systems approaching continuous phase transitions5,6. Entropy production is not necessarily maximized under such conditions and can display a continuous change26, or discontinuities in its derivative27. However, a number of simple nonequilibrium systems maximize their entropy production at a critical point. Examples are the entropy production of an Ising model with an oscillatory field and a mean-field majority vote model28,29,30. It is therefore important to investigate the case of the kinetic Ising system as a general model of physical and biological networks. We previously showed that the entropy production of the stationary asymmetric SK model with finite size takes a maximum around a critical point by applying mean-field approximations preserving fluctuations in the system31. However, this result (and the aforementioned references) relies on approximations and numerical simulations. Therefore, the assumption that entropy production is maximized near continuous phase transitions has not yet been ratified by exact solutions of spin models in the thermodynamic limit. In this study, we show analytically that the entropy production is locally maximized at critical phase transition points, representing a potentially useful phase transition correlate for systems without a globally defined free energy or heat capacity. Nevertheless, we also show that entropy production can take larger values for largely heterogeneous couplings in low-temperature regimes exhibiting disordered but nearly deterministic dynamics. Thus, entropy production must be examined carefully, as its increase does not necessarily indicate that the system is approaching a critical state. Instead, combining the entropy rate and entropy production yields a more precise picture of the irreversible processes.

The paper is organized as follows. First, we introduce maximum entropy Markov processes, their entropy production, and a generating functional method used to compute the system’s moments and the entropy production in both discrete and continuous time. Next, we describe the asymmetric SK model with synchronous and asynchronous updates and a path integral method calculating the configurational average of the generating functional. This yields an exact solution of the entropy production, magnetization, and correlations in an infinite system. We employ our theoretical results to draw phase diagrams of the order parameters and entropy production for synchronous and asynchronous dynamics with and without randomly sampled external fields. The theoretical predictions are then corroborated by numerical simulation. We also examine the critical line of nonequilibrium phase transitions, the temporal structure of the dynamics, and their relations to the entropy production. Finally, we conclude the paper by discussing the implications of our results for the study of biological systems.

Results

Maximum entropy Markov chains

The principle of maximum entropy is a foundation of equilibrium statistical mechanics32. The principle has been later generalized for treating time-dependent phenomena, as the principle of maximum caliber or maximum path entropy33,34. Under consistency requirements preserving causal interactions, the maximum caliber principle yields a Markov process35. To see this, we start with a discrete-time stochastic process with N discrete-state elements defined at time u as su = {s1,u, …, sN,u} for discrete-time trajectories of length t + 1 defined by a path probability p(s0:t). Later, we will show this discrete-time formulation can be generalized to an equivalent continuous-time formulation under appropriate assumptions.

Path entropy is defined as

$${S}_{0:t}=-\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{0:t}}p({{{{{{{{\bf{s}}}}}}}}}_{0:t})\log p({{{{{{{{\bf{s}}}}}}}}}_{0:t}).$$
(1)

Maximizing Eq. (1), subject to constraints, yields the least structured distribution p(s0:t) consistent with observations36. In causal network models, entropy maximization has to be constrained with a set of temporal consistency requirements35, as was first established by Kolmogorov37. Specifically, for any positive integer u(≤t), we impose

$$\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{u}}{p}_{0:u}({{{{{{{{\bf{s}}}}}}}}}_{0:u})={p}_{0:u-1}({{{{{{{{\bf{s}}}}}}}}}_{0:u-1}),$$
(2)

where p0:u(s0:u) is given by

$${p}_{0:u}({{{{{{{{\bf{s}}}}}}}}}_{0:u})=\underset{p({{{{{{{{\bf{s}}}}}}}}}_{0:u})}{\arg\max}\,{S}_{0:u}.$$
(3)

That is, we impose consistency between the marginal distribution for the maximum entropy path s0:u−1 in p0:u(s0:u) and the maximum entropy distribution of path s0:u−1, p0:u−1(s0:u−1). This constrains path distribution dependencies between consecutive states. We will drop the subscript in the path probability when not needed.

Maximizing Eq. (1) with constraints fn(su, su−1) = Cn,u (where Cn,u is a constant for the n-th constraint at time u), an initial distribution p(s0), and Eq. (2) results in a Markovian process35

$$p({{{{{{{{\bf{s}}}}}}}}}_{0:t})= p({{{{{{{{\bf{s}}}}}}}}}_{0})\mathop{\prod }\limits_{u=1}^{t}p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})\\ \propto p({{{{{{{{\bf{s}}}}}}}}}_{0})\mathop{\prod }\limits_{u=1}^{t}\exp \left[\mathop{\sum}\limits_{n}{\lambda }_{n}{f}_{n}({{{{{{{{\bf{s}}}}}}}}}_{u},\,{{{{{{{{\bf{s}}}}}}}}}_{u-1})\right].$$
(4)

The path entropy can be then decomposed into

$${S}_{0:t}= -\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{0:t}}p({{{{{{{{\bf{s}}}}}}}}}_{0:t})\left(\mathop{\sum}\limits_{u}\log p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})+\log p({{{{{{{{\bf{s}}}}}}}}}_{0})\right)\\= \mathop{\sum}\limits_{u}{S}_{u|u-1}+{S}_{0},$$
(5)

where S0 is the entropy of the initial distribution and Su∣u−1 is a conditional entropy, defined as

$${S}_{u|u-1}=-\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1}}p({{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1})\log p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1}),$$
(6)

which, at the steady state described in the following, corresponds to the Kolmogorov-Sinai entropy or entropy rate, \(\mathop{\lim }\nolimits_{t\to \infty }\frac{1}{t}{S}_{0:t}\).

Nonequilibrium steady state

A Markov chain converges to a unique stationary distribution if the system is irreducible (all states are accessible from any state in finite time) and aperiodic (the greatest common divisor of the number of steps for returning to the same state with non-zero probability is one38). We can confirm that these requirements are satisfied by Eq. (4) with finite transition probabilities, thus warranting the existence of a steady-state distribution π(su), which can be either in or out of thermodynamic equilibrium, as explained in the following.

For a discrete-time Markov chain, the evolution of the state probability distribution follows a master equation:

$${p}_{u}({{{{{{{{\bf{s}}}}}}}}}_{u})= \mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{u-1}}p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1}){p}_{u-1}({{{{{{{{\bf{s}}}}}}}}}_{u-1})\\= {p}_{u-1}({{{{{{{{\bf{s}}}}}}}}}_{u})+\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{u-1}}{j}_{{{{{{{{{\bf{s}}}}}}}}}_{u-1}\to {{{{{{{{\bf{s}}}}}}}}}_{u}}^{u}.$$
(7)

Here pv(su) is a marginal probability distribution of a state su calculated for the distribution at time v. For simplicity, we will omit the subscript and write p(su) when v = u. \({j}_{{{{{{{{{\bf{s}}}}}}}}}_{u-1}\to {{{{{{{{\bf{s}}}}}}}}}_{u}}^{u}\) are the system’s probability fluxes:

$${j}_{{{{{{{{{\bf{s}}}}}}}}}_{u-1}\to {{{{{{{{\bf{s}}}}}}}}}_{u}}^{u}\equiv p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})p({{{{{{{{\bf{s}}}}}}}}}_{u-1})-p({{{{{{{{\bf{s}}}}}}}}}_{u-1}|{{{{{{{{\bf{s}}}}}}}}}_{u}){p}_{u-1}({{{{{{{{\bf{s}}}}}}}}}_{u}).$$
(8)

In the limit of small probability fluxes, the system can be described by an equivalent continuous-time process:

$$\frac{dp({{{{{{{\bf{s}}}}}}}},t)}{dt}=\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }}{j}_{{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\to {{{{{{{\bf{s}}}}}}}}}(t)$$
(9)
$${j}_{{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\to {{{{{{{\bf{s}}}}}}}}}(t)\equiv w({{{{{{{\bf{s}}}}}}}}|{{{{{{{{\bf{s}}}}}}}}}^{{\prime} })p({{{{{{{{\bf{s}}}}}}}}}^{{\prime} },t)-w({{{{{{{{\bf{s}}}}}}}}}^{{\prime} }|{{{{{{{\bf{s}}}}}}}})p({{{{{{{\bf{s}}}}}}}},t),$$
(10)

where t refers to the continuous time and \(w({{{{{{{\bf{s}}}}}}}}|{{{{{{{{\bf{s}}}}}}}}}^{{\prime} })\) are transition rates.

The system is stationary or in a steady state if the sum of all probability fluxes is zero for all su, i.e., \({\sum }_{{{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1}}\,{j}_{{{{{{{{{\bf{s}}}}}}}}}_{u-1}\to {{{{{{{{\bf{s}}}}}}}}}_{u}}^{u}=0\). In addition, this will be an equilibrium steady state if \({j}_{{{{{{{{{\bf{s}}}}}}}}}_{u-1}\to {{{{{{{{\bf{s}}}}}}}}}_{u}}^{u}=0\) for all pairs su−1, su, resulting in the detailed balance condition

$$p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})\pi ({{{{{{{{\bf{s}}}}}}}}}_{u-1})=p({{{{{{{{\bf{s}}}}}}}}}_{u-1}|{{{{{{{{\bf{s}}}}}}}}}_{u})\pi ({{{{{{{{\bf{s}}}}}}}}}_{u}),$$
(11)

where π(su) is the steady-state distribution. When detailed balance is broken under the stationary condition, i.e., some \({j}_{{{{{{{{{\bf{s}}}}}}}}}_{u-1}\to {{{{{{{{\bf{s}}}}}}}}}_{u}}^{u}\ne 0\) but their sum is equal to zero, the stationary system is in a nonequilibrium steady state.

Steady-state entropy production

Stochastic thermodynamics describes a link between the time-irreversible stochastic trajectories with surroundings in the form of heat (entropy) dissipation. As the system evolves, it experiences an entropy change \({\sigma }_{u}^{{{{{{{{\rm{sys}}}}}}}}}\):

$${\sigma }_{u}^{{{{{{{{\rm{sys}}}}}}}}}={S}_{u}-{S}_{u-1}=\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1}}p({{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1})\log \frac{p({{{{{{{{\bf{s}}}}}}}}}_{u-1})}{p({{{{{{{{\bf{s}}}}}}}}}_{u})}.$$
(12)

Nonequilibrium systems maintain irreversible dynamics by continuously dissipating heat (entropy) to their environments. Under local detailed balance21,22,39 in a system coupled to a heat bath, the entropy change results from subtracting the entropy dissipated to the heat bath \({\sigma }_{u}^{{{{{{{{\rm{bath}}}}}}}}}\) from the (total) entropy production σu:

$${\sigma }_{u}^{{{{{{{{\rm{sys}}}}}}}}}={\sigma }_{u}-{\sigma }_{u}^{{{{{{{{\rm{bath}}}}}}}}},$$
(13)

where the entropy change of the heat bath is given as

$${\sigma }_{u}^{{{{{{{{\rm{bath}}}}}}}}}=\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1}}p({{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1})\log \frac{p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})}{p({{{{{{{{\bf{s}}}}}}}}}_{u-1}|{{{{{{{{\bf{s}}}}}}}}}_{u})}.$$
(14)

Here p(su−1∣su) is a transition probability (from Eq. (4)) but evaluated by the reverse trajectory21,40, that is, we define it using the transition function at time u, but switch su and su−1. This equation relates the system’s time asymmetry with the entropy change of the reservoir.

The entropy production σu at time u is then given as

$${\sigma }_{u}=\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1}}p({{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1})\log \frac{p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})p({{{{{{{{\bf{s}}}}}}}}}_{u-1})}{p({{{{{{{{\bf{s}}}}}}}}}_{u-1}|{{{{{{{{\bf{s}}}}}}}}}_{u})p({{{{{{{{\bf{s}}}}}}}}}_{u})},$$
(15)

which is the Kullback–Leibler divergence between the forward and backward trajectories8,18,20,41. Due to the non-negativity of the divergence, the entropy production is non-negative, σu ≥ 0. This entropy production vanishes if the probability of forward trajectories is identical to a posterior of past states given the future state20, i.e., when the process loses time-asymmetry in prediction and postdiction42.

Alternatively, the dissipation function8,43,44 quantifies the difference between incoming and outgoing fluxes in Eq. (8):

$${\widetilde{\sigma }}_{u}=\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1}}p({{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1})\log \frac{p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})p({{{{{{{{\bf{s}}}}}}}}}_{u-1})}{p({{{{{{{{\bf{s}}}}}}}}}_{u-1}|{{{{{{{{\bf{s}}}}}}}}}_{u}){p}_{u-1}({{{{{{{{\bf{s}}}}}}}}}_{u})},$$
(16)

which directly assesses a violation of the detailed balance. The entropy production σu and dissipation function \({\widetilde{\sigma }}_{u}\) are equivalent under steady-state conditions. Furthermore, both quantities become equivalent in the continuous-time limit8,42,44 and converge to the entropy production rate39:

$$\frac{d\sigma (t)}{dt}=\frac{1}{2}\mathop{\sum}\limits_{{{{{{{{\bf{s}}}}}}}},{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }}{j}_{{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\to {{{{{{{\bf{s}}}}}}}}}(t)\log \frac{w({{{{{{{\bf{s}}}}}}}}|{{{{{{{{\bf{s}}}}}}}}}^{{\prime} })p({{{{{{{{\bf{s}}}}}}}}}^{{\prime} },t)}{w({{{{{{{{\bf{s}}}}}}}}}^{{\prime} }|{{{{{{{\bf{s}}}}}}}})p({{{{{{{\bf{s}}}}}}}},t)}.$$
(17)

In a steady state, the entropy production is caused by dissipation only and becomes equivalent to the house-keeping entropy production caused by the non-conservative forces under a steady state45,46. Both σu and \({\widetilde{\sigma }}_{u}\) result in:

$${\sigma }_{u}={\widetilde{\sigma }}_{u}={\sigma }_{u}^{{{{{{{{\rm{bath}}}}}}}}}=-{S}_{u|u-1}+{S}_{u|u-1}^{r}.$$
(18)

Here \({S}_{u|u-1}^{r}\) is the entropy of the time-reversed conditional distribution:

$${S}_{u|u-1}^{r}\equiv -\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1}}p({{{{{{{{\bf{s}}}}}}}}}_{u},{{{{{{{{\bf{s}}}}}}}}}_{u-1})\log p({{{{{{{{\bf{s}}}}}}}}}_{u-1}|{{{{{{{{\bf{s}}}}}}}}}_{u}).$$
(19)

In this paper, we study the steady-state entropy production in Eq. (18), which is critical for evaluating the interaction of the nonequilibrium processes with their environment.

Generating functional

Consider a maximum caliber path probability (Eq. (4))

$$p({{{{{{{{\bf{s}}}}}}}}}_{0:t})=\mathop{\prod }\limits_{u=1}^{t}p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})p({{{{{{{{\bf{s}}}}}}}}}_{0}),$$
(20)

For simplicity, we will assume \({p}_{0}({{{{{{{\bf{s}}}}}}}})={\prod }_{i}\delta \left[{s}_{i},{s}_{i,0}\right]\) – the initial distribution is a Kronecker delta with a unique initial state – and ignore the term. However, the following steps are general for any p0(s0).

In equilibrium systems, the partition function retrieves their statistical moments. A nonequilibrium equivalent function is a generating functional or dynamical partition function. To obtain not only the statistical properties averaged over trajectories, but also the forward/time-reversed conditional entropies (Eqs. (6), (19)), we define the following generating functional:

$${Z}_{t}({{{{{{{\bf{g}}}}}}}})=\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{0:t}}p({{{{{{{{\bf{s}}}}}}}}}_{0:t})\exp \left[\Gamma ({{{{{{{\bf{g}}}}}}}},{{{{{{{{\bf{s}}}}}}}}}_{0:t})\right],$$
(21)
$$\Gamma ({{{{{{{\bf{g}}}}}}}},{{{{{{{{\bf{s}}}}}}}}}_{0:t})= \mathop{\sum}\limits_{i,u}{g}_{i,u}{s}_{i,u}-\mathop{\sum}\limits_{u}{g}_{u}^{S}\epsilon ({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})\\ -\mathop{\sum}\limits_{u}{g}_{u}^{{S}^{r}}\epsilon ({{{{{{{{\bf{s}}}}}}}}}_{u-1}|{{{{{{{{\bf{s}}}}}}}}}_{u}),$$
(22)

where \(\epsilon ({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})\equiv -\log p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})\). In the limit t → ∞, the logarithm of the generating functional converges to the large deviation function47,48,49,

$$\varphi ({{{{{{{\bf{g}}}}}}}})=\mathop{\lim }\limits_{t\to \infty }\frac{1}{t}\log {Z}_{t}({{{{{{{\bf{g}}}}}}}}),$$
(23)

which plays the role of a free-energy function for nonequilibrium trajectories50. The vector g is composed of parameters gi,u (i = 1, …, N, u = 1, …, t) and \({g}_{u}^{S},{g}_{u}^{{S}^{r}}\) (u = 1, …, t) retrieving the system’s statistical properties. The parameters gi,u recover the statistical moments of the systems like the average rates and correlations:

$$\mathop{\lim }\limits_{{{{{{{{\bf{g}}}}}}}}\to {{{{{{{\bf{0}}}}}}}}}\frac{\partial {Z}_{t}({{{{{{{\bf{g}}}}}}}})}{\partial {g}_{i,u}}=\mathop{\lim }\limits_{{{{{{{{\bf{g}}}}}}}}\to {{{{{{{\bf{0}}}}}}}}}{\left\langle {s}_{i,u}\right\rangle }_{{{{{{{{\bf{g}}}}}}}}}=\left\langle {s}_{i,u}\right\rangle,$$
(24)
$$\mathop{\lim }\limits_{{{{{{{{\bf{g}}}}}}}}\to {{{{{{{\bf{0}}}}}}}}}\frac{{\partial }^{2}{Z}_{t}({{{{{{{\bf{g}}}}}}}})}{\partial {g}_{i,u}\partial {g}_{j,v}}=\mathop{\lim }\limits_{{{{{{{{\bf{g}}}}}}}}\to {{{{{{{\bf{0}}}}}}}}}{\left\langle {s}_{i,u}{s}_{j,v}\right\rangle }_{{{{{{{{\bf{g}}}}}}}}}=\left\langle {s}_{i,u}{s}_{j,v}\right\rangle,$$
(25)

where angle brackets are defined as

$${\left\langle f({{{{{{{{\bf{s}}}}}}}}}_{0:t})\right\rangle }_{{{{{{{{\bf{g}}}}}}}}}=\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{0:t}}f({{{{{{{{\bf{s}}}}}}}}}_{0:t})\exp \left[\mathop{\sum}\limits_{i,u}{g}_{i,u}{s}_{i,u}\right]p({{{{{{{{\bf{s}}}}}}}}}_{0:t}),$$
(26)
$$\left\langle f({{{{{{{{\bf{s}}}}}}}}}_{0:t})\right\rangle=\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{0:t}}f({{{{{{{{\bf{s}}}}}}}}}_{0:t})p({{{{{{{{\bf{s}}}}}}}}}_{0:t}).$$
(27)

In addition, \({g}_{u}^{S},{g}_{u}^{{S}^{r}}\) retrieve the conditional and time-reversed conditional entropy terms, \({S}_{u|u-1},{S}_{u|u-1}^{r}\):

$${S}_{u|u-1}= -\mathop{\lim }\limits_{{{{{{{{\bf{g}}}}}}}}\to {{{{{{{\bf{0}}}}}}}}}\frac{\partial {Z}_{t}({{{{{{{\bf{g}}}}}}}})}{\partial {g}_{u}^{S}}\\= \mathop{\lim }\limits_{{{{{{{{\bf{g}}}}}}}}\to {{{{{{{\bf{0}}}}}}}}}{\left\langle \epsilon ({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})\right\rangle }_{{{{{{{{\bf{g}}}}}}}}}=\left\langle \epsilon ({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})\right\rangle,$$
(28)
$${S}_{u|u-1}^{r}= -\mathop{\lim }\limits_{{{{{{{{\bf{g}}}}}}}}\to {{{{{{{\bf{0}}}}}}}}}\frac{\partial {Z}_{t}({{{{{{{\bf{g}}}}}}}})}{\partial {g}_{u}^{{S}^{r}}}\\= \mathop{\lim }\limits_{{{{{{{{\bf{g}}}}}}}}\to {{{{{{{\bf{0}}}}}}}}}{\left\langle \epsilon ({{{{{{{{\bf{s}}}}}}}}}_{u-1}|{{{{{{{{\bf{s}}}}}}}}}_{u})\right\rangle }_{{{{{{{{\bf{g}}}}}}}}}=\left\langle \epsilon ({{{{{{{{\bf{s}}}}}}}}}_{u-1}|{{{{{{{{\bf{s}}}}}}}}}_{u})\right\rangle,$$
(29)

and thus the steady-state entropy production (Eq. (18)):

$${\sigma }_{u}=\mathop{\lim }\limits_{{{{{{{{\bf{g}}}}}}}}\to {{{{{{{\bf{0}}}}}}}}}\left(\frac{\partial {Z}_{t}({{{{{{{\bf{g}}}}}}}})}{\partial {g}_{u}^{S}}-\frac{\partial {Z}_{t}({{{{{{{\bf{g}}}}}}}})}{\partial {g}_{u}^{{S}^{r}}}\right).$$
(30)

Synchronous and asynchronous, asymmetric Sherrington-Kirkpatrick model

We consider N interacting elements su (spins or neurons), taking each element i at time u a binary state si,u = { − 1, 1}. Constraints take the form of delayed pairwise couplings (i.e., fij(su, su−1) = si,usj,u−1 in Eq. (4)). This results in the dynamics:

$$p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})=\mathop{\prod}\limits_{i}\frac{\exp \left[\beta {s}_{i,u}{h}_{i,u}\right]}{2\cosh \left[\beta {h}_{i,u}\right]},$$
(31)
$${h}_{i,u}={H}_{i,u}+\mathop{\sum}\limits_{j}{J}_{ij}{s}_{j,u-1},$$
(32)

where β is the inverse temperature. The system’s state at time u depends on the previous time-step (Fig. 1a).

Fig. 1: Asymmetric kinetic SK model.
figure 1

a The asymmetric kinetic Ising model describes a Markov chain where states at time su depend on pairwise couplings to states su−1. This model shows disordered dynamics for large coupling variance both at high and low temperatures (b and c), ordered dynamics for low temperatures and a low coupling variance (d), and critical dynamics at the phase transition (e).

Full size image

The equation above is a general formulation of a kinetic Ising model with time-dependent fields Hi,u. The dynamics can include both synchronous and asynchronous Ising systems by introducing a set of independent Bernoulli random variables: τi,u = 0, 1 with probabilities 1 − α and α (i.e., \({\tau }_{i,u} \sim {{{{{{{\rm{Bernoulli}}}}}}}}\left(\alpha \right)\)) and making Hi,u stochastic processes:

$${H}_{i,u}={\Theta }_{i,u}+(1-{\tau }_{i,u})K{s}_{i,u-1}.$$
(33)

Note that in the limit of K → ∞, the state si,u is tightly coupled to the previous state si,u−1. Therefore, the state changes only if τi,u = 1. We have the following transition probability in the K → ∞ limit:

$$p({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})= \mathop{\prod}\limits_{i}\left({\tau }_{i,u}w({s}_{i,u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})\right.\\ \left.+(1-{\tau }_{i,u})\delta \left[{s}_{i,u},{s}_{i,u-1}\right]\right),$$
(34)

with the transition rate

$$w({s}_{i,u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})=\frac{\exp \left[\beta {s}_{i,u}{h}_{i,u}^{1}\right]}{2\cosh \left[\beta {h}_{i,u}^{1}\right]},$$
(35)

where \({h}_{i,u}^{1}={\Theta }_{i,u}+{\sum }_{j}{J}_{ij}{s}_{j,u-1}\) and Θi,u is an external field. With α = 1 we have a kinetic Ising system under parallel or synchronous updates. In the limit α → 0, we have in turn a kinetic Ising system with asynchronous updates (i.e., at most one spin is updated each time step), converging to a continuous-time master equation.

The generating functional of the kinetic Ising system (Eq. (21)) is defined by the functions

$$\epsilon ({{{{{{{{\bf{s}}}}}}}}}_{u}|{{{{{{{{\bf{s}}}}}}}}}_{u-1})=-\mathop{\sum}\limits_{i}\left(\beta {s}_{i,u}{h}_{i,u}+\log \left[2\cosh \left[\beta {h}_{i,u}\right]\right]\right),$$
(36)
$$\epsilon ({{{{{{{{\bf{s}}}}}}}}}_{u-1}|{{{{{{{{\bf{s}}}}}}}}}_{u})=-\mathop{\sum}\limits_{i}\left(\beta {s}_{i,u-1}{h}_{i,u}^{r}+\log \left[2\cosh \left[\beta {h}_{i,u}^{r}\right]\right]\right),$$
(37)

where \({h}_{i,u}^{r}={H}_{i,u}+{\sum }_{j}{J}_{ij}{s}_{j,u}={h}_{i,u+1}+{H}_{i,u}-{H}_{i,u+1}\).

The equilibrium Ising model with symmetric random Gaussian couplings is referred to as the SK model. In the fully-asymmetric SK model, the couplings Jij are quenched independent variables, each following a Gaussian distribution

$$p({J}_{ij})=\frac{1}{\sqrt{2\pi \Delta {J}^{2}/N}}\exp \left[\frac{-1}{2\Delta {J}^{2}/N}{\left({J}_{ij}-\frac{{J}_{0}}{N}\right)}^{2}\right],$$
(38)

with mean J0/N and variance ΔJ2/N scaled by N.

The asymmetric SK model shows a variety of population dynamics. Figure 1 shows exemplary dynamics under asynchronous updates without fields (Θi,u = 0). It shows disordered dynamics for large coupling variance ΔJ2 both at high and low temperatures (i.e. low and large β, Fig. 1b, c), ordered dynamics for low temperatures and low ΔJ2 (Fig. 1d), and critical dynamics at the phase transition (Fig. 1e).

Solution of the asymmetric Sherrington-Kirkpatrick model

The solution of the kinetic version of the SK model with asymmetric and quenched couplings can be obtained by computing the generating functional averaged over the couplings (referred to as the configurational average):

$${[{Z}_{t}({{{{{{{\bf{g}}}}}}}})]}_{{{{{{{{\bf{J}}}}}}}}}=\int\,d{{{{{{{\bf{J}}}}}}}}{Z}_{t}({{{{{{{\bf{g}}}}}}}})\mathop{\prod}\limits_{i,j}p({J}_{ij}).$$
(39)

This integral cannot be solved directly because of the \(\log \left[2\cosh \left[\cdot \right]\right]\) terms in Eqs. (36) and (37), which depend nonlinearly on Jij. A path integral method51 to find a solution introduces a delta integral representing βhi,u with auxiliary variables θi,u = β(Hi,u + ∑jJijsj,u−1) as well as \(\beta {h}_{i,u}^{r}\) with an auxiliary variable ϑi,u = θi,u+1 + β(Hi,u − Hi,u+1). Let θ = {θi,u} (note u = 1, …, t + 1) and ϑ = {ϑi,u} (u = 0, …, t) denote a set of the auxiliary variables. Using conjugate variables \(\hat{{{{{{{{\boldsymbol{\theta }}}}}}}}}=\{\hat{\theta }_{i,u}\}\) to represent the delta function in the integral form, the configurational average is written as

$${[{Z}_{t}({{{{{{{\bf{g}}}}}}}})]}_{{{{{{{{\bf{J}}}}}}}}}= \frac{1}{{\left(2\pi \right)}^{N(t+1)}}\int\,d{{{{{{{\boldsymbol{\theta }}}}}}}}d\hat{{{{{{{{\boldsymbol{\theta }}}}}}}}}d{{{{{{{\bf{J}}}}}}}}\left(\mathop{\prod}\limits_{i,j}p({J}_{ij})\right)\\ \cdot \mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{1:t}}\exp \left[\mathop{\sum}\limits_{i,u}({s}_{i,u}{\theta }_{i,u}-\log \left[2\cosh {\theta }_{i,u}\right])\right.\\ \left.+\mathop{\sum}\limits_{i,u}{{{{{{{\rm{i}}}}}}}}{\hat{\theta }}_{i,u}\left({\theta }_{i,u}-\beta {H}_{i,u}-\beta \mathop{\sum}\limits_{j}{J}_{ij}{s}_{j,u-1}\right)+\Gamma ({{{{{{{\bf{g}}}}}}}},{{{{{{{{\bf{s}}}}}}}}}_{0:t},{{{{{{{\boldsymbol{\theta }}}}}}}},{{{{{{{\bf{\vartheta }}}}}}}})\right],$$
(40)

with

$$\Gamma ({{{{{{{\bf{g}}}}}}}},{{{{{{{{\bf{s}}}}}}}}}_{0:t},{{{{{{{\boldsymbol{\theta }}}}}}}},{{{{{{{\bf{\vartheta }}}}}}}})= \mathop{\sum}\limits_{i,u}{\Gamma }_{i,u}({{{{{{{\bf{g}}}}}}}},{{{{{{{{\bf{s}}}}}}}}}_{0:t},{\theta }_{i,u},{\vartheta }_{i,u})\\= \mathop{\sum}\limits_{i,u}\left({g}_{i,u}{s}_{i,u}+{g}_{u}^{S}\left({s}_{i,u}{\theta }_{i,u}-\log \left[2\cosh \left[{\theta }_{i,u}\right]\right]\right)\right.\\ \left.+{g}_{u}^{{S}^{r}}\left({s}_{i,u-1}{\vartheta }_{i,u}-\log \left[2\cosh \left[{\vartheta }_{i,u}\right]\right]\right)\right).$$
(41)

Note that the summation of \({\hat{\theta }}_{i,u}\) terms is performed over u = 1, …, t + 1 to retrieve the fields of both the forward and backward trajectories.

The integral over Jij can be now performed directly over linear exponential terms (see Supplementary Note 1). After integration, Eq. (40) incorporates quadruple-wise interactions among spins s0:t and conjugate variables \(\hat{{{{{{{{\boldsymbol{\theta }}}}}}}}}\) (Supplementary Eq. (S1.10)), similar to replica interactions in the equilibrium SK model12. These interactions are simplified by introducing Gaussian integrals and a saddle-point approximation in the thermodynamic limit (Supplementary Eq. (S1.26)). The saddle-point solution can be simplified by introducing four types of order parameters (Supplementary Eq. (S1.31)). In fully-asymmetric networks, two of these order parameters are found to be zero at g = 0, yielding a solution in terms of the order parameters mu and qu,v (see Supplementary Eqs. (S1.48) and (S1.49)):

$${m}_{u}=\frac{1}{N}\mathop{\sum}\limits_{i}{\left[\left\langle {s}_{i,u}\right\rangle \right]}_{{{{{{{{\bf{J}}}}}}}}},$$
(42)
$${q}_{u,v}=\frac{1}{N}\mathop{\sum}\limits_{i}{\left[\left\langle {s}_{i,u}{s}_{i,v}\right\rangle \right]}_{{{{{{{{\bf{J}}}}}}}}}.$$
(43)

Finally, conjugate variables \(\hat{{{{{{{{\boldsymbol{\theta }}}}}}}}}\) in the saddle-point solution can be substituted with a multivariate Gaussian integral (Supplementary Eq. (S1.57)), leading to a factorized generating functional

$${[{Z}_{t}({{{{{{{\bf{g}}}}}}}})]}_{{{{{{{{\bf{J}}}}}}}}}= \mathop{\prod}\limits_{i}\mathop{\sum}\limits_{{{{{{{{{\bf{s}}}}}}}}}_{i,1:t}}\int\,d{{{{{{{\boldsymbol{\xi }}}}}}}}p({{{{{{{\boldsymbol{\xi }}}}}}}})\exp \left[\mathop{\sum}\limits_{u}{s}_{i,u}{\overline{h}}_{i,u}({\xi }_{u})\right.\\ + \mathop{\sum}\limits_{u}\beta {s}_{i,u-1}{\widetilde{h}}_{i,u-1}-\mathop{\sum}\limits_{u}\log 2\cosh \left[\beta {\overline{h}}_{i,u}({\xi }_{u})\right]\\ \left.+\mathop{\sum}\limits_{u}{\Gamma }_{i,u}({{{{{{{\bf{g}}}}}}}},{{{{{{{{\bf{s}}}}}}}}}_{0:t},\beta {\overline{h}}_{i,u}({\xi }_{u}),\beta {\overline{h}}_{i,u}^{r}({\xi }_{u+1}))\right],$$
(44)

where the stochastic elements ξ = (ξ1, …, ξt+1) affecting each spin i follow a multivariate normal distribution \(p({{{{{{{\boldsymbol{\xi }}}}}}}})={{{{{{{\mathcal{N}}}}}}}}({{{{{{{\bf{0}}}}}}}},{{{{{{{\bf{q}}}}}}}})\) with q defining qu−1,v−1 as the covariance of each pair ξu, ξv for u, v ∈ 1, …, t + 1. Here, at g = 0, spin interactions are effectively substituted by same-spin temporal couplings in mean effective fields

$${\overline{h}}_{i,u}({\xi }_{u})={H}_{i,u}+{J}_{0}{m}_{u-1}+\Delta J{\xi }_{u},$$
(45)
$${\overline{h}}_{i,u}^{r}({\xi }_{u+1})={H}_{i,u}+{J}_{0}{m}_{u}+\Delta J{\xi }_{u+1},$$
(46)
$${\widetilde{h}}_{i,u-1}=0.$$
(47)

Applying Eqs. (24) and (25) to the configurational average in Eq. (44), we obtain the order parameters mu and qu,v:

$${m}_{u}=\frac{1}{N}\mathop{\sum}\limits_{i}\int\,{{{{{{{\rm{D}}}}}}}}z\tanh \left[\beta {\overline{h}}_{i,u}(z)\right],$$
(48)
$${q}_{u,v}=\frac{1}{N}\mathop{\sum}\limits_{i}\int\,{{{{{{{\rm{D}}}}}}}}x{y}^{({q}_{u-1,v-1})}\tanh \left[\beta {\overline{h}}_{i,u}(x)\right] \cdot \tanh \left[\beta {\overline{h}}_{i,v}(y)\right],$$
(49)

where the Gaussian stochastic terms are simplified to

$${{{{{{{\rm{D}}}}}}}}z=\frac{1}{\sqrt{2\pi }}\exp \left[-\frac{1}{2}{z}^{2}\right],$$
(50)
$${{{{{{{\rm{D}}}}}}}}x{y}^{({q}_{u,v})}=\frac{1}{2\pi \sqrt{1-{q}_{u,v}^{2}}}\exp \left[\frac{-{x}^{2}-{y}^{2}+2{q}_{u,v}xy}{2(1-{q}_{u,v}^{2})}\right].$$
(51)

In contrast with the equilibrium SK model, mu is independent of qu,v, resulting in the lack of a spin-glass phase as suggested by previous studies25.

The configurational average of Eqs. (28) and (29) results in the following conditional entropy and time-reversed conditional entropy

$${\left[{S}_{u|u-1}\right]}_{{{{{{{{\bf{J}}}}}}}}}= \mathop{\sum}\limits_{i}\int-{{{{{{{\rm{D}}}}}}}}z\left(\beta \left({H}_{i,u}+{J}_{0}{m}_{u-1}\right)\tanh \left[\beta {\overline{h}}_{i,u}(z)\right]\right.\\ + {\beta }^{2}\Delta {J}^{2}\left(1-{\tanh }^{2}\left[\beta {\overline{h}}_{i,u}(z)\right]\right)\\ \left.-\log \left[2\cosh \left[\beta {\overline{h}}_{i,u}(z)\right]\right]\right),$$
(52)
$${\left[{S}_{u|u-1}^{r}\right]}_{{{{{{{{\bf{J}}}}}}}}}= \mathop{\sum}\limits_{i}\int-{{{{{{{\rm{D}}}}}}}}z\left(\beta ({H}_{i,u}+{J}_{0}{m}_{u})\tanh \left[\beta {\overline{h}}_{i,u-1}(z)\right]\right.\\ + {\beta }^{2}\Delta {J}^{2}{q}_{u,u-2}\left(1-{\tanh }^{2}\left[\beta {\overline{h}}_{i,u-1}(z)\right]\right)\\ \left.-\log \left[2\cosh \left[\beta {\overline{h}}_{i,u}^{r}(z)\right]\right]\right).$$
(53)

Up to this point, our results are general for time-dependent fields Hi,u, covering synchronous and asynchronous updates by Eq. (33). We obtain the results for the synchronous SK model by setting α = 1 or, equivalently, Hi,u = Θi,u. For time-independent fields (Θi,u = Θi), the system converges to a steady state determined by the solution of the self-consistent equations given by Eqs. (48) and (49). Finally, using Eq. (30), the steady-state entropy production under the synchronous updates is obtained as

$${\left[{\sigma }_{u}\right]}_{{{{{{{{\bf{J}}}}}}}}}={\beta }^{2}\Delta {J}^{2}(1-{q}_{u,u-2})\mathop{\sum}\limits_{i}\int\,{{{{{{{\rm{D}}}}}}}}z(1-{\tanh }^{2}\left[\beta \left({\Theta }_{i}+{J}_{0}{m}_{u-1}+\Delta Jz\right)\right]),$$
(54)

with mu−1 and qu,u−2 given by their steady-state values (i.e., independent of u). Note that for the synchronous system the steady-state solution of qu,v is the same for all u, v. In the following, we will use m and q to represent these steady-state solutions.

To calculate the steady-state solutions for the asynchronous SK model, we calculate the generating functional \({\left[{Z}_{t}({{{{{{{\bf{g}}}}}}}})\right]}_{{{{{{{{\bf{J}}}}}}}},{{{{{{{\boldsymbol{\tau }}}}}}}}}\) that is additionally averaged over the independent random variables τi,u in Eq. (33). We show in Supplementary Note 2 that the resulting order parameters in continuous-time m(t) and \(q({t}^{{\prime} },t)\) are subject to the following dynamical equations:

$$\frac{dm(t)}{dt}=\frac{1}{N}\mathop{\sum}\limits_{i}\int\,{{{{{{{\rm{D}}}}}}}}z\tanh \left[\beta {h}_{i}^{*}(z,t)\right]-m(t).$$
(55)
$$\frac{dq({t}^{{\prime} },t)}{dt}={q}^{1}({t}^{{\prime} },t)-q({t}^{{\prime} },t),$$
(56)
$$\frac{d{q}^{1}({t}^{{\prime} },t)}{d{t}^{{\prime} }}= \frac{1}{N}\mathop{\sum}\limits_{i}\int\,{{{{{{{\rm{D}}}}}}}}x{y}^{(q({t}^{{\prime} },\,t))}\tanh \left[\beta {h}_{i}^{*}(x,\,{t}^{{\prime} })\right]\\ \cdot \tanh \left[\beta {h}_{i}^{*}(y,\,t)\right]-{q}^{1}({t}^{{\prime} },\,t)$$
(57)

with \({h}_{i}^{*}(z,t)={\Theta }_{i}+{J}_{0}m(t)+\Delta Jz\). Here \({q}^{1}({t}^{{\prime} },t)\) is the spin correlation conditioned on spins being updated at time t. The steady-state solutions of m(t) and \(q({t}^{{\prime} },t)\) (assuming \({t}^{{\prime} }\gg t\)) converge to the same steady-state values m and q found for the synchronous SK model (see Supplementary Note 2).

In the continuous-time limit, the steady-state entropy production converges to a steady-state entropy production rate (Eq. (17)) given by:

$${\left[\frac{d\sigma (t)}{dt}\right]}_{{{{{{{{\bf{J}}}}}}}},{{{{{{{\boldsymbol{\tau }}}}}}}}}= \mathop{\lim }\limits_{\alpha \to 0}{\beta }^{2}\Delta {J}^{2}(1-q(t+\alpha,t-\alpha ))\\ \cdot \mathop{\sum}\limits_{i}\int\,{{{{{{{\rm{D}}}}}}}}z\left(1-{\tanh }^{2}\left[\beta {h}_{i}^{*}(z,t)\right]\right).$$
(58)

The delayed-self correlation \(q({t}^{{\prime} },t)\) is discontinuous at \({t}^{{\prime} }=t\) (i.e., \(\mathop{\lim }\limits_{\alpha \to 0}q(t+\alpha,t-\alpha )\ne q(t,t)=1\), see Supplementary Fig. S1), warranting that the entropy production rate can be non-zero for appropriate parameters.

Given the analytical solutions of the system, we will now study the phase diagram of the SK model. In contrast with the naive replica-symmetric solution of the equilibrium SK model, the equations above are exact in the model with asymmetric couplings in the thermodynamic limit under both synchronous and asynchronous updates.

The SK model without external fields

Figure 2a, b displays the phase diagram of the steady-state order parameters, m and q, for both synchronous and asynchronous updates, respectively derived from Eqs. (48) and (49) as a function of the inverse temperature β and the width of the coupling distribution ΔJ, when the external fields are fixed at zeros (Θi,u = 0) and the mean coupling is J0 = 1. In this setting, the inverse temperature β controls the magnitude of the couplings. The phase diagram shows two distinct regions, one in which the order parameters are fixed at zero (zero magnetization and zero self-correlations, m = 0 and q = 0) – indicating disordered states – and the other in which the order parameters become positive (m > 0 and q > 0) – indicating ordered states. Therefore, the system exhibits a nonequilibrium analogue of the paramagnetic-ferromagnetic (disorder-order) phase transition controlled by the parameters, β and ΔJ. The dashed line in each panel shows the critical values of ΔJ as a function of β, which is obtained by solving the following equation (see Supplementary Note 3),

$$\frac{1}{\beta {J}_{0}}=\int\,{{{{{{{\rm{D}}}}}}}}z\left(1-{\tanh }^{2}\left[\beta \left(\Delta Jz\right)\right]\right).$$
(59)

The solution will be denoted as ΔJc(β). As studied in Supplementary Note 3, this critical phase transition corresponds to the mean-field universality class, as in the order-disorder phase transition of the equilibrium SK model (note that the spin-glass phase has different exponents52).

Fig. 2: Order parameters of the asymmetric SK model with zero fields under synchronous and asynchronous updates.
figure 2

The average magnetization m and the average delayed self-coupling q are shown in the phase diagram of the inverse temperature β and coupling heterogeneity ΔJ using a model with fixed parameters J0 = 1, ΔH = 0. The dashed line represents the critical line separating ordered and disordered phases. The dotted line represents the critical value at zero temperature (β → ∞).

Full size image

Depending on the coupling variance ΔJ, the dynamics do or do not undergo the nonequilibrium phase transition by varying the inverse temperature β. The critical ΔJc(β) at β → ∞ is given as ΔJc(∞) = 0.79501 (dotted horizontal line). If the distribution is narrower than the critical value ΔJc(∞), the process undergoes the phase transition by changing β. If the distribution is wider than the critical value, the order parameters are fixed at zeros (m = 0, q = 0) for any β. Note that, for β → ∞ (zero temperature), the activation function approaches the threshold nonlinearity given by the sign function; therefore, the process becomes deterministic. That is, for the large values of β, the process approaches deterministic dynamics yielding either ordered or disordered states for smaller or larger ΔJ, respectively. We remark that the disordered state with m = 0 and q = 0 at high β (low temperature) does not indicate the spin-glass phase as expected for the equilibrium Ising system (see Supplementary Note 4). We confirmed the non-existence of a spin-glass phase for the asymmetric kinetic SK model by finding that the system decays exponentially in this region (Supplementary Note 5).

The reduction in uncertainty at higher β is indicated by the reduction of the conditional entropy (the path entropy) Su∣u−1 by increasing β (Fig. 3a). This figure additionally shows that the conditional entropy decreases slowly with increasing β along the critical line of the phase transitions. This means that strong couplings and diverse patterns co-exist along the critical line. In contrast, the time-reversed conditional entropy \({S}_{u|u-1}^{r}\) (Fig. 3b) displays opposite dependency on β for the broader or narrower coupling distributions. Time-reversed conditional entropy quantifies how surprising the reverse process is under the forward model. With coupling distributions narrower than the critical value ΔJc(∞), the time-reversed conditional entropy diminishes by increasing β, indicating that the reverse processes takes place with increasingly high probabilities. This is because the spin state is fixated at all up or down under the ferromagnetic-like state for all time, losing temporal asymmetry. In contrast, the reverse process becomes less likely to happen as the dynamics becomes more deterministic by increasing β yet remains disordered as long as the coupling distribution is broader than ΔJc(∞). This distinct behavior between the conditional entropy and its time-reversed version found at the wider coupling distributions and high inverse temperatures yields the strong time-asymmetry in this regime.

Fig. 3: Steady-state entropy rate and entropy production of the asymmetric SK model under synchronous updates.
figure 3

a The phase diagram of the conditional entropy \({\left[{S}_{u|u-1}\right]}_{{{{{{{{\bf{J}}}}}}}}}\) (equivalent to the entropy rate) as a function of the inverse temperature β and the coupling heterogeneity ΔJ. b The conditional entropy of the reverse dynamics \({\left[{S}_{u|u-1}^{r}\right]}_{{{{{{{{\bf{J}}}}}}}}}\). c The entropy production at a steady state. The white dashed line is a critical line for the nonequilibrium phase transitions. d (inset) The horizontal sections of the entropy production (ΔJ = 0.4, 0.5, 0.6, 0.7, and 0.7950), showing that it peaks at the critical line. All figures are based on a model with fixed parameters Hi,u = 0 and J0 = 1.

Full size image

The entropy production under the steady-state condition quantifies the difference between the conditional and time-reversed conditional entropy. Figure 3c displays the phase diagram of the entropy production for the synchronous Ising model (the asynchronous Ising model has a similar behavior but different scaling, see Supplementary Note 2 or Fig. 4). The entropy production is maximized at the high β under the broader coupling distributions, where we find a significant difference between these two conditional entropies. Namely, strong time-asymmetry appears when the dynamics are disordered, nearly deterministic processes. The entropy production increases with β if the coupling distribution is wider than ΔJc(∞). In contrast, the entropy production is locally maximized at the critical point (white dashed line) with the coupling distribution being narrower than ΔJc(∞) (see also Fig. 3d). For the narrowly distributed couplings, the process exhibits a paramagnetic-like (randomized or disordered) phase at lower β and a ferromagnetic-like (ordered) phase at higher β (Fig. 2), neither of which can exhibit adequately asymmetric dynamics in time. Time-asymmetry appears between the ordered and disordered phases, namely at the critical point. As a consequence, the steady-state entropy production can be a measure of the criticality in this regime. However, more importantly, the magnitude of the entropy production is far more significant in the regime of large ΔJ and β than near the critical states, due to the strong time-asymmetry caused by the combination of disordered and quasi-deterministic dynamics.

Fig. 4: Verification of the exact mean-field solutions by simulating the kinetic Ising systems with synchronous and asynchronous updates.
figure 4

We repeated 400, 000 simulations of systems under synchronous (top) and asynchronous (bottom) updates with Θi,u = 0 and ΔJ = 0.5: a, d Sample estimates of the mean activation rate \(\frac{1}{N}{\sum }_{i}{\left[\left\langle {s}_{i,u}\right\rangle \right]}_{{{{{{{{\bf{J}}}}}}}},{{{{{{{\boldsymbol{\tau }}}}}}}}}\) compared with the theoretical order parameter m. b, e Sample estimates of the average delayed self-covariances \(\frac{1}{N}{\sum }_{i}{\left[\left\langle {s}_{i,u}{s}_{i,u-d}\right\rangle \right]}_{{{{{{{{\bf{J}}}}}}}},{{{{{{{\boldsymbol{\tau }}}}}}}}}- \frac{1}{N}{\sum }_{i}{\left[\left\langle {s}_{i,u}\right\rangle \right]}_{{{{{{{{\bf{J}}}}}}}},{{{{{{{\boldsymbol{\tau }}}}}}}}} \frac{1}{N}{\sum }_{i}{\left[\left\langle {s}_{i,u-d}\right\rangle \right]}_{{{{{{{{\bf{J}}}}}}}},{{{{{{{\boldsymbol{\tau }}}}}}}}} \) (d = 1 for the synchronous system and d = 10N for the asynchronous one) computed from samples compared with the theoretical order parameter q-m2. c, f Sample estimates of the entropy production and entropy production rate (Supplementary Eqs. (S6.6) and (S6.7)) compared with its mean-field value at the thermodynamic limit \(\frac{1}{N}{\left[\sigma \right]}_{{{{{{{{\bf{J}}}}}}}}},\frac{1}{N}{\left[\frac{d\sigma }{dt}\right]}_{{{{{{{{\bf{J}}}}}}}},{{{{{{{\boldsymbol{\tau }}}}}}}}}\) (Eqs. (54) and (58)).

Full size image

To verify our theoretical predictions for the order parameters and steady-state entropy production, we compared them with the values computed from the sample trajectories by numerically simulating the kinetic SK models (see Supplementary Note 6 for the details). We constructed the kinetic Ising model with parameters Θi,u = 0 and randomly generated Jij with ΔJ = 0.5 and J0 = 1 while changing the inverse temperature β. We ran simulations of the model for t = 128 time steps and repeated the simulation 106 times at each β. We computed the mean activation rate \(\frac{1}{N}{\sum }_{i}{\left[\left\langle {s}_{i,u}\right\rangle \right]}_{{{{{{{{\bf{J}}}}}}}},{{{{{{{\boldsymbol{\tau }}}}}}}}}\), the average delayed self-correlations \(\frac{1}{N}{\sum }_{i}{\left[\left\langle {s}_{i,u}{s}_{i,u-1}\right\rangle \right]}_{{{{{{{{\bf{J}}}}}}}},{{{{{{{\boldsymbol{\tau }}}}}}}}}\), and the normalized entropy production and entropy production rates \(\frac{1}{N}{\left[\sigma \right]}_{{{{{{{{\bf{J}}}}}}}}},\frac{1}{N}{\left[\frac{d\sigma }{dt}\right]}_{{{{{{{{\bf{J}}}}}}}},{{{{{{{\boldsymbol{\tau }}}}}}}}}\) from trajectory and parameter sampling. We used the values at the last time step (u = t), where we confirmed that the statistics approached their steady-state values.

Figure 4 compares the theoretical order parameter m and q with the mean activation rate and delayed self-correlations computed from the simulated trajectories for system size N = 32, …, 1024 under synchronous (Fig. 4a, b) and asynchronous (Fig. 4d, e) updates. The simulated values approach the theoretical prediction as the size increases, albeit the convergence speed slows down near the critical temperature as it is expected. Similarly, we confirm in Fig. 4c, f that the entropy production from the sample trajectories for synchronous and asynchronous systems converges to the mean-field value at the thermodynamic limit as we increase the system size. Note that entropy production for synchronous updates differs from the entropy rate in the asynchronous update in continuous-time limit due to different values for the delayed correlation term q in Eqs. (54) and (58). These results corroborate our theoretical predictions that the steady-state entropy production peaks at the critical nonequilibrium phase transitions. We further verified by simulations with ΔJ = 1 that the steady-state entropy production increases when the significantly heterogeneous systems approach the quasi-deterministic regime (Supplementary Fig. S4).

Finally, knowing the order parameters of the system under the configurational average, we can investigate the structure of the patterns emerging from the dynamics of a sufficiently large but finite system under certain conditions. We calculate the probability Ω(n) of observing a state su+n again for the first time after n steps, starting from the same pattern su = su+n (Supplementary Note 7). For zero temperature (β → ∞) and synchronous updates, Ω(n) describes the probability of observing a periodic pattern of length n since transitions become deterministic and can result in periodic patterns. In general, the distribution of these patterns depends on the higher-order correlations between spins across time steps. However, we observe that for the disordered phase (m = 0) as well as in the deep ordered phase (m ~ 1), these correlations disappear for the configurational average in the asynchronous model or in the synchronous model for large n (see Supplementary Fig. S5). In these regions of the phase diagram, we can approximate the probability of observing a pattern n as \(\left[{\Omega }^{(n)}\right]\approx \exp \left[(1-n)\lambda \right]\) with \(\lambda={\left(\frac{1+m}{2}\right)}^{N}\) (see Supplementary Note 7). The expected length until a repeated pattern is observed is then

$$\mathop{\sum}\limits_{n}n\left[{\Omega }^{(n)}\right]\approx \frac{1}{\lambda }={\left(\frac{2}{1+m}\right)}^{N}.$$
(60)

At the disordered phase (m = 0), the average length of observed patterns exhibits a maximum value, growing exponentially at a rate of 2N. In contrast, when the system enters the ordered phase (m ~ 1), the growth rate decreases as m increases. In the limit m = 1, the system reaches a static equilibrium of average length 1, where the same pattern is repeated indefinitely. These results are consistent with expected dynamics under order-disorder phase transitions. Thus, bringing a quasi-deterministic system to a more stochastic regime by decreasing β to the critical value βc (with ΔJ smaller than ΔJc(∞)) increases the diversity of irreversible patterns and hence entropy production. However, further reduction of β makes the system more random (i.e., less irreversible transitions), leading to a decrease in entropy production. In contrast, the large entropy production found at the disordered phase at large β and ΔJ (wider than ΔJc(∞)) is caused by diverse oscillatory dynamics whose average pattern length is 2N as in the random dynamics (β = 0). Adding stochasticity to the dynamics by decreasing β in this regime reduces the entropy production monotonically.

The SK model with uniformly distributed external fields

Next, we apply non-zero external fields to the spins. For simplicity, we will consider the synchronous Ising system with unchanging fields Hi,u = Hi, assuming a uniform distribution \({H}_{i} \sim U\left(-\Delta H,\Delta H\right)\). Figure 5a, b shows the β − ΔJ phase diagram for the order parameters. With this change, we observe non-zero correlation q in the area where we previously saw the disordered states (m = 0 and q = 0, Fig. 2b). Figure 5c, d displays the order parameters as a function of the inverse temperature and ΔH, where we examine the effect of heterogeneity in the external fields while fixing the coupling variability, ΔJ = 0.2. The critical line of ΔHc(β) is obtained in this case as a solution of the following self-consistent equation (Supplementary Eq. (S3.9)):

$$\frac{\Delta H}{{J}_{0}}=\int\,{{{{{{{\rm{D}}}}}}}}z\,\tanh \left[\beta \left(\Delta H+\Delta Jz\right)\right].$$
(61)

Again, as studied in Supplementary Note 3, this critical phase transition corresponds to the mean-field universality class. Since the right-hand side term is less than or equal to 1 regardless of β and ΔJ, the phase transition occurs only when ΔH < J0 is satisfied. Intuitively, there is a competition between the dispersion induced by the field diversities ΔH and the cohesion induced by the mean coupling strength J0. The ordered phase takes place only if J0 counteracts the dispersion induced by the heterogeneity of external fields. More precisely, the critical ΔHc(β) at the low temperature limit (β → ∞) is obtained by solving \(\Delta H/{J}_{0}=\int\,{{{{{{{\rm{D}}}}}}}}z\,{{{{{{{\rm{sign}}}}}}}}\left[\Delta H+\Delta Jz\right]\). Here we have ΔHc(∞) = 1. We observe the phase transition by varying β if ΔH < ΔHc(∞), and no phase transition if ΔH > ΔHc(∞). Note that q increases monotonically with β even for ΔH > ΔHc(∞) when the distributed fields are introduced.

Fig. 5: Order parameters of the asymmetric SK model with heterogeneous fields.
figure 5

a, b The average magnetization m and average delayed self-coupling q are shown as a function of ΔJ and β. Fixed parameters are J0 = 1, ΔH = 0.5. The dashed line represents the critical line separating ordered and disordered phases. The horizontal dotted line represents the critical ΔJ at zero temperature (β → ∞). c, d The phase diagram of order parameters as a function of ΔH and β. Fixed parameters are J0 = 1 and ΔJ = 0.2 and variable ΔH. The dashed line is a critical ΔH at zero temperature.

Full size image

We now examine the conditional entropy, its reverse, and entropy production for the synchronous system with distributed fields, using the β vs. ΔH phase diagram. Similarly to the observation in the model without fields, the conditional entropy decreases with higher β (it becomes more deterministic processes, see Fig. 6a). The time-reversed conditional entropy also decreases with increasing β for all ΔH, indicating that the reverse process is more and more likely to happen regardless of ΔH (Fig. 6b). As seen previously, the time-reversed conditional entropy diminishes under the ferromagnetic-like states (ΔH < ΔHc(∞)). In contrast, we also observe the reduction of the time-reversed conditional entropy at higher β for ΔH > ΔHc(∞). Note that we observe increased correlations q at higher β for ΔH > ΔHc(∞) when we introduce non-zero external fields (see Fig. 5d), which results in a reduction of the reversed entropy similarly to ferromagnetic-like states. Both conditional and time-reversed conditional entropies decrease much slower along the critical line than in other regions, although with different magnitudes. As a result, we see the maximization of the entropy production around critical points more clearly than the β–ΔJ phase diagram (Fig. 6c, d). Finally, at the zero temperature limit (β → ∞), the entropy production peaks at ΔHc(∞) = 1 (Fig. 6e).

Fig. 6: Conditional entropies and entropy production of the asymmetric SK model with heterogeneous fields under synchronous updates.
figure 6

a The normalized conditional entropy Su∣u−1 (equivalent to the entropy rate under a steady state). b the normalized conditional entropy of the reverse dynamics \({S}_{u|u-1}^{r}\). c The normalized entropy production at a steady state. d Horizontal sections of the entropy production (ΔH = 0.2, 0.4, 0.6, 0.8, 1.0, and 1.2). It peaks at the critical line. e A vertical section of the entropy production at zero temperature (β → ∞). All plots are based on a model with fixed parameters J0 = 1, ΔJ = 0.2 and variable ΔH and β.

Full size image

Discussion

In this paper, we studied in detail the nonequilibrium thermodynamics of the kinetic, asymmetric SK model for both synchronous and asynchronous dynamics. As expected, the order parameters reveal that the model exhibits order-disorder nonequilibrium phase transitions analogous to the paramagnetic-ferromagnetic phase transitions in the equilibrium Ising model. There are, however, no phase transition akin to a spin-glass (which does not emerge due to coupling asymmetry, as previously reported for continuous-time asymmetric SK models25). In addition, we show that the steady-state entropy production is maximized near nonequilibrium phase transition points, being its first derivative discontinuous at these points (Fig. 3d). This result is similar to previously reported critical behavior of the entropy production caused by external stimulation or inertial dynamics (via self-coupling) in homogeneous systems with asynchronous updates by means of naive mean-field approximations or numerical simulations28,29,30. Nevertheless, our result is novel in that it provides the critical behavior of the entropy production caused by asymmetric, heterogeneous couplings using an exact analytical solution for such complex systems. In addition, the studied model displays a region with disordered oscillations in its phase diagram, where the entropy production takes even larger values than in the critical regime. This phase takes place for disordered systems with low entropy rates, i.e., the heterogeneous connections are strong enough to make the dynamics disordered but quasi-deterministic (Fig. 3c, top-right). In contrast, the entropy production does not increase when we increase the heterogeneity of external fields (Fig. 6c).

Taken together, our results indicate that the behavior of entropy production peaking at a critical point is more general than the simple mean-field, homogeneous models, therefore a non-smooth change of the steady-state entropy production (or entropy dissipated to an external reservoir) can be a useful indicator of a number of nonequilibrium phase transitions. At the same time, our results demonstrate that an increase in entropy production does not necessarily mean that the system is approaching a phase transition point. Instead, combining the order parameters, entropy rate, and entropy production yields a more precise picture of the complex systems and their phase transitions.

Typically, solutions of the symmetric (equilibrium) SK model involve the replica trick to calculate the configurational average of the logarithm of the partition function12. This method introduces an integer number of replicas of a system for averaging the disorder and then recovers the solution using a continuous number of replicas in the zero limit under the replica symmetry or replica-symmetry breaking ansatz. This treatment forces researchers to check the validity of solutions before reaching correct solutions16,53. As an alternative to the replica methods, the path integral methods have been widely used in analyzing the symmetric SK model54,55. However, for partially- or fully-symmetric SK models, the path integral method does not give a definite analytical solution but needs to be computed with Monte Carlo approaches56. Fortunately, the path integral method derives an exact analytical solution for the case of the fully asymmetric nonequilibrium SK model25, which we extended to cover synchronous and asynchronous updates, and theoretically underpinned their nonequilibrium properties by deriving the exact solution of the steady-state entropy production and entropy rates of the system.

Nonequilibrium properties of biological and adaptive systems have received the attention of neuroscience and biological science communities. For example, increased entropy production in macroscopic neural activity was suggested as a signature of physically and cognitively demanding tasks4, conscious activity5,6 or neuropsychiatric diseases like schizophrenia, bipolar disorders, and ADHD7. While it is not easy to contrast their findings based on the coarse-grained analysis of ECoG or fMRI data with the present results, our precise characterization of the entropy production of the prototypical system sheds light on what kind of behaviors we might expect from these complicated systems. Most importantly, our results indicate two scenarios to increase entropy production by controlling the connection heterogeneity (ΔJ) and neuron’s nonlinearity (β). These global changes in the model parameters can be achieved in the brain as gain modulation often mediated by neuromodulators57. One scenario to increase entropy production is that the system approaches a critical state as seen in the low ΔJ in Fig. 2 or Fig. 5. The other scenario is to make the system more heterogeneous and sensitive by increasing ΔJ and β. A significant difference is that the former process maintains stochastic nature while the latter yields quasi-deterministic disorder, as indicated by the high or reduced entropy rate. Therefore, the results suggest that it is crucial to investigate the multiple possibilities of nonequilibrium states to underpin the unconscious (sleep or anesthesia), awake, and engaged states more precisely.

Neuroscientists have often discussed the role of temporal patterns in spiking activities of neurons in computation or in memory consolidation and retrieval. One central topic of the debate is whether neurons, e.g., cortical or hippocampal ones, exhibit precise sequential patterns in a repeated manner58,59,60,61. Such precise sequences should result in a large entropy production similar to the low-temperature disordered phase of the kinetic Ising system. Alternatively, one may explain a broad range of irreversible temporal patterns, including avalanche dynamics62,63, by the dynamics near the nonequilibrium phase transitions without the precise sequential structure. As suggested by Eq. (60), the same network that shows simple periodic patterns at zero temperature can retrieve the diverse patterns yielding large entropy production by being poised near the critical phase transition point. The current study highlights the need to dissociate the two scenarios, characterizing different temporally irreversible spiking patterns to understand the distinct roles in neural computation using multiple thermodynamic quantities.

Finally, our analytical solutions offer a benchmark for – the aforementioned and other – methods for estimating thermodynamic quantities. For example, characterizing entropy production from brain imaging data requires methods for coarse-graining the phase diagram4,5. The kinetic SK model can serve as a test bench for such methods as it is an analytically tractable system with a well-known phase diagram. Moreover, we can use them to examine both established and novel mean-field theories in estimating the thermodynamic properties of large-scale systems. For example, one can directly fit the Ising model to neuronal spiking data using mean-field methods for finite-size networks, from which one can estimate various thermodynamic quantities of the system31,64. Accurately estimating these quantities in large networks gives deeper insights into the nonlinear computations of cortical circuitries. The exact solutions provided here serve to evaluate the accuracy of these approximation methods applied to large-scale networks and provide a benchmark of the thermodynamics quantities of infinitely large networks.