Flows, scaling, and the control of moment hierarchies for stochastic chemical reaction networks

Stochastic chemical reaction networks (CRNs) are complex systems that combine the features of concurrent transformation of multiple variables in each elementary reaction event and nonlinear relations between states and their rates of change. Most general results concerning CRNs are limited to restricted cases where a topological characteristic known as deficiency takes a value 0 or 1, implying uniqueness and positivity of steady states and surprising, low-information forms for their associated probability distributions. Here we derive equations of motion for fluctuation moments at all orders for stochastic CRNs at general deficiency. We show, for the standard base case of proportional sampling without replacement (which underlies the mass-action rate law), that the generator of the stochastic process acts on the hierarchy of factorial moments with a finite representation. Whereas simulation of high-order moments for many-particle systems is costly, this representation reduces the solution of moment hierarchies to a complexity comparable to solving a heat equation. At steady states, moment hierarchies for finite CRNs interpolate between low-order and high-order scaling regimes, which may be approximated separately by distributions similar to those for deficiency-zero networks and connected through matched asymptotic expansions. In CRNs with multiple stable or metastable steady states, boundedness of high-order moments provides the starting condition for recursive solution downward to low-order moments, reversing the order usually used to solve moment hierarchies. A basis for a subset of network flows defined by having the same mean-regressing property as the flows in deficiency-zero networks gives the leading contribution to low-order moments in CRNs at general deficiency, in a $1/n$ expansion in large particle numbers. Our results give a physical picture of the different informational roles of mean-regressing and non-mean-regressing flows and clarify the dynamical meaning of deficiency not only for first-moment conditions but for all orders in fluctuations.

Figure 1

Graphical representation for a minimal CRN model and the first example in the sequence that will be developed in Sec. 6: one chemical species, two complexes, a single linkage class, and deficiency zero. These terms are defined in the remainder of the section.

Figure 2

Additional state added relative to the model of Fig. 1. This model is not weakly reversible.

Figure 3

Variant one-species model in which the complex graph is weakly reversible. The steady-state concentrations are the same as those in Fig. 2, which appear as a regular limit at $\epsilon \to 0$. Because $\delta =1$, the distribution at the steady state is no longer Poissonian.

Figure 4

Asymptotic expansions for moments and moment ratios for the model of Fig. 3. (a) Asymptotic expansion for $ln{\widehat{\mathrm{\Phi }}}_k$ descaled with $\xi =\sqrt{\alpha /\beta }$ under the recursion relation (68), upward from $k=1$. Five traces are generated by starting with $\langle N\rangle \equiv \langle N^3\rangle$ fixed and $\langle N^2\rangle$ values spaced by $1\times {10}^{-10}$ around the stable value. The group of trajectories become uncontrollably divergent by $k=25$. Asterisks are evaluations of the corresponding moments from a Gillespie simulation and the value of $\langle N\rangle$ was used to supply the unspecified normalization $c_0\approx 2.67$ in Eq. (68) for the recursion series. (b) and (c), which differ only in the plotted range, show ratios ${\widehat{\mathrm{\Phi }}}_{k+1}/{\widehat{\mathrm{\Phi }}}_k$ (for which $c_0$ appears only in the lowest term $\langle N\rangle /1$). Curves computed by recursion downward from $k=200$ with the starting approximation (69), shown as the red dash-dotted curve in (b). Gillespie simulation results overlaid as asterisks. The Poisson expectation $\xi$ from MFT and the large-$k$ asymptotic limit $\epsilon /\beta$ are shown as green dashed lines for reference.

Figure 5

This network, for some rate constants, can have two nonequilibrium steady states in the mean-field approximation. Even when this is the case, however, because the particle number is finite, the stochastic system always has only one ergodically sampled long-term steady distribution. Some subgraphs are common with Fig. 3, but the complex graph has two linkage classes, again giving $\delta =1$.

Figure 6

Solution $ln{\widehat{\mathrm{\Phi }}}_k$ to the recursion relations (33), descaled with ${\xi }^{\left(3\right)}=9$, extended downward from $k=200$ using the asymptotic approximation (77). The Poisson basis elements for ${\xi }^{\left(1\right)}=1$ and ${\xi }^{\left(3\right)}=9$ corresponding to the stable solutions in MFT are shown, respectively, as red dash-dotted and green dotted lines, for reference. The values ${\widehat{\mathrm{\Phi }}}_k$ for $k\in 0,...,3$ match an expansion (54) with the coefficients (79) and ${\langle {\mathrm{\Psi }}_Y\left(N\right)\rangle }^{\prime }\equiv 0$. Symbols are from a direct Gillespie simulation.

Figure 7

Effective potential $\mathrm{\Xi }\left(\overline{n}\right)$ from Eq. (81).

Figure 8

(a) Time series for the particle number $N$ in the four-state model of Fig. 5. Dotted lines label the minima of the effective potential (81) in the stationary-point expansion. (b) Histogram of the stationary distribution for $N$ from the simulations, showing a monotonic decrease and a shoulder with a mode around $N\approx 15$.

Figure 9

This CRN uses two species in a cross-catalytic configuration to produce the same potential for bistability that the network of Fig. 5 produces through one-species autocatalysis. In this CRN $\delta =2$ and the scaling behavior of each species separately is similar in many respects to one-species scaling of the network from Fig. 5.

