Entropy production selects nonequilibrium states in multistable systems

Far-from-equilibrium thermodynamics underpins the emergence of life, but how has been a long-outstanding puzzle. Best candidate theories based on the maximum entropy production principle could not be unequivocally proven, in part due to complicated physics, unintuitive stochastic thermodynamics, and the existence of alternative theories such as the minimum entropy production principle. Here, we use a simple, analytically solvable, one-dimensional bistable chemical system to demonstrate the validity of the maximum entropy production principle. To generalize to multistable stochastic system, we use the stochastic least-action principle to derive the entropy production and its role in the stability of nonequilibrium steady states. This shows that in a multistable system, all else being equal, the steady state with the highest entropy production is favored, with a number of implications for the evolution of biological, physical, and geological systems.

In this Supporting Information we derive equations used in the main text, and give additional explanations and results.
Since this equation is valid for all X, the detailed balancing equation follows W + (X)P (X) = W − (X + 1)P (X + 1).
Hence, the 1D nonequilibrium problem can be mapped onto an effective equilibrium problem, which can be solved by iteration for the steady-state probability distribution and effective potential. Specifically, the steady-state probability distribution is and using concentrations instead of copy numbers, i.e. x = 1/Ω ∞ X=0 XP (X) and scaling p(x) = ΩP (X), this leads to Here, the stochastic potential defined in the large-volume limit is given by with scaling W ± (X) = Ωγ ± (x) and the definitions γ + = w +1 + w −2 and γ − = w −1 + w +2 .
The rates are now given by where a = A/Ω and b = B/Ω indicate concentrations. The remaining integral in Eq. S10 can be calculated analytically, producing [3,4] Φ Fig. S1, along with other potentials used in the main text). The volume-independent prefactor in Eq. S9 is given by with Z a normalization constant. Stochastic potential Φ has indeed minima at the stable steady states of the deterministic solution and a local maximum at the unstable deterministic steady state since dΦ dx = ln γ − γ + (S14) Parameters: k +1 a = 0.5, k −1 = 3, k +2 = k −2 = 1, b = 4, and Ω = 10.
is zero for γ + = γ − and is > 0 (< 0) for stable (unstable) steady states with ∆γ = dγ − /dx − dγ + /dx. The sign of the second order derivative of Φ can easily be understood when considering standard linear stability analysis. Linearizing the ordinary differential equation of the Schlögl model (Eq. 14 in the main text) with x = x 0 + δx around steady state x 0 produces with perturbations decaying exponentially for ∆γ > 0.

II. TRANSITION RATES
Next, we consider the transition rates between the stable low (x * 1 ) and high (x * 2 ) steady states. Using a modified Fokker-Planck equation, known to correctly describe transport laws at large volume, the mean first-passage time is estimated from the splitting of the degenerate smallest eigenvalue [3,5]. The transition rate for switching from the low to the high state via the unstable state (x 0 ) can then be calculated via (S17b) Similarly, the transition rate for switching from the high (x * 2 ) to the low (x * 1 ) state is given by Their ratio is given by the simpler expression (S19) Rates scale with system size and become exponentially small with increasing Ω. However, unlike the conventional Fokker-Planck equation, these rates depend on exact stochastic potential Φ(x).

III. DERIVING THE SCHLÖGL MODEL FROM GENERIC BISTABLE MODEL
The Schlögl model based on the ordinary differential equation can be derived from the generic bistable model often used in similar form in biology [6], ecology [7], and climate [8,10] research. For example, in gene regulation x is the expression level of a protein (e.g. permease LacY in the bistable lactose utilization system [6]). In this case, Eq. S21 is interpreted as follows: the first term on the right-hand side represents basal expression, the second term is self-induced autoactivation by positive feedback, and the third term describes protein degradation. In another example from ecology [7], x represents the fraction of the lake surface covered by charophyte vegetation. Dramatic state shifts can occur through the sudden loss of transparency and vegetation observed in shallow lakes subject to human-induced eutrophication.
Returning to Eq. S21, performing a Taylor expansion of the denominator of the second term for k +2 x/(k −2 b) < < 1 to first order, we readily obtain the Schlögl model Eq. S20.
However, in the Schlögl model, x can take any non-negative value, i.e. it is not restricted to small values.

IV. DERIVATION OF ENTROPY PRODUCTION FORMULA
In the main text we used the entropy production formula with k B the Boltzmann's contant (omitted in the main text). In Eq. S22, the first term represents the net flux while the second term describes the ratio of the affinity (or 'thermodynamic force') and temperature T given by the ratio of the chemical potential difference of the reaction and T , i.e. ∆µ/T . In equilibrium both the net flux and the chemical potential difference are zero. Here, we derive forumla Eq. S22 following [9], using the general chemical where α, β, γ, and δ are stoichiometric coefficients, and A, B, C, and D are molecular species with concentrations a, b, c, and d. The rate constants for forward and backward rates are k + and k − , respectively. Assuming the law of mass action, the forward and reverse rates are respectively given by R + = k + a α b β and R − = k − c γ d δ , so that the net flux (or the reaction velocity) is given by What is the chemical potential difference in terms of the reaction rates? In general, the chemical potential difference (difference in Gibbs free energy) is given by where the chemical potential of some species X can be expressed as where µ 0 X is the chemical potential at a reference condition and x is the concentration of X (where assume an ideal solution with activity coefficient equal one). With this in mind, the chemical potential difference Eq. S24 can be written as The first term on the right-hand side of Eq. S26 stems from the chemical potential difference at the reference condition ∆µ is the equilibrium constant, only depending on temperature T but not on the chemical concentrations (or, more precisely, activities). The second term originates from the concentrationdependent log-terms of the chemical potentials. Upon further simplification, we finally which proves Eq. S22.

