Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the Schlögl model revisited

Abstract

Schlögl's model is the canonical example of a chemical reaction system that exhibits bistability. Because the biological examples of bistability and switching behaviour are increasingly numerous, this paper presents an integrated deterministic, stochastic and thermodynamic analysis of the model. After a brief review of the deterministic and stochastic modelling frameworks, the concepts of chemical and mathematical detailed balances are discussed and non-equilibrium conditions are shown to be necessary for bistability. Thermodynamic quantities such as the flux, chemical potential and entropy production rate are defined and compared across the two models. In the bistable region, the stochastic model exhibits an exchange of the global stability between the two stable states under changes in the pump parameters and volume size. The stochastic entropy production rate shows a sharp transition that mirrors this exchange. A new hybrid model that includes continuous diffusion and discrete jumps is suggested to deal with the multiscale dynamics of the bistable system. Accurate approximations of the exponentially small eigenvalue associated with the time scale of this switching and the full time-dependent solution are calculated using Matlab. A breakdown of previously known asymptotic approximations on small volume scales is observed through comparison with these and Monte Carlo results. Finally, in the appendix section is an illustration of how the diffusion approximation of the chemical master equation can fail to represent correctly the mesoscopically interesting steady-state behaviour of the system.

Figure 1

Plots of dx/dt versus x for the cases of chemical equilibrium (dashed curve) and bistability (solid curve) in the ODE model. The reaction rates are k1=3, k2=0.6, k3=0.25 and k4=2.95 for both plots. The pump parameters are a=0.5, b=29.5 for chemical equilibrium and a=1, b=1 for the bistable case.

Figure 2

Regions of monostability and bistability in the ab-plane with reaction rates fixed at the values k1=3, k2=0.6, k3=0.25 and k4=2.95.

Figure 3

Plots of the steady state in the CME model as b changes, with reaction rates fixed at k1=3, k2=0.6, k3=0.25, k4=2.95 and a=1.

Figure 4

Plots of the steady state in the CME model as V changes (with the same values of k _i as in figure 3; a=1 and b=2). (a) a=1, b=2, V=20; (b) a=1, b=2, V=40 and (c) a=1, b=2, V=80.

Figure 5

The Markov chain representation of Schlögl's model. The cycle moves forward in the clockwise direction. The steady-state flux is equal to the forward and backward cycle fluxes, which are the same in the steady state.

Figure 6

The entropy production rate from the master equation (divided by the volume) matches the entropy production rate of the ‘more stable’ steady state in the ODE model.

Figure 7

The stochastic entropy production rate versus the barrier height between the left FCA and right FCA . The lines have been darkened when the FCA is more stable to illustrate that the entropy production rate decreases as the height of the more stable barrier increases. The parameter values used are k1=3, k2=0.6, k3=0.25, k4=2.95, b=0.5 and 0.9≤a≤2.4. The point at which the FCA exchange stability is a=1.105, is indicated by the intersection of the bold and dashed lines. Grey dashed line, when n ₋ is less stable; grey solid line, when n ₊ is less stable; black dashed line, when n ₋ is more stable; black solid line, when n ₊ is more stable.

Figure 8

(a–d) The eigenvectors v _0,1 and w _0,1 using parameters in equation (4.4) and V=40. The deterministic steady states are indicated by dots on the n axis. Note that v ₀ is plotted on a logarithm scale so that both peaks are visible.

Figure 9

Relative residuals in the first four eigenvalue/eigenvector pairs, plotted on a logarithm scale (solid curve, λ ₀; dashed curve, λ ₁; dotted curve, λ ₂; dot-dashed curve, λ ₃).

Figure 10

Plot of the probability steady state (solid curve), calculated using detailed balance and as the eigenvector v ₀ in Matlab (dashed curve). Note that the drastic errors begin to occur in the Matlab's solution for values of p _ss less than machine epsilon, 10⁻¹⁶.

Figure 11

Comparison of Matlab's calculated λ ₁ with that obtained through Monte Carlo simulations. As the volume size increases, a λ ₁ decays exponentially. Solid line, Matlab, y=−0.22×−1.4; dot-dashed line, asymptotic, y=−0.27×−2.8; dashed line, Monte Carlo, y=−0.21×−1.8.

Figure 12

Convergence of the height differences to the asymptotic predictions Δu _± as V increases.