Stationary solutions of the master equation for single and multi-intermediate autocatalytic chemical systems

Abstract

In this work, we test a hypothesized form for the stationary solution Ps(X,Y) of the stochastic master equation for a reacting chemical system with two reactive intermediates X and Y, and multiple steady states. Thermodynamic analyses and the exact results for nonautocatalytic or equilibrating systems suggest an approximation of the form Pas(X,Y)=𝒩 exp(−φ/kT), where the function φ is a line integral of a differential ‘‘excess’’ work Fφ, which depends on species-specific affinities. The differential Fφ is inexact. In a preceding paper, we have given an analytic argument for the use of the deterministic kinetic trajectory, connecting (X,Y) to the steady state (Xs,Ys) as the path of integration for Fφ. Here, we show that use of the deterministic trajectories leads to a potential φdet which is continuous across the separatrix between the domains of attraction of the two stable steady states in the model studied. We compare the approximate form of Ps(X,Y) thus generated with numerical solutions of the time-dependent master equation in the limit of attainment of a stationary distribution. Because the time required for convergence to the stationary distribution scales as eN with the particle number N in cases with two stable steady states, the numerical work is limited to systems with 𝒪(10–100) X and Y particles. System size affects the accuracy of the approximation. To isolate system-size effects, we compare numerical solutions and the corresponding approximations to Ps(X) for two single-intermediate master equations, since the approximation becomes exact in the limit of large particle number for such equations. Based on these comparisons, for the systems with two intermediates, the agreement between the approximation and the numerical solutions is reasonable. The agreement improves as the number of particles increases in those test cases where it has thus far been possible to vary the system size over an order of magnitude. The results obtained by integrating along deterministic trajectories are better than those from straight-line or line-segment paths. The numerical work on small, single-variable systems with two stable steady states leads to two new observations: (i) the relative heights of the steady state peaks in the exact stationary distribution may invert as the system size increases and (ii) an approximation used commonly for particle counting may give results inconsistent with the exact stationary distribution when the particle number is small, while an alternative approximation improves the agreement.