V. INTERCHANGEABLE POTENTIALS
In Eqs. 24a-c of the main text, the gradients of the ODE potentials appear. Here, we derive the corresponding Fokker-Planck (FP) potential. The FP equation is another approximation of the master equation and constitutes a partial differential equation for the probability density function p(a, b, x), containing contributions from the rate of change of the probability distribution due to currents (advection) and diffusion (noise) [11]. From the solution, the gradient of the potential would follow, e.g. Φ FP (a) = ∂Φ FP (a, b, x)/∂a.
However, for sufficiently small noise at steady state, the three species A, B, and X can be considered to fluctuate independently, and so we can derive the potential for each species separately. For instance, for species A the FP equation is (now without flux F ) where we assume that b and x are fixed at certain values (such as the steady-state values b * and x * , respectively). The steady-state solution is given by The gradient of the potential for species A is which contrasts the exact large-Ω solution of the master equation Φ (a) = ln(w +1 /w −1 ) (in analogy to Eq. 12 of the main text). A similar potential can be obtained for species B. In Eq. 29 of the main text, we use first use Eq. S30 and subsequently replace the FP potential by the exact potential from the master equation for large Ω.
As outlined in the main text, we make the assumption that we can replace time averages by ensemble averages. In the main text we used ẋΦ ODE (x) = 0 for the 1D Schlögl model at steady state, where Φ ODE (x) acts as a force from a conservative potential. While at steady stateẋ = 0 for the deterministic ordinary differential equation,ẋ is not well defined in the Langevin equation due to the noise (jumps). Heuristically, we can argue for the ensemble but here we will present a formal proof. Using Eq. 24a of the main text, we can proceed, now with an overall minus sign, with ξ x (t) white noise multiplied by amplitude √ g * x (i.e. constant effective temperature). Next, we need to calculate − √ g * x ξ x (t)Φ ODE (x) = η(t)G(x) , where we redefined noise η(t) = √ g * x ξ x (t) and flux G(x) = −Φ ODE (x). This calculation requires specification of the rule of discretization. Following [12] we Taylor expand the left-hand side of Using Stratonovich (midpoint) discretization for the second term on the right-hand side of with η(t) = 0 and η(t)η(t ) = (g * x ) 2 δ t,t , we obtain for Eq. S33 Hence, the quantity of interest from below Eq. S32b becomes Note, since G(x(t)) and η(t) are independent, η(t)G (x(t))G(x(t)) = η(t) G (x(t))G(x(t)) = 0 (and similarly η(t)G (x(t))G(x(t + τ )) = 0) and η(t)G (x(t))η(t + τ ) = η(t)η(t + τ ) G (x(t)) = 0 due to white-noise property.
Put into Eq. S32b this results in due to an equipartition-type theorem for this effective equilibrium system when close to the potential energy minimum. To illustrate this, assume for simplicity a harmonic potential Φ ODE (x) = Kx 2 /2. Now, we obtain relation K x 2 = T eff /2, where T eff = (g * x ) 2 is an effective temperature. Hence, species X appears to effectively be at equilibrium without producing entropy.
which now contains imposed flux F . We further obtain √ g * a ξ a Φ ODE (a) = 1/2 (g * a ) 2 Φ ODE (a) , which we used in the main text. Now there is entropy production due to ȧΦ ODE (a) = 0 from flux F , and similarly for species B.

VII. ENTROPY PRODUCTION DOES NOT DETERMINE PROBABILITY OF TRAJECTORY ALONE
Here we use simulations to show that trajectories are not solely selected thermodynamically based on maximal entropy production in line with the main-text result that the classical action needs to minimized as well. To investigate this, we use the entropy production defined by the log-likelihood ratio ∆S Γ = ln(P Γ /P −Γ ) with P Γ and P −Γ the probabilities of the forward and backwards (time-reversed) trajectories Γ and −Γ, respectively. Specifically, we have [2] for n events during time t (see also Eq. 8 in the main text). Based on a trajectory, e.g.
as obtained from exact Gillespie simulations for the Schlögl model, the numerator in the logarithm represents the product of all the transition rates, while the denominator contains the product of the transition rates in the time-reversed direction. As can be seen in Fig.   S2, the slope fluctuates between low and high entropy production rates in line with maintext Fig. 2a,d. Ensemble averaging over many trajectories produces the average entropy production ∆S Γ [2]. If plotted for increasingly longer trajectory durations t, then the slope represents the average entropy production rate (up to a potential transient at the beginning if not initiated at steady state) [2,[13][14][15]. At equilibrium both numerator and denominator are equal and ∆S Γ vanishes as expected.
simulations for different b mod values and then evaluated ∆S Γ using a specific b value (here b = 4 inside the bistable regime). This way we are able to use perturbed trajectories, different from the correct trajectories that solve the master equation for a certain set of parameters. Fig. S3 shows that the ensemble-averaged entropy production ∆S Γ (t) is an increasing function of b mod . In particular, the correct b mod = b does not correspond to a maximum in entropy production. This can be understood by considering the log-likelihood ratio in Eq.
S39. The entropy production is large whenever the probability of the forward trajectory is much larger than the probability of the backward trajectory. However, this statement is independent of how likely the forward trajectory is in the first place, i.e. whether it is a solution of the master equation. Hence, the maximum entropy production principle cannot be the sole determinant for selecting the most likely trajectory in line with the main text result.