Figure 10

Graphic representation of the action of the generator $\widehat{\mathrm{\Lambda }}$ given by Eq. (51) for the CRN of Fig. 9. The grid represents adjacent values of the index pair $(k_a,k_b)$, $(k_a\pm 1,k_b)$, $(k_a,k_b\pm 1)$, and $(k_a\pm 1,k_b\pm 1)$, as indicated by axis labels in the first panel. The closed circle in the first panel represents multiplication by the moment at $(k_a,k_b)$. Solid arrows indicate multiplication of moments by the projection vector $\left[\text{1}-1\right]$, acting on the moments at the tip and tail of the arrow. The labels on the lines indicate the function of parameters and $k$ values that multiplies each such projector. Thus the first panel is simply the term $-\eta (k_a+k_b){\widehat{\mathrm{\Phi }}}_{(k_a,k_b)}$, etc. The first and second panels, and the lower right pair of terms in the third panel, come from the terms at order $j_a+j_b=1$ in Eq. (51), given in Eq. (87). The terms running antidiagonally through the center in the third panel come from the terms at order $j_a+j_b=2$ and the terms in the upper left corner of the third panel come from the terms at order $j_a+j_b=3$, shown in Eq. (89). The large-$k$ asymptotic Poisson $\widehat{\mathrm{\Phi }}\to \mathcal{N}\underline{1}$ is annihilated identically by both panels with arrows.

Figure 11

Comparison of Gillespie simulations (circles) for the CRN model of Fig. 9 to the leading solution ${\phi }_{\kappa }^{\left(0\right)}$ (solid curve and crosses) from Eq. (94). (a) Ratios ${\widehat{\mathrm{\Phi }}}_{(k_a,k_b)}/{\widehat{\mathrm{\Phi }}}_{(k_a-1,k_b)}={\widehat{\mathrm{\Phi }}}_{(k_a,k_b)}/{\widehat{\mathrm{\Phi }}}_{(k_a,k_b-1)}$ for $k_a=k_b$ (black $\times$ and adjacent circles) and ${\widehat{\mathrm{\Phi }}}_{(k_a-1,k_b)}/{\widehat{\mathrm{\Phi }}}_{(k_a-1,k_b-1)}={\widehat{\mathrm{\Phi }}}_{(k_a,k_b-1)}/{\widehat{\mathrm{\Phi }}}_{(k_a-1,k_b-1)}$ for $k_a=k_b$ (red $+$ and adjacent circles). (b) and (c) Comparisons between these adjacent-moment ratios in the two-species model and moment ratios ${\widehat{\mathrm{\Phi }}}_k/{\widehat{\mathrm{\Phi }}}_{k-1}$ in the one-species model of Fig. 5 (red dash-dotted curve and asterisks) that has an equivalent MFT solution. (b) Plot of $k_a+k_b$ directly against one-species $k$, showing that for the lowest moments the total particle number controls similar moment values in the two models. (c) Plot of $(k_a+k_b)/2$ (equivalently, $k_a$ or $k_b$) against one-species $k$, showing that the transition to scaling dominated by the autocatalytic reactions is governed independently by $k_a$ and $k_b$.

Figure 12

Ratios of moments ${\widehat{\mathrm{\Phi }}}_{(k_a,k_b)}/{\widehat{\mathrm{\Phi }}}_{(k_a-1,k_b-1)}$ plotted versus $k_a+k_b\equiv k$, along the diagonal $k_a=k_b$. The series (94) is truncated at five successive orders of approximation ${\phi }^{\left({\alpha }_{\max }\right)}$ for ${\alpha }_{\max }=0,...,4$ (color sequence blue, green, red, cyan, magenta). Symbols are the ratios obtained directly from sampled moments of a stationary Gillespie simulation, showing that the order ${\alpha }_{\max }=0$ provides a good approximation to the diagonal moments, as shown in Fig. 11. Divergence of higher-order terms, which occurs with opposite sign for odd versus even ${\alpha }_{\max }$, reflects instability of the asymptotic expansion downward from large $\kappa$.

Figure 13

Graph of the residual error measure $\partial {\widehat{\mathrm{\Phi }}}_{(k_a,k_b)}/\partial \tau$ across the contour $k_a+k_b=49$, with the series (94) truncated at five successive orders of approximation ${\phi }^{\left({\alpha }_{\max }\right)}$ for ${\alpha }_{\max }=0,...,4$ (color sequence blue, green, red, cyan, magenta). Under perfect cancellation, the residual error at each order ${\alpha }_{\max }$ would scale as $q^{2({\alpha }_{\max }+1)}$. This scaling is very closely approximated for ${\alpha }_{\max }=1$ and degrades due to imperfect control of asymptotic expansions at higher orders, though the approximate behavior is attained. Crossing of the error curves suggests a finite radius of convergence in $q^2$ of the series (94), for $\left|q\right|/\kappa \sim 0.38$. The asymptotic expansion remains this good or better for all larger $\kappa$.

Figure 14

Linear color map of the error function $\partial {\widehat{\mathrm{\Phi }}}_{(k_a,k_b)}/\partial \tau$ for ${\alpha }_{\max }=5$, evaluated as the left-hand side of Eq. (95), which would equal zero for a stationary distribution. The deviation from zero is raised to the 0.1 power to produce a linear cross section if the true residuals scale as $q^{10}$.