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

TL;DR

This paper develops a comprehensive mathematical framework to analyze the behavior of stochastic chemical reaction networks at all deficiency levels, revealing new insights into their moment hierarchies, scaling regimes, and steady-state properties.

Contribution

It derives equations for all-order moments in CRNs with arbitrary deficiency, extending previous results limited to deficiency 0 or 1, and introduces a finite representation of the generator acting on factorial moments.

Findings

Moment hierarchies interpolate between scaling regimes at steady state.

Boundedness of high-order moments aids in solving for low-order moments in multi-stable systems.

A $1/n$-expansion characterizes the contribution of certain network flows to moments.

Abstract

Stochastic chemical reaction networks (CRNs) are complex systems which 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 value 0 or 1. Here we derive equations of motion for fluctuation moments at all orders for stochastic CRNs at general deficiency. We show, for the case of 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 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-0 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-0 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.

Figures

Equations

$$\mathrm{A}$$

$$\mathrm{A}\stackrel{\alpha }{\rightharpoonup }2\mathrm{A}.$$

$$\mathbb{A}$$

$$=\sum _i{w}_i\sum _j\left(k_{ji}{w}_j^T-k_{ij}{w}_i^T\right)$$

$${\Psi }_Y^i\left(\mathrm{n}\right)$$

$$\frac{\partial \rho }{\partial \tau }$$

$$\text{mbox}shorthandfor\frac{\partial {\rho }_{\mathrm{n}}}{\partial \tau }$$

$$\varphi \left(z\right)\equiv \sum _{\mathrm{n}}\left(\prod _p{z}_p^{{\mathrm{n}}_p}\right){\rho }_{\mathrm{n}}.$$

$$\frac{\partial \varphi }{\partial \tau }$$

$$\text{mbox}shorthandfor\frac{\partial }{\partial \tau }\varphi \left(z\right)$$

$$⟨O\left(\mathrm{n}\right)⟩\equiv \sum _{\mathrm{n}}O\left(\mathrm{n}\right){\rho }_{\mathrm{n}}$$

$$n\equiv ⟨\mathrm{n}⟩=\sum _{\mathrm{n}}\mathrm{n}{\rho }_{\mathrm{n}},$$

$$n={\frac{\partial }{\partial z}\log \varphi \left(z\right)|}_{z\equiv 1},$$

$${\mathrm{n}}_p^{\underline{k_p}}$$

$$\equiv 0$$

$${\Phi }_k\equiv ⟨\prod _p{\mathrm{n}}_p^{\underline{k_p}}⟩.$$

$$\frac{\partial \Phi }{\partial \tau }$$

$$\text{mbox}shorthandfor\frac{\partial {\Phi }_k}{\partial \tau }$$

$$z_p$$

$$\left[a_p,a_q^{†}\right]={\delta }_{pq},$$

$$1$$

$$\prod _{p=1}^P{z}_p^{{\mathrm{n}}_p}\times 1$$

$$a_p^{†}{a}_p\left|\mathrm{n}\right)={\mathrm{n}}_p\left|\mathrm{n}\right).$$

$$a_p^k\left|\mathrm{n}\right)$$

$${a_p^{†}}^k{a}_p^k\left|\mathrm{n}\right)$$

$$\varphi \left(z\right)=\sum _{\mathrm{n}}\prod _{p=1}^P{z}_p^{{\mathrm{n}}_p}{\rho }_{\mathrm{n}}$$

$$\left(0\right|e^{{\sum }_p{a}_p}\left|\mathrm{n}\right)=1,\forall \mathrm{n},$$

$$\left(0\right|e^{{\sum }_p{a}_p}\left|\varphi \right)=\sum _{\mathrm{n}}{\rho }_{\mathrm{n}}=1.$$

$${\Psi }_Y^i\left(\mathrm{n}\right)$$

$${\psi }_Y^i\left(a\right)$$

$${{\psi }^{†}}_Y^i\left(a^{†}\right)\equiv \prod _p{a^{†}}_p^{y_p^i}$$

$${\psi }_Y^i\left(e^{\partial /\partial \mathrm{n}}\right)\equiv \prod _p{e}^{y_p^i\partial /\partial {\mathrm{n}}_p}=e^{{y^i}^T\partial /\partial \mathrm{n}};$$

$$\mathrm{T}$$

$$=\sum _{\left(i,j\right)}\left[{\psi }_Y^j\left(e^{-\partial /\partial n}\right)-{\psi }_Y^i\left(e^{-\partial /\partial n}\right)\right]k_{ji}{\psi }_Y^i\left(e^{\partial /\partial n}\right){\Psi }_Y^i\left(\mathrm{n}\right).$$

Full text

Keywords: Chemical Reaction Network, deficiency, stochastic processes, moment hierarchies, Anderson-Craciun-Kurtz theorem

I Introduction

Random walks of multiple independent particles on ordinary graphs are simple processes in several fundamental senses Lovasz:survey:93 ; Luo:spectral_embed:03 ; Lee:spectral:10 . Each step involves a change of a single degree of freedom, and for the base case of proportional sampling without replacement, the rate law for fluxes is linear in concentrations. Moreover, identification of topological properties of ordinary graphs which determine characteristics of random-walk dynamics can generally be carried out in polynomial time Lovasz:survey:93 ; Hordijk:autocat_alg:04 .

Just the opposite is true of the stochastic processes associated with Chemical Reaction Networks (CRNs). Each elementary reaction event (for a general network) can involve the concurrent conversion of multiple inputs into multiple outputs  Danos:rule_based_modeling:08 ; Harmer:info_carriers:10 , making the elementary network on which the reactions occur not an ordinary graph, but a directed multi-hypergraph  Andersen:comp_rules:13 ; Andersen:generic_strat:14 . For the base case of proportional sampling without replacement over reactants, the rate law becomes nonlinear, leading generically to possibilities for complex dynamics and multiple (stable or metastable) steady states Feinberg:notes:79 ; Gunawardena:CRN_for_bio:03 ; Baez:QTSM:17 . Moreover, identification of key topological properties such as shortest reaction sequences connecting inputs to outputs or complete sets of self-amplifying cycles, which affect the character of CRN dynamics, are known to be NP-hard problems for hypergraphs Berge:hypergraphs:73 ; Andersen:NP_autocat:12 .

As a consequence, although stochastic processes with the essential features of CRNs are ubiquitous in biochemistry Metzler:BC:03 , systems biology Palsson:systems_bio:06 , ecology Smith:pop:11 , and epidemiology Allen:stoch_process:03 , and are thus of large practical and theoretical interest, few results exist for these systems Feinberg:def_01:87 ; Craciun:multistability:05 ; Craciun:bistability:06 ; Ji:PhD:11 ; Joshi:Survey:15 compared with the large literature that exists for random walks on ordinary graphs Lovasz:survey:93 . In addition to systems which clearly have CRN structure because the underlying processes obey constraints of stoichiometry (the source of concurrency), the CRN framework is flexible enough to furnish a representation for systems of broad interest to non-equilibrium thermodynamics such as the zero-range process Evans:zero_range:05 , where topological characteristics that are known to lead to simple steady states for CRNs can be used to sieve for exactly solvable cases. Indeed, some of the key results that are known for stochastic CRNs Anderson:product_dist:10 were motivated in part by earlier work of Kelly on a related class of queuing networks Kelly:rev_stoc_networks:79 .

Partitioning off the simple sub-architecture in the

core of complex CRNs

A representation scheme for CRNs was made standard by the work of Feinberg Feinberg:notes:79 ; Feinberg:def_01:87 , which separates the structure of independent reaction events from the stoichiometric relations that determine their action in the chemical state space. We will emphasize that, well beyond the use made by its original authors, this decomposition defines the fundamental partition in CRN architecture between a sub-system with the same simplicity as the random walk on an ordinary graph, and the remainder of the CRN constraints responsible for concurrency, nonlinearity, and their resulting complexity. Identifying the subsystem isomorphic to a simple process is the key to decomposing the scaling behaviors in the moment hierarchies of CRNs, and to distinguishing the roles of different classes of flows in the dynamics and steady states.

The emphasis in the work of Feinberg, and in closely-related work by Horn and Jackson Horn:mass_action:72 , was the existence and uniqueness of solutions to the mass-action rate equations with strictly positive concentrations: a limited, deterministic, and static problem. We will show below that the Feinberg decomposition is even more useful in the analysis of the stochastic processes associated with CRNs, where in addition to exact results in deterministic limits, it can serve as a foundation for systematic approximation methods in the general case. A host of results follow, including novel representations of the generator of the stochastic process acting on the moment hierarchy, duality relations of the kind explored in Stochastic Thermodynamics Seifert:stoch_thermo_rev:12 , and scaling relations that suggest solution methods using matched asymptotic expansions.

An expanded role for deficiency

A new dimensional property termed deficiency was introduced by Feinberg Feinberg:notes:79 , which was central to both the results on existence and uniqueness of steady states, and to a limited but important corollary about the form of some of the distributions associated with such states. Deficiency can be computed as a topological index from the graphical structure associated with a CRN Feinberg:notes:79 ; Gunawardena:CRN_for_bio:03 ; Baez:QTSM:17 , but its importance comes from its meaning as a count of the dimensionality of chemical flows that can proceed in a steady state without being subject to mean-regression due to changes in chemical concentrations. When there are no such regression-free flows – when deficiency equals zero – the network is guaranteed to have unique strictly positive steady states at general parameters. A remarkable result due to Anderson, Craciun, and Kurtz (ACK) Anderson:product_dist:10 is that under the same conditions, the steady-state distributions have a simple factorial form under proportional sampling rules, and a class of related forms under more general rules (e. g., Michaelis kinetics), as long as they sample chemical species independently of one another.

The concept expressed by deficiency remains key to organizing the stochastic processes for general CRNs, even when their deficiency is nonzero, and the clues for why this should be so are already latent in the ACK theorem  Anderson:product_dist:10 . For the simple case of proportional sampling without replacement (which gives rise to mass-action kinetics, and which we will assume in the remainder of the article), the ACK solutions are either products of Poisson distributions, or hypersurfaces within such product-Poissons constrained by particle conservation laws (with no loss of generality for the claims below). In a Poisson distribution, all higher moments are universal functions of the mean value, so in a product-Poisson, the entire moment hierarchy is controlled by the set of first-moments. In a deficiency-zero network at steady state, the ACK theorem effectively states that first-moment values carry all the “information” in the distribution. When deficiency is nonzero, the mean-regressing flows are no longer the exclusive dynamical entities, but as we will show they remain the dominant entities governing low-order moments, in an asymptotic expansion where the small parameter is the order of any moment of the distribution relative to its mean particle number. The approach by which we will construct this result also shows how to extract other scaling regimes associated with the remaining flows in networks with non-zero deficiency, and the way these control complementary asymptotic expansions for high-order moments relative to mean particle numbers.

Coupling CRN theory to Doi operator algebras for

stochastic processes

Our approach to the stochastic processes associated with CRNs grows out of a set of linear-algebra methods due to Masao Doi Doi:SecQuant:76 ; Doi:RDQFT:76 , for the treatment of generating functions for general discrete-state stochastic processes. The Doi operator algebra provides a starting point for numerous solution methods,111One of the best-known of these is the coherent-state expansion of generating functionals due to Luca Peliti Peliti:PIBD:85 ; Peliti:AAZero:86 . It is particularly useful for semiclassical approximations and other stationary-point methods, which we mention but do not pursue in depth here. but its simplifying effect is particularly elucidating for CRNs. In the classical mass-action theorems of Feinberg, the rate equations for CRNs take an awkward and not-very-perspicuous form, in which stoichiometry is expressed asymmetrically in nonlinear rate laws and concurrency constraints on inputs and outputs. In the Doi representation of the full stochastic process, all formal asymmetry between inputs and outputs disappears. It can be seen that the asymmetry of the Feinberg problem reflects the way the particular projection operator to the first-moment equations of motion interacts with the formally-symmetric generator of all fluctuations. Within the framework of the Doi algebra, we formulate the projection operator that gives the equations of motion for arbitrary moments, and show how the Feinberg equations generalize to a new representation of the generator of the stochastic process acting on the moment hierarchy, which directly expresses the scaling influence of different network flows. The underlying symmetry of the Doi representation of the stochastic process generator remains, and can serve as a point of departure for the derivation of duality relations for stochastic CRNs, which we mention here but develop in a separate publication Smith:DP_duality:17 .

Organization of the presentation

The presentation is organized as follows:

In Sec. II we introduce the general concepts and notation associated first with Chemical Reaction Networks, then with general discrete-state stochastic processes, and finally for the particular forms of stochastic processes associated with CRNs. The culmination of this section is Eq. (27), the Liouville-form expression for the generator of the stochastic process of a CRN motivated by the Feinberg decomposition, which is the basis for all other results in the paper.

Sec. III introduces the factorial moments which are the natural observables for CRNs with simple mass-action rate laws, and shows that a dynamical equation for all orders of fluctuations is closed and has a finite-order generator in this set of moments. The culmination of this section is the representation of the generator in Eq. (57).

Sec. IV then shows how topological characteristics of the CRN are linked to dynamical properties of flows, and reviews the Feinberg deficiency-0 theorem and the associated Anderson-Craciun-Kurtz theorem. We introduce what we term the stoichiometric decomposition of the Liouville operator in Eq. (63), which separates two dynamically different classes of mean-regressing and non-mean-regressing flows.

Sec. V uses the representation of the generator on the lattice of factorial moments to show how different combinations of rate constants from a CRN govern scaling properties of moments in different ranges of the moment order. This section shows how matched asymptotic expansions can be used to solve for steady-state moment hierarchies recursively, and shows the (related) sense in which the deficiency-0-like subset of flows in the stoichiometric decomposition dominate low-order moments. Sec. VI then provides a sequence of worked examples of ascending complexity to introduce each of the concepts above and show its effects in a solution.

We have chosen to develop all main results in their general forms in Sections II – V, in the interest of economy and continuity of the argument, postponing examples to Sec. VI where they may be directly compared. The simplest case and the starting point involves one species and deficiency zero (VI.1), illustrating the way the ACK theorem is recovered in the Doi algebra with a one-line proof. We then introduce nonzero deficiency while keeping the same mass-action equations of motion (VI.2), to show how the subspace of flows resembling a zero-deficiency network can be extracted (but also why, in the general case, this cannot be expressed as a zero-deficiency sub-process of the full process), and how it controls the scaling of low-order moments. Next we hold deficiency fixed but change the network topology to one in which autocatalytic feedback produces multiple steady states in the mass-action approximation (VI.3). This case introduces the first non-trivial role for the asymptotic expansion and shows, counter-intuitively, how a criterion of boundedness for asymptotically high-order moments anchors the recursion downward to specify the low-order moments of the ergodic distribution over the two steady states. All effects up to this point are illustrated with single-species networks. At the end we introduce a two-species network in which cross-catalysis replaces the single-species autocatalysis (VI.4), yielding an equivalent bistable classical system if the species are not distinguished. This case demonstrates a non-trivial use of the stochastic process generator acting on the moment hierarchy, and shows how the factorability of the ACK theorem for multi-species distributions is lost at non-zero deficiency.

Readers who prefer to alternate general notation and instances are encouraged to browse the examples in Sec. VI in parallel with reading the formal development in the earlier sections. A more direct track to the main results on recursions in the moment hierarchy is also provided in Krishnamurthy:CRN_moments:17 .

II Chemical Reaction Networks and their associated stochastic

processes

The next three sub-sections review the standard concepts for CRNs and two aspects of stochastic-process algebras – representations of the generator and then the Doi operator formalism – and introduce the notation in which we will represent them in this paper. In Sec. II.4 these are then brought together to obtain the Liouville-operator representation of the generator for a stochastic process CRN that will be the basis for all further constructions.

II.1 The elements of a CRN

Our decomposition of CRNs follows that of Feinberg Feinberg:notes:79 , but we use a more complete graphic representation, illustrated for our simplest example process in Fig. 1. We will present CRNs using both this graphic form, and a corresponding reaction-scheme form such as

[TABLE]

The following list explains the fundamental division of Feinberg, Horn and Jackson that separates the so-called complex network222Here the stress is on the first syllable – 'com,plex network – to be distinguished from references to systems that are com'plex. from the stoichiometric relations that interpret its action in the state space of chemical species:

The chemical species: In examples we will denote explicit species names with single capital Roman letters such as ${\rm A}$. Where a set of species is indicated, we index them with subscripts $p\in 1,\ldots,P$, the number of distinct chemical species. In network diagrams a species is denoted with a filled dot:

Complex: a multiset of species, which is the input or output of a reaction. It may be written as a sum rather than a set. Examples might be ${\rm A}$, $2{\rm A}$ (or equivalently ${\rm A}+{\rm A}$). In network diagrams, a complex is denoted by an open circle with one or more (labeled) dashed line-stubs indicating the reaction(s) in which it participates and (labeled) solid line-stubs indicating the participating species: (In one-species or one-reaction networks, or after the network diagram has been assembled, we may suppress the labels to improve readability.) In set notation, we index complexes with subscripts $i$ or $j$.

(Directed) reaction: an ordered pair of complexes with an associated rate constant. For complexes indexed $i$ and $j$, respectively, the reaction from $i$ to $j$ would correspond to the ordered pair $\left(i,j\right)$, and the associated rate constant is denoted $k_{ji}$. (For simple models, rate constants may be given simplifying labels such as $\alpha$.)

An expression with one or more reaction is known as a reaction scheme, such as

[TABLE]

In network diagrams, a reaction is represented with a dashed arrow between the input and output complexes (optionally labeled with the rate constant):

Chemical Reaction Network (CRN): a collection of reactions. By default we adopt the finest-grained description in which all reactions are unidirectional; bi-directional reactions are indicated with pairs of directed arrows. In reaction schemata we may also condense notation for bidirectional reactions, as above in Scheme 1.

Graphically a CRN is a well-formed doubly-bipartite network (two kinds of nodes and two kinds of links), in which all reactions terminate in two complexes, and all line-stubs from a complex are filled by the appropriate reactions or links to chemical species333Strictly speaking, a CRN corresponds to a directed multi-hypergraph, in which the reactions correspond to directed hyperedges, and the input and output complexes are their vertex sets. The graphic depiction used here and elsewhere is a bipartite representation of the underlying hypergraph. The computational complexity of numerous search and optimization problems on hypergraphs, which are simple on ordinary graphs, is one consequence of the concurrency of inputs and outputs on a hyperedge., as in Fig. 1.

Adjacency/rate matrix: In Feinberg’s representation of CRNs Feinberg:notes:79 ; Feinberg:def_01:87 , only the complexes and the reactions are denoted explicitly, and they form an ordinary, directed graph. In a stochastic formulation, reactions can occur as independent events on the links, analogous to the steps in a simple random walk. Complexes are treated as if they have activities, and the set of rate constants map these activities to reaction rates. In this way, both stoichiometric constraints and the determination of activities of complexes are cordoned off as separate information from the rate- and connectivity-structure of the complex-network. The latter is given by an ordinary adjacency/rate matrix, identical in form to the graph Laplacian for a simple random walk. Following the notation made standard by Feinberg, we denote this matrix by $\mathbb{A}$.

Arranging complexes in a column vector indexed by $i$, let $w_{i}$ be the indicator function that is nonzero on complex $i$ only (so the $j$th component ${\left(w_{i}\right)}_{j}\equiv{\delta}_{ij}$), and $w_{i}^{T}$ its transpose which can act as a projection operator. Then the adjacency/rate matrix can be written as a sum of dyadics, arranged in a variety of ways. Two of these we call the reaction representation and the complex representation, written

[TABLE]

In expressions such as the canonical mass-action rate law – which will appear as Eq. (64) below, after its context has been properly introduced – the row vectors $w^{T}$ select the activities determining reaction rates, and the column vectors $w$ identify the net flux into and out of complexes. The reaction representation is often more intuitive for constructing generators of the stochastic process, while the complex representation, in accumulating all flows into or out of a complex, more directly reflects the cause of mean-regression that underlies the concept of deficiency.

Stoichiometric matrix: In order to connect complexes to species, treat ${\rm n}\equiv\left[{{\rm n}}_{p}\right]$ as a column vector, and introduce a matrix $Y$ with rows $y_{p}\equiv{\left[y_{p}^{i}\right]}^{T}$. $y_{p}^{i}$ are the stoichiometric coefficients indicating the number of instances of species $p$ in complex $i$. In graphs such as Fig. 1, this is the number of solid lines from the species node $p$ to the complex node $i$. The product $Y\mathbb{A}$ connects fluxes at complexes (the row index on $\mathbb{A}$) to fluxes at species (the row index on $Y\mathbb{A}$).

Interpretation of complex activities: Complexes are not the same as chemical species, and their (virtual) activities in the Feinberg complex-graph must be fixed in terms of the (actual) activities of the species. We will call any such dependence an interpretation of the complexes.444It is an interesting question, which we leave for other work, whether the random walk on the complex network can be represented formally in terms of a linear algebra of pseudo-particles, which are formal proxies for the products of operators representing real particles of the chemical species.

To give an interpretation of complex activities in terms of species activities that is convenient to use with the adjacency/rate matrix, let ${\Psi}_{Y}\equiv\left[{\Psi}_{Y}^{i}\right]$ be a column vector with components given by

[TABLE]

${\Psi}_{Y}$ defines the activity products corresponding to proportional sampling without replacement on the discrete indices ${{\rm n}}_{p}$. We will return below to the relation between discrete sampling and mass-action rate laws, after we have introduced notation and concepts for the stochastic process that governs mass action and all higher moments.

II.2 Representations of a stochastic process

A certain standard machinery underlies all discrete-state stochastic processes, including those associated with CRNs. Before deriving the particular forms for the network decomposition of Sec. II.1, we introduce general notation and constructions in this section.

The point we wish to emphasize is that, while the stochastic process for a CRN has a uniquely-defined generator, that generator may have many representations, depending on whether we solve the stochastic process for its probability density, or the moment-generating functional of that density, or the hierarchy of moments evaluated directly.

When the activities of complexes are defined by proportional sampling without replacement as in Eq. (4) (the simplest case, corresponding to ideal gases or solutions), the moments that appear naturally in all rate equations are what we term factorial moments, and it is for these that the equations of motion take the most compact form.

The essential components in the stochastic-process description are then the following:

Probability density function and transfer matrix: States of the CRN are indexed by values of the vector ${\rm n}$, and reactions are treated as instantaneous changes of state.555This level of coarse-graining in the description of reaction events is the standard assumption also in Stochastic Thermodynamics Seifert:stoch_thermo_rev:12 . Our starting point in describing the stochastic process is a probability density function ${\rho}_{{\rm n}}$ indexed on the values of ${\rm n}$.

The density $\rho$ evolves on a time coordinate $\tau$ under a master equation

[TABLE]

The matrix ${\rm T}\equiv\left[{\rm T}_{{\rm n}{{\rm n}}^{\prime}}\right]$ is called the transfer matrix, and is one representation of the generator of the stochastic process.

Moment-generating function and Liouville operator: The moment-generating function is formed from ${\rho}_{{\rm n}}$ with the introduction of a vector $z\equiv\left[z_{p}\right]$ of complex coefficients, as the Laplace transform

[TABLE]

The generating function evolves under a Liouville equation of the form

[TABLE]

$\mathcal{L}$ is called the Liouville operator. Its form is defined from the transfer matrix in Eq. (5), and it provides an alternative representation of the generator of the stochastic process, acting on Laplace transforms.

Expectations and moments and their time-dependence: The expectation of an arbitrary function $O\!\left({\rm n}\right)$ (for “Observable”) of the components of ${\rm n}$, in the background $\rho$, is denoted

[TABLE]

Since $\rho$ may be a continuous-valued quantity whereas ${\rm n}$ is discrete, we introduce a particular short-hand with math-Italic font for the first moment

[TABLE]

which may vary continuously if $\rho$ does. The mass-action rate equations are expressed entirely in terms of $n$. Depending on the form of the underlying distribution $\rho$ and whether higher-order correlations can be expressed as functions of $n$, the mass-action equations may be exact or they may involve a (generally-unregulated) approximation known as the mean-field approximation.

$n$ may also be obtained from the generating function as

[TABLE]

and expectations of more complex observables can be built up by acting on $\phi$ appropriately with higher-order derivatives in $z$.

The time-dependence of $n$, which is the object of classical first-moment equations or chemical rate equations, can be obtained by acting on either $\rho$ or $\phi$ with its corresponding generator, by Equations (5) or (7).

Factorial moments for CRNs: In many applications, the natural moments for which to study dynamics are either ordinary powers ${{\rm n}}^{k}$, obtained by acting on $\phi$ with higher derivatives in $z$, or cumulants vanKampen:Stoch_Proc:07 , obtained by acting on $\log\phi$ with higher derivatives in $z$. For CRNs with reaction rate laws corresponding to proportional sampling without replacement, however, the simplest dynamical relations are obtained for the factorial moments, for which we therefore introduce a specific notation. For a single component ${{\rm n}}_{p}$ and power $k_{p}$,

[TABLE]

(We adopt the convention in the second line of Eq. (11) for $k_{p}>{{\rm n}}_{p}$ because it allows a simplification in sum notations below.) We will use the form of Eq. (11) as a general formula for truncated factorials, allowing other integer-valued arguments such as stoichiometric coefficients to take the place of ${{\rm n}}_{p}$ as the argument.

Moment hierarchy and a new representation for the generator: For a vector $k\equiv\left[k_{p}\right]$ of powers, we introduce the factorial moment hierarchy indexed by $k$, as the expectation

[TABLE]

We show in Sec. III that $\Phi$ evolves in time as

[TABLE]

The matrix $\Lambda\equiv\left[{\Lambda}_{kk^{\prime}}\right]$ has finitely many nonzero entries determined by the stoichiometric coefficients, and its form may be derived from the Liouville operator $\mathcal{L}$. $\Lambda$ provides yet a third representation of the generator of the stochastic process, which is particularly well-suited to the study of CRNs because it exposes different scaling regimes controlled by combinations of the rate parameters corresponding to different flows. We develop the implications of scaling in Sec. V.

II.3 The Doi operator algebra for a discrete-state stochastic

process

For most of its purposes as a generating function, it is not necessary that $\phi\!\left(z\right)$ be an analytic function of a vector $z$ of complex-valued arguments. Often only the formal power series in $z$, and its algebra with the derivative $\partial/\partial z$, is required.

The abstraction of the linear algebra of generating functions in terms of formal raising and lowering operators follows a procedure due to Masao Doi Doi:SecQuant:76 ; Doi:RDQFT:76 . We have elaborated the details of the mapping and its interpretation extensively elsewhere Smith:Signaling:11 ; Smith:LDP_SEA:11 ; Smith:evo_games:15 , and here we only summarize the notation, which by now is standard.666Indeed, most treatments open directly with the Doi algebra Mattis:RDQFT:98 ; Cardy:FTNEqSM:99 . We have used the two-step introduction by way of conventional analytic generating functions because it clarifies the meaning of some terms in the Doi algebra that can be obscure when presented without introduction.

The Doi algebra denotes the argument-variables and their derivatives as abstract raising and lowering operators,

[TABLE]

because partial differentiation then imposes on these operators the conventional commutation algebra

[TABLE]

where ${\delta}_{pq}$ is the Kronecker $\delta$.

Generating functions, which are polynomials multiplying the number 1, are written as the action of the raising operators on a formal right-hand null state, while the projection operator that takes the trace of a generating function with an integral is written as a left-hand null state:

[TABLE]

where ${\delta}^{P}\!\left(z\right)$ is the Dirac $\delta$ in $P$ dimensions, and the inner product of the null states is normalized: $\left(0\mid 0\right)=1$.

The basis for generating functions is the set of number states which are elementary monomials. For any vector ${\rm n}$,

[TABLE]

number states are eigenstates of the set of number operators $a^{\dagger}_{p}a_{p}$:

[TABLE]

In particular for use with CRNs, we note the role of lowering operators in extracting the truncated factorials of number arguments: for any non-negative integer $k$,

[TABLE]

where $\left|{\rm n}-k_{p}\right)$ is the number state with $k$ subtracted from ${{\rm n}}_{p}$ and all ${{\rm n}}_{q}$ for $q\neq p$ unchanged.

With these steps the generating function becomes a vector in a linear space:

[TABLE]

All number states are normalized with respect to the Glauber inner product, defined by

[TABLE]

and the Glauber inner product with a generating function is simply the trace of the underlying probability density:

[TABLE]

The above conventions define the standard representation in which we will work with generating functions in the remainder of this article.

II.4 The stochastic process for a CRN

We now apply the above operator formalism to the particular forms of transfer matrices and generating functions produced by CRNs, and show how the reaction and complex representations decompose the resulting representations of the generators.

Corresponding to the column vector ${\Psi}_{Y}$ of activities from Eq. (4), we introduce a second column vector ${\psi}_{Y}\equiv\left[{\psi}_{Y}^{i}\right]$ which takes as its argument the Doi lowering operators $a_{p}$ that extract the truncated factorials of ${{\rm n}}_{p}$ according to Eq. (19). Making use of the notation (11) for these truncated factorials, we may write the coefficients in these two vectors as

[TABLE]

There is a corresponding row vector for the adjoints: ${\psi}_{Y}^{\dagger}\equiv{\left[{{\psi}^{\dagger}}_{Y}^{i}\right]}^{T}$ is a row vector of components defined on the complex indices $i$, where the notation means

[TABLE]

(As we did with the notation (11) for factorial moments, we take the definitions (23,24) for $\psi$ and ${\psi}^{\dagger}$ as general forms, in which other arguments besides $a$ and $a^{\dagger}$, such as the first-moment value $n$ can appear. This is convenient when expressing the approximation made in the mass-action rate law, and exhibiting its relation to the exact equation of motion as an operator expression.)

We now use the above notations to group the sample numbers and index shifts that describe the concurrent conversion of reactants into products in the elementary reaction events of a CRN. We begin with the transfer matrix, and then show the simplifications afforded by working with the Liouville operator.

For proportional sampling without replacement, the number dependence of the probability for the reaction $\left(i,j\right)$ from state ${\rm n}$ is simply given by the component ${\Psi}_{Y}^{i}\!\left({\rm n}\right)$ from Eq. (23). The way the master equation acts on indices is slightly more complicated: For the delivery of probability into state ${\rm n}$, the master equation must sample both $\rho$ and ${\Psi}_{Y}^{i}$ at a value shifted from ${\rm n}$ by the stoichiometric coefficients that are consumed at complex $i$ minus those that are produced at complex $j$. The simplest way to express such shifts is to let the vector of shift operators $e^{\partial/\partial{\rm n}}$ (with the exponential evaluated component-wise on ${\rm n}$) serve as an argument to ${\psi}_{Y}^{i}$, thus:

[TABLE]

with this convention for both the rate constant and the shift operator, the matrix ${\rm T}$ from Eq. (5) can be written

[TABLE]

Here the dot-product ($\cdot$) between ${\psi}_{Y}\!\left(e^{\partial/\partial n}\right)$ and ${\Psi}_{Y}\!\left({\rm n}\right)$ in the first line indicates that these two vectors are to be multiplied component-wise with respect to the complex index $i$, so that their product is then extracted by the indicator functions $w_{i}^{T}$ in $\mathbb{A}$ from Eq. (3). The cross-term in ${\psi}_{Y}^{j}$ and ${\psi}_{Y}^{i}{\Psi}_{Y}^{i}$, made explicit in the second line of Eq. (26), performs the required index shift on ${\rm n}$ in both $\rho$ and ${\Psi}_{Y}$ to account for the particles lost from the system’s state through complex $i$ and those gained by the system’s state through complex $j$. The other cross-term, with ${\psi}_{Y}^{i}\!\left(e^{-\partial/\partial n}\right)$, simply cancels the shift operators in ${\psi}_{Y}^{i}\!\left(e^{\partial/\partial n}\right)$ and represents the loss of probability from state ${\rm n}$ with rate ${\Psi}_{Y}^{i}\!\left({\rm n}\right)$.

Working with the Liouville operator from Eq. (7) is much more straightforward, because the lowering operators $a_{p}$ both extract sample numbers and shift indices according to Eq. (19), so these do not need to be separately tracked as they are in the transfer matrix. Using the exact definitions (23,24) with the arguments $a$ and $a^{\dagger}$ implicit, $\mathcal{L}$ takes the form

[TABLE]

Eq. (27) is one of the central equations of this paper, and underlies many of the simplifications we present here. In this expression, all the formal asymmetry of the standard first-moment rate equations from the CRN literature has disappeared, and particle consumption and creation are now treated symmetrically. This is the first of many simplifications gained by working with the Laplace transform and the Doi operator algebra.

III The dynamics of moment hierarchies

From the foregoing constructions we can directly compute the equations of motion for arbitrary moments of the density ${\rho}_{{\rm n}}$. These equations are finitely generated if we work in a basis of factorial moments, the demonstration of which is the main result of this section.

III.1 Dynamics of factorial moments for a single species

For a non-negative integer $k$, the operator that extracts the truncated factorial ${{\rm n}}_{p}^{\underline{k}}$ from number-states, which are the basis for a general state vector $\left|\phi\right)$, is $a_{p}^{k}$. Therefore the expectation of ${{\rm n}}_{p}^{\underline{k}}$ in the state $\left|\phi\right)$ is given by

[TABLE]

and from Eq. (7) and the form (27) for $\mathcal{L}$, its time-dependence is given by

[TABLE]

In order to obtain a recursion relation for the time dependence of $\left<{{\rm n}}_{p}^{\underline{k}}\right>$ in terms of the values of other factorial moments, we must commute the product of lowering operators $a_{p}^{k}$ through all powers of raising operators, which are gathered in the coefficients of the row-vector ${\psi}^{\dagger}_{Y}$. The result of the commutation is a finite series with descending powers of $a_{p}$ and $a_{p}^{\dagger}$. For positive integers $k$ and $y$, the evaluation of operator products of powers of raising and lowering operators is given by777The proof is by induction. If $k\geq y$, start with an elementary evaluation of $a_{p}^{k}{a^{\dagger}}_{p}$ and then induct on $y$. If $k\leq y$, start with an elementary evaluation of $a_{p}{a^{\dagger}}_{p}^{y}$ and induct on $k$.

[TABLE]

The first line emphasizes the symmetric roles of $k$ and $y$ in the combinatorial coefficient. In the second line we have used the definition (11) applied to $y^{\underline{j}}$ (rather than ${{\rm n}}^{\underline{j}}$) to simplify the index of summation in the case that $k>y$.

Now we may expand the evaluation appearing in Eq. (29), using the sum (32)

[TABLE]

as an operator identity acting on general states. Here we have used the property of the Doi operator algebra that commutation through the exponential $e^{\sum_{q}a_{q}}$ shifts all $a_{p}^{\dagger}\rightarrow\left(a_{p}^{\dagger}+1\right)$, after which all factors of $a_{p}^{\dagger}$ annihilate the right ground state $\left(0\right|$. Thus we have eliminated all factors of $a^{\dagger}$, along the way extracting from ${\psi}_{Y}^{\dagger}$ the factorial moments ${\left(y_{p}^{i}\right)}^{\underline{j}}$ of the stoichiometric coefficients.

The powers of the lowering operator $a_{p}^{k-j}$ in Eq. (38) are the same form as terms already in ${\psi}_{Y}$, so we can absorb them into ${\psi}_{Y}$ by shifting the stoichiometric coefficients in row $p$, which we denote as

[TABLE]

$Y+{\left(k-j\right)}_{p}$ is the matrix in which the ${}_{p}^{i}$ component is $y_{p}^{i}+k-j$, $\forall i$, and $y_{q}^{i}$ is unchanged for $q\neq p$.

From these evaluations we can re-express Eq. (29) as

[TABLE]

Here $Y_{p}^{\underline{j}}$ is a row vector in which the $i$th component is the factorial moment ${\left(y_{p}^{i}\right)}^{\underline{j}}$. Note that $Y_{p}^{\underline{0}}$ is the row vector of 1s, $Y_{p}^{\underline{1}}=Y_{p}$, the $p$th row of $Y$, etc. The first line of Eq. (48) contains the full sum over $j$ from Eq. (38), and the second line uses the fact that $1^{T}\mathbb{A}\equiv 0$ to eliminate the $j=0$ term. In the third line, we have expressed the shift of coefficients in the $p$th row of $Y$ again using an exponential shift operator denoted $e^{\partial/\partial Y_{p}}$, which acts on all components $y_{p}^{i}$. This will be a convenient notation for working with moments involving multiple species.

III.2 The generator of the stochastic process acting on the

moment hierarchy

For multiple species, we generalize the notation to an integer-valued vector of powers $k\equiv\left[k_{p}\right]$, and arrange summation indices similarly in vectors $j\equiv\left[j_{p}\right]$. Recalling the definition (12) of the factorial moment hierarchy $\Phi$, we write its time derivative as

[TABLE]

The notation ${\dot{\prod}}_{p}$ denotes a product over species $p$ within each index $i$ of the row vectors $Y_{p}^{\underline{j_{p}}}$. (Note that in the sums over $j_{p}$, we must now retain the $j_{p}=0$ entries, because even if one index $Y_{p}^{\underline{j_{p}}}={\left[1\right]}^{T}$, there may be others in the sum where $j_{p^{\prime}}\neq 0$, and the product $\left({\dot{\prod}}_{p}Y_{p}^{\underline{j_{p}}}\right)\mathbb{A}$ is only assured to vanish when all $j_{p}=0$.) Now the shift operator $e^{{\left(k-j\right)}_{p}\,\partial/\partial Y_{p}}$ in the first line of Eq. (57) offsets all coefficients $y_{p}^{i}$ in the row $p$ by ${\left(k-j\right)}_{p}$, so we write the matrix $Y$ in the second line with its rows shifted uniformly by the entries in the column vector $\left(k-j\right)$.

The third line of Eq. (57) gives the definition of the matrix $\Lambda$ introduced in Eq. (13). Each entry in the matrix ${\Psi}_{Y+\left(k-j\right)}\!\left({\rm n}\right)$ is itself a truncated factorial, so the expression is closed on $\Phi$. Although the orders $k_{p}$ may be arbitrarily large, the matrix $\Lambda$ has only finitely-many nonzero entries, limited by the largest values of $j_{p}$ for which the rows $Y_{p}^{\underline{j_{p}}}$ have non-vanishing entries (recalling the definition (11) for truncated factorials).

Eq. (57) is the main result with which we will work in this paper. It contains the rate equations for the species numbers $n$ studied by Feinberg, Horn and Jackson, and extends these to give a compact representation for the dynamics of all higher-order moments as well. In Sec. VI we illustrate graphic methods for representing $\Lambda$.

IV Network topology, deficiency, and the classification and role

of different network flows

The results up to this point are true for a general CRN. Thus they say nothing directly about the topological properties that may afford simplifications such as uniqueness and positivity of steady states or factorability of distributions. In this section we shift to a consideration of topology and its implications, including the concept of deficiency and the Feinberg Feinberg:def_01:87 and Anderson-Craciun-Kurtz (ACK) theorems Anderson:product_dist:10

IV.1 Deficiency, and a basis to decompose network flows

The connectivity of the adjacency matrix $\mathbb{A}$, together with the stoichiometric matrix $Y$, determines the linear subspace of ${\rm n}$-values that can be accessed through any flow on the network, called the stoichiometric subspace. The number and character of steady states depends on whether the ${\rm n}$-dependence of the mass-action rate equations within this subspace admits a Lyapunov function Horn:mass_action:72 ; Baez:QTSM:17 . That, in turn, depends on whether all flows that transport net matter into or out of any complex must also transport some net matter into or out of some chemical species, thus increasing its chemical potential in a direction that opposes the flow. If so, then all the flows are mean-regressing and strictly positive steady states are unique. If not, there are net fluxes at the complexes that do not lead to net fluxes of species, and for these there is no force leading to mean-regression. In the latter case, multiple steady states or non-strictly positive steady states888These would be boundary solutions where some concentrations equal zero. cannot be ruled out Craciun:multistability:05 ; Craciun:bistability:06 .

The character of the dynamical steady states therefore depends on the relative dimension of $\ker\mathbb{A}$ (the flows that absorb or emit no material at the complexes) and $\ker Y\mathbb{A}$ (the flows that absorb or emit no material at the chemical species), within the stoichiometric subspace. A sketch of the demonstration that this is a topological characteristic, following Feinberg Feinberg:notes:79 ; Feinberg:def_01:87 , Horn and Jackson Horn:mass_action:72 , follows:

For any CRN, the stoichiometric subspace corresponds to

[TABLE]

The dimension of $S$, and of the subspace of flows through complexes that do not produce motions within $S$, are denoted

[TABLE]

$\delta$ is the deficiency of the CRN.

The following relations hold (as identities) among dimensions in $Y\mathbb{A}$ and $\mathbb{A}$:

[TABLE]

If $C$ is the total number of complexes, then it follows that

[TABLE]

The expression for $\dim\left(\ker\left(\mathbb{A}\right)\right)$ is simple, and follows from the fact that $\mathbb{A}$ functions as an ordinary graph Laplacian on each connected component of the complexes. The argument involves the following observations:

Weak reversibility and linkage classes: Connected components in the simple graph that includes only complexes and reactions are termed linkage classes in the CRN literature. Weak reversibility is the condition that any node in a linkage class can be reached from any other by some sequence of reactions. The subset of complexes in a linkage class, which can be reached starting from any complex, and which subsequently are never exited, is called a strong terminal linkage class. Weak reversibility of the whole CRN ensures that each linkage class is a strong terminal linkage class, and (as in Gunawardena Gunawardena:CRN_for_bio:03 ) we will limit to this case for simplicity.999The more general case differs only by superficial book-keeping to exclude complexes that are exited permanently.

The counting rule for deficiency: Weakly reversible processes are ergodic on each linkage class, so by the Perron-Frobenius theorem or an equivalent argument Baez:QTSM:17 , there is one basis vector for $\ker\left(\mathbb{A}\right)$ for each linkage class. Let $l$ denote the number of linkage classes. Then $l=\dim\left(\ker\left(\mathbb{A}\right)\right)$, and the Feinberg counting result that

[TABLE]

follows from Eq. (61).

Complex-balanced steady states: Flows in $\ker\left(A\right)$ are termed complex-balanced, because they require no net transport of flux to or from any complex to the species that make it up. All steady states must (tautologically) be in $\ker\left(Y\mathbb{A}\right)$. If $\delta=0$ the two spaces have the same dimension and thus are the same. In other words, the steady-state condition that there be no sources or sinks at species nodes entails the condition that there be no sources or sinks at complexes. In this case a convexity argument Feinberg:notes:79 ; Baez:QTSM:17 implies that there is a unique steady state in the positive orthant for any value of the rate constants.

Nonzero deficiency: If $\delta>0$ there are species-balancing flows that are not complex-balancing, and $\dim\left(\ker\left(\mathbb{A}\right)\right)$ is larger than the number of constraints from $\partial n/\partial\tau=0$ (which is only $s$). Steady states then generally exist out of the subspace of $\ker\left(\mathbb{A}\right)$.

IV.1.1 Using mean-regression as a basis to decompose flows

beyond the ($\delta=0$)-condition

We noted above that a basis for $\ker\left(A\right)$ has one vector on each linkage class. At the level of individual linkage classes, these are profiles in $\Psi$ proportional to the maximal eigenvectors of the Perron-Frobenius theorem for simple diffusion under $\mathbb{A}$. To characterize the remainder of the space of activities on the complexes, we require a basis that decomposes the pre-image of the stoichiometric subspace from the remaining flows.

Let ${\left\{e_{\alpha}\right\}}_{\alpha=1}^{s}$ be a basis for ${\ker\left(Y\mathbb{A}\right)}^{\perp}\subseteq{\mathbb{R}}^{C}$. Changes in $\Psi$ along these directions lead to changes in flows $\partial n/\partial\tau$ within $S$, which we term $s$-flows.

Let ${\left\{{\tilde{e}}_{\beta}\right\}}_{\beta=1}^{\delta}$ be a basis for $\ker\left(Y\mathbb{A}\right)/\ker\left(\mathbb{A}\right)$. Changes in $\Psi$ along these directions lead to changes by species-balanced but not complex-balanced flows, which do not alter $\partial n/\partial\tau$ and which we term $\delta$-flows.

If follows that jointly $\left\{{\left\{e_{\alpha}\right\}}_{\alpha=1}^{s},{\left\{{\tilde{e}}_{\beta}\right\}}_{\beta=1}^{\delta}\right\}$ form a basis for ${\ker\left(\mathbb{A}\right)}^{\perp}\subseteq{\mathbb{R}}^{C}$.

This basis leads to a decomposition of $\mathbb{A}$ and therefore of $\mathcal{L}$, different from either the reaction or the complex representations in Eq. (3), which we might term the “stoichiometric” representation:

[TABLE]

The $s$-flows and $\delta$-flows distinguished in the stoichiometric representation turn out to make dimensionally different contributions to the moment equations (57). The original Feinberg result, expressed in this language as saying that the $s$-flows completely determine the first-moment equations of motion – and that these are the only such flows for CRNs of deficiency zero – extends to a claim that the $s$-flows dominate the behavior of moments at all orders lower than the mean particle numbers, for general CRNs where they are not the only flows. We explore the consequences of these scaling dimensions in Sec. V.

IV.2 The Feinberg deficiency-0 theorem and the

Anderson-Craciun-Kurtz theorem

Feinberg’s Deficiency-0 theorem was originally framed Feinberg:notes:79 as a result in the convex analysis of the mass-action rate law, which is generally assumed to be a mean-field approximation. The Anderson-Craciun-Kurtz theorem Anderson:product_dist:10 uses the Feinberg existence proof for solutions to the mass-action rate law, but in proving that the underlying distributions have product-Poisson form,101010under the sampling model assumed throughout this paper it actually strengthens Feinberg’s original theorem: because factorial moments of Poisson distributions have exactly the relation to the first moment assumed in mean-field formulations, the Feinberg solution to the mass-action rate equations need not be framed as a result in (unregulated) mean-field approximation, but can rather be seen as an exact result.

The Doi operator algebra, together with the stoichiometric decomposition in Eq. (63), provides an elegant way to both prove the ACK theorem and see that it makes the Feinberg theorem exact when $\delta=0$, and also see why and how these results no longer hold when $\delta>0$.

We begin with the set of cases $\sum_{p}k_{p}=1$ of Eq. (48), in which the sum on $j$ contains the single term $j=k$; these are the equations of motion for the set of first moments. For species $p$:111111Note how, from the symmetric form of $\mathcal{L}$ in Eq. (27), the asymmetry of the usual rate equations has resulted from the projection onto the dynamics of a particular moment.

[TABLE]

So far we work at general deficiency, but because $Y\mathbb{A}{\tilde{e}}_{\beta}\equiv 0,\forall\beta$ by construction, only the basis elements corresponding to $s$-flows are non-zero.

Eq. (64) is exact for this stochastic process, and it is almost the same as the standard expression for the mass-action rate law, except that it involves an expectation of the observables ${\Psi}_{Y}\!\left({\rm n}\right)$, which may include higher-order correlations in ${\rm n}$. Arbitrarily ignoring these correlations, and replacing $\left<{\Psi}_{Y}\!\left({\rm n}\right)\right>$ with ${\psi}_{Y}\!\left(\left<{\rm n}\right>\right)$, defines the mean-field approximation.

The one case where the mean-field form is exact is when the state $\left|\phi\right)$ equals some coherent state $\left|\xi\right)$. The coherent states are the generating functions of Poisson distributions, constructed in the Doi algebra as

[TABLE]

$\xi\equiv\left[{\xi}_{p}\right]$ is a vector of the mean particle numbers ${\xi}_{p}=n_{p}$ for the $P$ chemical species.

Coherent states are eigenstates of the Doi lowering operator, and thus

[TABLE]

giving the mean-field form as an exact result:

[TABLE]

Using the basis $\left\{e_{\alpha}\right\}$ to handle the counting of dimensions in the stoichiometric subspace (not all linear combinations of $n_{p}$ are necessarily dynamic in a particular CRN), the condition that $\partial n/\partial\tau=0$ in Eq. (64) becomes

[TABLE]

These are the set of equations proved by Feinberg to have a unique, strictly-positive solution when $\delta=0$, for all non-degenerate values of the rate constants. Thus, if the distribution $\left|\phi\right)$ is a coherent state, the Feinberg result is an exact solution.

The coherent state identified by Eq. (68) is the ACK solution if it is a solution at all. However, whereas any mean-field equation can be solved for some coherent-state parameter, the corresponding state is only a solution to the whole moment hierarchy if it also leads to stasis of all higher moments. To see why this is assured at $\delta=0$ and not otherwise, we insert the stoichiometric decomposition (63) into the general moment hierarchy (57) to obtain:

[TABLE]

On a coherent state, the vector of factorial moments $\left<{\Psi}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>$ differs from the value $\left<{\Psi}_{Y}\!\left({\rm n}\right)\right>$ for the first-moment condition only by an overall factor $\prod_{p}{\xi}_{p}^{{\left(k-j\right)}_{p}}$. If $\delta=0$, only $s$-flows are present in the decomposition, and the set of projections $\left\{e_{\alpha}^{T}\left<{\Psi}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>\right\}=\prod_{p}{\xi}_{p}^{{\left(k-j\right)}_{p}}\times\left\{e_{\alpha}^{T}\left<{\Psi}_{Y}\!\left({\rm n}\right)\right>\right\}$ vanishes exactly when Eq. (68) holds, for all values of $k$ and $j$. Although the truncated factorials of stoichiometric coefficients $\left[\dot{\prod_{p}}Y_{p}^{\underline{j_{p}}}\,\right]$ generally differ from $Y$, this changes only the weight of the inner products with the $e_{\alpha}$ and thus the strength with which each $s$-flow contributes to the rate equation away from the steady state. Hence the coherent-state solution identified by the first-moment condition (68) is a steady-state solution for the whole moment hierarchy, proving the ACK theorem.

If $\delta>0$, then $\delta$-flows also exist in the sum (73). However, whereas $Y\mathbb{A}{\tilde{e}}_{\beta}\equiv 0,\forall\beta$, $\left[{\dot{\prod}}_{p}Y_{p}^{\underline{j_{p}}}\,\right]$ will not generally project out $\mathbb{A}{\tilde{e}}_{\beta}$ in the higher-order terms $\sum_{p}j_{p}>1$, except possibly in special cases of fine-tuning of the rate parameters. Therefore, in general, $k$-dependent linear combinations of ${\tilde{e}}_{\beta}^{T}\left<{\Psi}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>$ and $e_{\alpha}^{T}\left<{\Psi}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>$ will be required to vanish at a steady state, obviating any simple Poisson solution.

V Scaling regimes, matched asymptotic expansions, and the

controlling role of mean-regressing flows

The series expansion in factorials of the stoichiometry from Eq. (57), together with the interpretation of each term as a projection operator for some product-Poisson distribution that we used to prove the ACK theorem in Sec. IV.2, provides a way to associate different regions in the lattice of factorial moments with control by different subsets of the network flows. Both a set of practical solution methods for steady states, based on matched asymptotic expansions, and the concept of approximating a complex distribution locally by a product-Poisson distribution, follow when we recognize that different terms in the series expansion are associated with different scaling behaviors because they capture distinct combinations of the rate constants from the network.

The next three sub-sections cover the following topics in order:

Descaling: We first introduce the descaling of the moment hierarchy with coherent-state parameters. Since, for Poisson distributions these are the only scale parameters, descaling turns the moment hierarchy for a coherent state into a vector of 1s.121212This is true as long as the descaling is done with the coherent state’s own $\xi$ values. More generally the moment hierarchy becomes a geometric progression.

Matched asymptotic expansions: Generalizing the Feinberg steady-state condition to all orders, as we did above to prove the ACK theorem, allows us to use the condition $\Lambda\Phi=0$ as a recursion relation on $k$ to solve the moment hierarchy, much like the solution of any Laplace equation, with $\Lambda$ serving as the Laplacian on the lattice of moments. When we do this with the descaled moment hierarchy, it becomes easy to show that recursion upward in any component of $k$ produces a convergent power-series expansion on the range $k/\xi\ll 1$ (where $\xi$ stands for whichever particle number corresponds to the component of $k$ being incremented), whereas recursion downward in $k$ is convergent for $k/\xi\gg 1$. This suggests a general method of solution for moment hierarchies using matched asymptotic expansions, where the matching conditions are imposed in the region $k\sim\xi$.

$1/n$-expansion about Poisson backgrounds: The same small-parameter recursion that controls the asymptotic expansion for $k/\xi\ll 1$ also shows a sense in which the mean-regressing flows (the $s$-flows) define a leading-Poisson approximation to low-order moments for a general CRN. This is true even when the projection onto the basis $\left\{e_{\alpha}\right\}$ does not define a zero-deficiency sub-network of the original CRN. We can construct a linear combination of the coherent-state solutions to the first-moment steady-state conditions for which the remainder term that must be added to obtain an exact solution makes a contribution that is $\mathcal{O}\!\left(k/\xi\right)$ smaller than the contribution of the Poisson backgrounds for the low-order moments (those with $k/\xi\ll 1$).

V.1 Descaling the moment equations with coherent-state

parameters

Just as $Y_{p}$ projects out the ACK product-Poisson for $\delta=0$ networks, each of the projectors $\left[{\dot{\prod}}_{p}Y_{p}^{\underline{j_{p}}}\,\right]$ in the moment-recursion equation (57) projects out some product-Poisson distribution if values for the corresponding coherent-state parameters can be found. (If they are not unique, it can project out more than one such solution.) We may choose to reference exact solutions for ${\Phi}_{k}$ to locally-chosen Poisson distributions in different regions of $k$ corresponding to different terms $j\equiv\left[j_{p}\right]$, and let the recursion equations solve for the (smaller) deviations from these reference-Poissons.

For any vector $\xi\equiv\left[{\xi}_{p}\right]$ of mean values, we may descale the activities ${\Psi}_{Y}^{i}$ from Eq. (23) as

[TABLE]

If we normalize the ${\Psi}_{Y}$ vectors in this way, a corresponding counter-normalization of the adjacency matrix can be defined as

[TABLE]

In cases where we wish to use the stoichiometric representation (63), a similar descaling of the corresponding projection vectors is

[TABLE]

for all $\alpha$ and $\beta$. All these are chosen so that

[TABLE]

Applying the scale transformations (74,75) to the equations of motion (57) for the moment hierarchy $\Phi$ gives the equation for a descaled hierarchy $\hat{\Phi}$, in which the leading geometric dependence on $\xi$ has been factored out:

[TABLE]

Following Feinberg but extending his consideration to all moments, we try to construct steady states, for which Eq. (83) can be used as a recursion relation among the moments of ${\hat{\Phi}}_{k}$.

V.2 Matched asymptotic expansions for the steady-state

condition

The terms that govern the behavior of recursions in the components of $k$, if Eq. (83) is used to (exactly or approximately) solve for ${\hat{\Phi}}_{k}$, are combinations of the form

[TABLE]

(with $k$, $j$, $Y$, and $\xi$ carrying indices for each $p$, which we suppress to reduce clutter). The ratios $Y^{\underline{j}}/j!$ are fixed parameters of $\hat{\Lambda}$ and in any case only finite in number. The ratios that govern scaling behavior across the moment hierarchy are the terms $k!/\left[\left(k-j\right)!\,{\xi}^{j}\right]$. Because the stoichiometric coefficients and therefore the limits in the sums over $j$ are finite, it is possible to consider a range of typical particle number $n\sim\xi\gg\max\left(j\right)$, in which $k!/\left[\left(k-j\right)!\,{\xi}^{j}\right]\sim{\left(k/\xi\right)}^{j}$ in the ranges that govern the transition between scaling regions. These scale factors govern the stability of asymptotic expansions as follows:

For any fixed value of $k$, increasing $j$ in the sum (83) lowers the order of all moments in $\left<{\hat{\Psi}}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>$, at the same time multiplying the corresponding term by a coefficient $\sim{\left(k/\xi\right)}^{j}$. Let $\underline{j}$ be the smallest value at which ${\dot{\prod}}_{p}\left(Y_{p}^{\underline{j_{p}}}/{\xi}_{p}^{j_{p}}\right)\hat{\mathbb{A}}$ does not vanish. (In general, this occurs when $\sum_{q}j_{q}=1$, so exactly one of the terms $Y_{p}^{\underline{j_{p}}}=Y_{p}$ and $Y_{q}^{\underline{j_{q}}}=1$ for all other $q\neq p$.)

To extend the recursion upward by one order, we must increment $k$ while holding $j$ fixed at $\underline{j}$. The new moments appearing at order $k$ are referred to those at the immediately preceding order in the recursion by higher-order terms $j>\underline{j}$ in the sum at the current $k$. The relative magnitude of the preceding terms to the new terms scales as $\sim{\left(k/\xi\right)}^{j-\underline{j}}$. For $k\ll\xi$, successively higher-order terms are expressed as sums of lower-order terms with positive powers of $\left(k/\xi\right)$, consistent with both a non-zero radius of convergence, and with damping-out of uncertainties in the initial conditions of the recursion. (The latter property is important for an asymptotic expansion to provide a robust solution algorithm.)

For $k\gg\xi$ the opposite is true: the lower-order terms must be solved as functions of the higher-order terms, which are multiplied by positive powers of $\left(\xi/k\right)$. Thus in this range the downward recursion is consistent with a non-zero radius of convergence, and damps out uncertainties in the starting conditions assumed at large $k$.

This argument is the basis for a solution in terms of matched asymptotic expansions, where stable recursions are carried out starting respectively from $k=0$ (up-going) and from asymptotically large $k$ (down-going), and matching conditions are imposed in the overlap region $k\sim\xi$, which are marginally stable for both series. An interesting feature of this solution is that, for CRNs where the mean-field approximation predicts multiple steady states, there may still be unique large-$k$ asymptotic behaviors required to ensure boundedness of moments at all orders. In such cases, it is the downward recursion from large $k$ that anchors the solution to the moment hierarchy. This is a counter-intuitive result given the conventional mean-field approach to moment closure, which attempts to anchor all higher-order moments in solutions to the first moments, but as a consequence cannot obtain the ergodic sum over multiple steady states which is a property of the exact all-orders solution.

We do not offer a formal proof that these asymptotic expansions can be consistently performed for all CRNs and all dimensionalities of the moment hierarchy, which is an exercise beyond the scope of the current paper. However, in one dimension, the recursion is elementary to define, and for higher-dimensional systems we offer examples of decompositions of the solution for which numerical simulation suggests that a similar expansion can be used.

V.3 Leading Poisson approximations to nonzero-deficiency CRNs

The above analysis of the scaling of terms in an asymptotic expansion for solutions to $\hat{\Lambda}\hat{\Phi}=0$ has an immediate corollary: in general steady-state solutions, a basis of product-Poisson distributions associated with $s$-flows dominates the low-order moments. The kernel of the argument is that, although the projection operators associated with both $s$- and $\delta$-flows share in the same scaling, at the lowest order where the upward-going recursion begins, the $\delta$-flow contributions are projected out – recall that their absence from the first-moment conditions is their defining feature – therefore only $s$-flow contributions serve as seeds for the polynomial expansion in $\left(k/\xi\right)$. We now demonstrate that relation:

The steady-state condition for Eq. (83), with the stoichiometric decomposition inserted from Eq. (63), becomes

[TABLE]

The lowest-$k$ conditions that the moment hierarchy must satisfy are the first-moment conditions, which are the set of terms $\sum_{q}k_{q}=1$ and $j=k$, for which exactly one $Y_{p}^{\underline{j_{p}}}=Y_{p}$ and $Y_{q}^{\underline{j_{q}}}=1$ for all other $q\neq p$, as noted in the previous section.

Write the state vector $\left|\phi\right)$ for a general solution to Eq. (89) as a sum

[TABLE]

in which $\left\{{\xi}^{\left(\gamma\right)}\right\}$ is the set of all mean-field solutions to the first-moment steady state conditions, and $c_{\gamma}$ are coefficients to be determined. By construction $e_{\alpha}^{T}\sum_{\gamma}c_{\gamma}{\psi}_{Y}\!\left({\xi}^{\left(\gamma\right)}\right)=0$, $\forall\alpha$.131313Note that $\left|{\phi}^{\prime}\right)$, as the generating function for a difference of distributions, will not generally be derived from any distribution with all positive values. Refer to the corresponding expectations as

[TABLE]

and likewise for higher-order moments and descaled moments $\hat{\Psi}$.

Under the decomposition (90), the steady-state condition (89) becomes

[TABLE]

If there is a unique steady-state solution, and if the descaled moment hierarchy $\hat{\Phi}$ is descaled with this $\xi$, then by construction it will be the case that ${\left<{\hat{\Psi}}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>}^{\left(\gamma\right)}=c_{\gamma}\left[1\right]$ (the vector of all 1s), $\forall k,j$ on the right-hand side of Eq. (101), which contains a single term in the sum on $\gamma$.141414Note that, even in the case of a unique steady state, we cannot presume that $c_{\gamma}=1$ unless $\left<{\Psi}_{Y}\right>$ includes a term proportional to ${\Phi}_{0}\equiv\left<1\right>$, because the first-moment condition does not otherwise fix the normalization of the geometric sequence within the total distribution. In the more general case, we can choose the descaling parameters so that ${\left<{\hat{\Psi}}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>}^{\left(\gamma\right)}=c_{\gamma}\left[1\right]$ for a particular $\gamma$ of our choice. (Generally this means descaling whichever term makes the largest contribution on the right-hand side, to remove $k$- and $j$- dependence in that term.) The argument that the term $\left|{\phi}^{\prime}\right)$ is sub-leading is then made in two steps:

For some set of coefficients $\left\{c_{\gamma}\right\}$ we can ensure that $e_{\alpha}^{T}{\left<{\Psi}_{Y}\!\left({\rm n}\right)\right>}^{\prime}=0$. This is because the set of steady-state solutions for ${\xi}^{\left(\gamma\right)}$ form a basis for the set of all solutions to the first-moment steady-state conditions. In general (even if there is only one solution for $\xi$), the required coefficients $\left\{c_{\gamma}\right\}$ may need to be determined by matching conditions to a large-$k$ asymptotic expansion.151515The mono-stable solution in Sec. VI.2 illustrates the need for a non-trivial normalization in the case of a unique solution, and the bistable solution in Sec. VI.3 illustrates the case of solution for a linear combination of Poisson backgrounds. (Note that the values of ${\tilde{e}}_{\beta}^{T}{\left<{\Psi}_{Y}\!\left({\rm n}\right)\right>}^{\prime}$ are unconstrained at order $\sum_{p}k_{p}=1$, and must be determined as part of the recursion on $k$.) 2. 2.

For $\sum_{p}k_{p}>1$, the leading-order dependence on ${\left<{\hat{\Psi}}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>}^{\left(\gamma\right)}$ on the right-hand side of Eq. (101) comes when $\sum_{p}j_{p}=2$, and by Eq.(84) this term is $\mathcal{O}\!\left[{\left(k/\xi\right)}^{2}{\left<{\hat{\Psi}}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>}^{\left(\gamma\right)}\right]$. The leading term on the left-hand side arises where $\sum_{p}j_{p}=1$, involves only the $s$-flows, and is $\mathcal{O}\!\left[\left(k/\xi\right){\left<{\Psi}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>}^{\prime}\right]$. Hence we conclude that, at the level of counting naive scaling dimensions, there is a perturbative expansion in small $\left(k/\xi\right)$ about a sum of coherent states, in which ${\left<{\Psi}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>}^{\prime}\sim\mathcal{O}\!\left[\left(k/\xi\right){\left<{\hat{\Psi}}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>}^{\left(\gamma\right)}\right]$. The reason it is meaningful to make such a scaling comparison, when the $j$ values used to estimate the powers of scale factors are different on the left-hand and right-hand sides of Eq. (101), is that we have been free to choose the de-scaling parameter for $\hat{\Phi}$ so that for whichever $\gamma$ gives the largest contribution, ${\left<{\hat{\Psi}}_{Y+\left(k-j\right)}\!\left({\rm n}\right)\right>}^{\left(\gamma\right)}=c_{\gamma}\left[1\right]$, removing $k$- and $j$-dependence from the source on the right-hand side of the equation.

This completes the argument.

VI Worked examples

We now demonstrate the above results for moment hierarchies in a cascade of examples. Each successive example increases the generality of the problem and introduces a new feature of general solutions of CRN’s. Explicit constructions involving transfer matrices or Liouville operators are given in the main text where they first occur, and the corresponding forms that differ only by elaboration for later examples are removed to App. A.

VI.1 A CRN with 1-species, 2 states and no conserved quantities

This model is a minimal non-trivial form for a CRN, showing how uniqueness of positive steady states follows from deficiency-0, and exhibiting the proof of the ACK theorem in terms of coherent-state projection operators from Sec. IV.2. The CRN is given by the graph introduced in Fig. 1 with the associated reaction scheme (1). Its mean-field rate equation is

[TABLE]

The master equation, illustrating the decomposition (26) for the transfer matrix, is

[TABLE]

For the transfer matrix (103), the Liouville operator is

[TABLE]

The first line is a direct translation of the reaction-representation from Eq. (3) for the conversion of particles in each unidirectional reaction. The second line extracts the overall factor of the projection operator $\left(\alpha-\beta a\right)$ that vanishes on a coherent state with parameter $\xi=\alpha/\beta$ in all moment equations, which is the proof of the ACK theorem given in Sec. IV.2.

VI.2 A CRN with 1-species, 3 states and

no conserved quantities

In this section we introduce a one-parameter family of models which have the same rate equation over the entire family. Fig. 2 depicts a limiting member which lacks weak reversibility. Although the limiting case is formally outside the scope of the assumptions in the rest of the paper, it is useful to highlight the role of $\delta$-flows in driving solutions away from the Poisson form associated with the ACK theorem. Weak reversibility may be established without changing the rate equation by adding two reactions to the graph of Fig. 2, to obtain the family (over the rate parameter $\epsilon$) of graphs shown in Fig. 3, which we analyze below.

For this model, we demonstrate the stoichiometric decomposition of the Liouville operator from Eq. (63), and descaling of the moment recursion equation. For the parameters we will use in simulations, a single, downward-going asymptotic expansion (obtained in Krishnamurthy:CRN_moments:17 ) is sufficient to solve the entire moment hierarchy to arbitrary precision, starting from an analytically derived large-$k$ limiting form. However, we will also demonstrate the upward/downward matched asymptotic expansion to illustrate the stability properties of the recursion in small-$k$ and large-$k$ ranges. Exact solutions to moment hierarchies of this kind are only possible for birth-death vanKampen:Stoch_Proc:07 type CRNS with one species, of which this class of models is an example, or for $\delta=0$ CRN’s, whereas the asymptotic expansions have a much wider applicability, as we demonstrate in later sections.

The rate equation for the CRN’s in both Fig.  2 and Fig. 3 is

[TABLE]

which differs from Eq. (102) only in changing the activities that govern particle creation and destruction.

The simplest CRN graph with rate equation (105) is shown in Fig. 2, with associated reaction scheme

[TABLE]

The model introduces competing auto-catalysis at two orders: particle creation occurs in proportion to the density of existing particles, while particle destruction occurs in proportion to the cube of the density. Because the CRN in Fig. 2 has only three complexes, however, the rate equation has only two non-trivial roots, and therefore cannot support multiple positive steady states. (It does, however, have the marginally stable steady state $n\equiv 0$.)

The CRN of Fig. 3 with reaction scheme

[TABLE]

These may be checked from Eq. (62) to have deficiency $\delta=1$.

The master equation is provided in Eq. (239), and the associated Liouville operator is

[TABLE]

The limit $\epsilon\rightarrow 0$ is degenerate with the graph of Fig. 2. The important feature of this Liouville operator is that the two projection terms, $\left(\alpha-\epsilon a\right)$ and $\left(\epsilon-\beta a\right)$, are now multiplied by distinct non-trivial operators (respectively $1$ and $\left(a^{\dagger}\!a-1\right)$, and cannot both be made to vanish independently at a single Poisson solution at general values of $\epsilon$. This is the way in which deficiency-1 (attained by adding a complex within a linkage class) moves the CRN outside the scope of the ACK theorem.

VI.2.1 The stoichiometric decomposition

We can use the property suggested by the concept of deficiency – the categorization of flows as mean-regressing versus non-mean-regressing – to further clarify how the non-independence of the projection terms in the Liouville operator (108) results in deviations from Poisson steady-state form.

The three decompositions of the Liouville operator (108), from Eq. (3) and Eq. (63), are given by

[TABLE]

The top line is the reaction representation (since the reactions are bi-directional, we have combined both departure terms in the row vectors). The middle line is the complex representation. In this representation it is clear why, if $\epsilon\rightarrow 0$, the reaction cannot be complex-balanced: all terms in a given row vector have the same sign, so any positive density produces non-zero flows at some complexes. The bottom line is the stoichiometric representation. The first dyadic corresponds to the $s$-flow $\mathbb{A}ee^{T}$, and the second dyadic corresponds to the $\delta$-flow $\mathbb{A}\tilde{e}{\tilde{e}}^{T}$.

Note that the diagonalization of the $s$-flow couples activity in the complex $a$ to changes of probability across complexes ${a^{\dagger}}^{2}$ and ${a^{\dagger}}^{3}$, and vice-versa with activity at $a^{3}$ and change of probability across complexes $a^{\dagger}$ and ${a^{\dagger}}^{2}$. Thus, despite the similarity in form to the deficiency-0 projector in the Liouville operator from Eq. (104), the $s$-flow projection in Eq. (LABEL:eq:Lv_123_R) cannot be written as a stand-alone Liouville operator from a deficiency-0 sub-network of the current network.

To illustrate the way in which different combinations of $s$- and $\delta$-flows control the scaling of ${\Phi}_{k}$ in different regions, we note the forms of projection operators at different orders in the sum (57):

[TABLE]

The lowest-order term $Y^{\underline{1}}\mathbb{A}$ projects out solutions ${\Phi}_{k+2}/{\Phi}_{k}=\alpha/\beta$, while the highest-order term $Y^{\underline{3}}\mathbb{A}$ projects out the solution ${\Phi}_{k+2}/{\Phi}_{k+1}=\epsilon/\beta$. These turn out to be the two limiting moment ratios, respectively, in the limits $k=1$ (the moment recursion formula has no term at $k=0$) and $k\rightarrow\infty$, as we now demonstrate.

VI.2.2 Scaling behavior of the rate equation used as a

recursion relation

The Poisson background in the expansion (90), projected out by $Y^{\underline{1}}\mathbb{A}$ and ensuring vanishing of the $s$-flow contribution to the moment dynamics at each order $k$, is given by {\psi}_{Y}={\left[\begin{array}[]{ccc}\xi&{\xi}^{2}&{\xi}^{3}\end{array}\right]}^{T}, where the mean number $\xi$ satisfies

[TABLE]

This CRN has a unique steady state, so the terms appearing in Eq. (101) are:

[TABLE]

As the numerical evaluations below will show, this is a model in which, despite uniqueness of the Poisson solution matching the first and third moments, the overall normalization of the moment hierarchy is not anchored in the lowest term ${\Phi}_{0}$, so a relative normalization $c_{0}$ for the ${\psi}_{Y}\!\left(\xi\right)$ contribution is undetermined. In this way, the role of the Poisson background in an exact solution of the moment hierarchy is different from a mean-field approximation. MFT would require ${\left<{\rm n}\right>}^{2}=\alpha/\beta$ in place of Eq. (192), which requires ${\left<{\rm n}\right>}^{2}={\left(c_{0}\xi\right)}^{2}$. The freedom for $\xi$ to differ from $\left<{\rm n}\right>$ by the normalization $c_{0}$ is necessary, because even for this simple network, the MFT prediction for the mean is not valid.

The exact recursion relation for the deviations from Poisson moments is then the finite sum

[TABLE]

For parameters that produce suitably small mean particle number, the recursion relation implied by Eq. (199) may be solved for all $k$ to any desired precision, from the upper asymptotic behavior alone, as shown in Krishnamurthy:CRN_moments:17 . We are interested here, however, in understanding the small-$k$ and large-$k$ behavior of the moments. To this effect, we can see directly both aspects of the scaling presented in Sec. V. As an asymptotic expansion, the recursion relation specifies the higher-order difference $\left({\hat{\Phi}}_{k+2}-{\hat{\Phi}}_{k}\right)$ as a power series in $\left(k/\xi\right)$ with coefficients from values and differences of $\Phi$ at lower $k$ indices.

The seed for the expansion at orders ${\hat{\Phi}}_{4}^{\prime}$ and higher from the Poisson background is the $k$-independent value $c_{0}\left(\epsilon/\sqrt{\alpha\beta}-1\right)$, multiplied by a polynomial of $\mathcal{O}\!\left(k/\xi\right)$. Note, however, that the terms ${\hat{\Phi}}_{2}^{\prime}$ and $c_{0}$, which are undetermined by the first-moment condition $e^{T}\left<{\Psi}_{Y}\!\left({\rm n}\right)\right>=0$, enter the recursion relation according to the scaling of the overall asymptotic expansion, and are permitted to be $\mathcal{O}\!\left(1\right)$ relative to $\left(\epsilon/\sqrt{\alpha\beta}-1\right)$.

Finally, we observe that for a fine-tuned value of the rate parameters $\epsilon=\sqrt{\alpha\beta}$, the correction term can be made to vanish, and the ACK-like solution projected to zero by the $s$-flow term $\mathbb{A}ee^{T}$ in Eq. (LABEL:eq:Lv_123_R) becomes a steady-state solution, even though for this CRN $\delta=1$.

Next we illustrate how the asymptotic expansion with a requirement of boundedness at large $k$ anchors the moment hierarchy at all orders. We seed the downward-going asymptotic expansion with the leading non-constant approximation to the recursion relation around the limiting ratio projected out by $Y^{\underline{3}}\mathbb{A}$ from Eq. (191), which has the form

[TABLE]

The corresponding leading-order approximation for large-$k$ moments may be written about any reference value $k_{0}$ as

[TABLE]

where $\mathcal{N}$ is an overall normalization to be determined.

Fig. 4 shows a comparison of the numerical regression from Eq. (199), both for an upward-going recursion from $k=1$, and for a downward-going recursion with the large-$k$ asymptotic seed (201), to estimates of the first 10 moments from a Gillespie simulation of the underlying process. The parameters used in the demonstration are: $\alpha=100$, $\beta=10$; $\epsilon=70$. So the relevant parameters are $\xi=\sqrt{\alpha/\beta}=\sqrt{10}$; $\epsilon/\beta=7$; $\epsilon/\sqrt{\alpha\beta}=\epsilon/\beta\xi\approx 2.2136$. We have normalized the constant $c_{0}$, which is unspecified by the recursion relation to the simulated mean $\left<{\rm n}\right>$, and find close agreement to all other the moment ratios.

VI.3 A CRN with 1-species, 4 states, 2 linkage classes and no

conserved quantities

We can preserve the number of $s$-flows and $\delta$-flows from the previous model, but introduce the possibility for multi-stability, by increasing both the number of complexes and the number of linkage classes by one. The resulting graph for a minimal model with this elaboration is shown in Fig. 5, and its reaction scheme is given by

[TABLE]

The corresponding rate equation is

[TABLE]

We have labeled the reaction rate constants in this model to reflect a set of cases that are often important in biochemistry and industrial synthesis: a reaction that is (directly or indirectly) self-catalyzed by feedback through the synthetic network is the main channel for production and decay of the product (rate constants $k_{1}$ and ${\bar{k}}_{1}$), compared to an uncatalyzed pathway that has nonzero but small rate ($\epsilon$), while a significant rate ($k_{2}$) remains for spontaneous decay of the product.

Quadratic-order autocatalysis in this CRN comes from the same pair of reactions as it does in Fig. 3. The addition of a fourth state to the complex graph creates a cubic first-moment rate equation and thus the possibility for multiple steady states.

For the graph of Fig. 5, the Liouville operator is

[TABLE]

The stoichiometric decomposition of this operator is similar to that from the previous model, and is given in Eq. (LABEL:eq:Lv_one_species_bistable).

The factorials $Y^{\underline{j}}$ and projection operators $Y^{\underline{j}}\mathbb{A}$ appearing in Eq. (57) are given by

[TABLE]

The lowest-order (in $k/\xi$) projector in Eq. (213) is $Y^{\underline{1}}\mathbb{A}$, which is the projection operator corresponding to the $s$-flow in Eq. (LABEL:eq:Lv_one_species_bistable) in the appendix. If required to vanish on a coherent state, it gives

[TABLE]

For appropriate parameter choices, this may have either a unique stable solution, or it may have three solutions, two stable and one between them that is unstable. Unlike the model of Fig. 3, both ${\Phi}_{0}$ and ${\Phi}_{2}$ have non-zero coefficients in $\left<{\Psi}_{Y}\!\left({\rm n}\right)\right>$, so the normalization of the mean is fixed relative to ${\Phi}_{0}\equiv 1$. Solutions of the form (90) must satisfy $\sum_{\gamma}c_{\gamma}=1$, and therefore we may set ${\left<{\Psi}_{Y}\!\left({\rm n}\right)\right>}^{\prime}\equiv 0$.

In this model (contrasted with the result in Eq. (191)), both of the projection operators $Y^{\underline{2}}\mathbb{A}$ and $Y^{\underline{3}}\mathbb{A}$ cancel the same ratio ${\Phi}_{k+2}/{\Phi}_{k+1}=k_{1}/{\bar{k}}_{1}\equiv K_{1}$, so there are only two scaling behaviors expressed in the model, respectively at $k\rightarrow 0$ and $k\rightarrow\infty$.

This CRN, also being a birth-death type process, can be solved exactly Anderson:Lyapunov:15 for the steady state. As for the 3-state model, this moment hierarchy may also be solved by recursion from an upper asymptotic limit that is derivable analytically (though again the numerical calculation is stable only for sufficiently small mean particle numbers.) If we define $\left<{{\rm n}}^{\underline{k}}\right>/\left<{{\rm n}}^{\underline{k-1}}\right>\equiv R_{k}$, then the set of $R_{k}$ must obey the recursion

[TABLE]

When we solve the recursion (215) numerically, directly in terms of moments ${\Phi}_{k}$, we begin with a more refined large-$k$ approximate form than the first-order approximation used as a seed in Eq. (200). The second-order leading non-constant approximation, corresponding to the form (201) given for the previous model, is

[TABLE]

where $\mathcal{N}$ is an arbitrary normalization to be fixed by ${\Phi}_{0}=1$, and

[TABLE]

The large-$k$ asymptotic behavior of $R_{k}$ in Eq. (215) can likewise be solved for an expansion in $1/k$ about the leading fixed point, in the same manner as Eq. (200) for the previous model. In this case, the leading departure is $\mathcal{O}\!\left(1/k^{2}\right)$ rather than $\mathcal{O}\!\left(1/k\right)$ as in the 3-complex model. In general, it can be shown that that the first departure from whatever fixed point is dictated by the leading large-$k$ projection operator in Eq. (57) is determined by the highest power of $k$ appearing in the expansion in the sum. This order corresponds to the largest stoichiometric coefficients for that component of $k$ appearing in the CRN.

VI.3.1 Bistability in MFT and handling mixtures of Poisson

basis elements around $k=0$

The handling of multi-stability in CRNs with $\delta>0$ introduces several new interesting properties, both within the moment recursion relations and in their relation to mean-field theory. First, MFT will generally predict multi-stability for Poisson solutions, whether or not the mean particle number is large enough that trajectories in the stochastic process actually generate a multi-modal density of particle numbers. The meaning of MFT solutions in relation to the analytic structure of representations of the generating function is an interesting topic from which we briefly draw results below, but mostly refer to other developments Cardy:Instantons:78 ; Smith:Signaling:11 ; Smith:LDP_SEA:11 ; Smith:evo_games:15 (see also Coleman:AoS:85 , Ch. 7). Second and more important, the moment relations are exact, and we therefore expect them to possess unique solutions corresponding to the ergodic distribution, even when mean particle number is large enough that the MFT representation of multistability corresponds to a true incipient161616We say “incipient” because in conventional usage, breaking of ergodicity is an asymptotic property, here in the scaling variable $\left<{\rm n}\right>$ as $\left<{\rm n}\right>\rightarrow\infty$. Formally, ergodicity breaking can still be considered well-defined even in finite-size systems, to the extent that it is associated with stationary paths in a semiclassical approximation that are essential singularities with respect to the asymptotic expansion in fluctuations Cardy:Instantons:78 . breaking of ergodicity. The expansion in high-order moments becomes an important if cryptic representation of the trajectories responsible for first-passages between domains.

We illustrate some of these properties for the case of bistability with a numerical example at parameters: $\epsilon=36$, $k_{2}=49$; $k_{1}=14$; ${\bar{k}}_{1}=1$. The two stable solutions to Eq. (214) are ${\xi}^{\left(1\right)}=1$ and ${\xi}^{\left(3\right)}=9$, and an unstable solution exists at ${\xi}^{\left(2\right)}=4$.

Fig. 6 shows the recursive solution for $\left<{\hat{{\rm n}}}^{\underline{k}}\right>$, descaled with ${\xi}^{\left(3\right)}=9$, starting from the large-$k$ asymptotic form (216). The solution exactly matches the expansion (90), with the coefficients given by

[TABLE]

Mean-field theory suggests no natural interpretation of the mixture (218) with a negative coefficient on an unstable solution. What would normally be done instead in MFT is to express the mean $\left<{\rm n}\right>\approx 1.56$ directly as a mixture of the two MFT-stable values ${\xi}^{\left(1\right)}=1$ and ${\xi}^{\left(3\right)}=9$ with a mixing coefficient

[TABLE]

We will use the phenomenological description (219) to understand qualitatively how MFT and stationary-point expansions relate to the exact solution of the moment hierarchy.

VI.3.2 Interpretation with a Kramers approximation for

first-passage times

The interpretation of the ergodic solution in terms of a sum over naïve mean-field backgrounds can be compared to a stationary-point expansion using the method of instantons, which is developed in Coleman:AoS:85 ; Smith:Signaling:11 ; Smith:evo_games:15 . Stationary-point locations and probabilities are governed by the minima of a non-equilibrium effective potential, which we have computed for this particular network in Smith:geo13:16 (Ch.7), and which takes the form171717In Ref. Smith:geo13:16 , we use the notation $\Phi\!\left(\bar{n}\right)$ for the effective potential, which we change here to $\Xi\!\left(\bar{n}\right)$ to avoid a collision with the notation for the moment hierarchy.

[TABLE]

The extrema of the effective potential are exactly the values of the Poisson parameters ${\xi}^{\left(\gamma\right)}$. (Here we arbitrarily set the zero of the effective potential to $n=4$, the saddle point.) A plot of the effective potential versus $n$ is shown in Fig. 7.

The probability to occupy either minimum may be approximated by the Kramers formula Sjogren:Kramers:16 derived from $\Xi$,

[TABLE]

which follows from a semiclassical approximation to the escape rates by the non-trivial stationary trajectories known as instantons. Fig. 7 shows that for these parameters, the minima of $\Xi$ are $\approx-0.25$ and $\approx-0.19$, respectively at $\bar{n}=1$ and $\bar{n}=9$, so the Kramers approximation is not expected to be quantitatively accurate. The corresponding second derivatives ${\Xi}^{\prime\prime}$ take values $\approx 0.07$ and $\approx 0.02$.

Fig. 8 shows a timeseries for the particle number ${\rm n}$ in a Gillespie simulation, with dashed lines indicating the minima of the effective potential (220). Some features of the stationary-point approximation are reflected: a majority of the timeseries remains near ${\rm n}={\xi}^{\left(1\right)}=1$ with strong mean regression, while excursions with modest persistence and wider fluctuations occur out to ${\rm n}\sim{\xi}^{\left(3\right)}=9$. However, the excursions do not have the character of fully metastable equilibria. The log of the empirical stationary distribution from the simulations is monotone decreasing, with a visible shoulder (the signature that excursions are persistent) with a mode around ${\rm n}\approx 15$.

Eq. (221) gives, for the occupation probabilities of the two states, approximate values

[TABLE]

in which $p_{\bar{9}}$ corresponds roughly to the empirical mixing coefficient $c_{\rm eff}$ in Eq. (219). The Kramers formula captures the larger weight on ${\left<{\Psi}_{Y}\!\left({\rm n}\right)\right>}^{\left(1\right)}$, but over-estimates the admixture of ${\left<{\Psi}_{Y}\!\left({\rm n}\right)\right>}^{\left(3\right)}$ by about a factor of three.

This CRN was also studied by Anderson et al. Anderson:Lyapunov:15 . Using the fact that the steady state probability $\rho_{\rm ss}$ is known, they showed that the non-equilibrium potential, defined here as $-\log(\rho_{\rm ss})$, converges to the Lyapunov function for the corresponding deterministic dynamics, in an appropriate scaling limit. The resulting Eq. (38) in Anderson:Lyapunov:15 is our Eq. (220) obtained as a large-deviation function.

VI.4 Two species, cross-catalysis, and loss of factorability

The final model we will develop shows the loss, for $\delta>0$, of the factorability which characterizes the steady states of deficiency-0 CRNs under the ACK theorem. We retain the properties already developed, of deviation from Poisson statistics, and the capacity for multistability, by simply changing the autocatalytic feedback in the model of Fig. 5 to a cross-catalytic feedback between two symmetric chemical species.

The resulting multistable network, for two species ${\rm A}$ and ${\rm B}$, is shown in Fig. 9, and its reaction scheme is given by

[TABLE]

The (now vector-valued) rate equation takes the form

[TABLE]

The cross-catalytic CRN from Fig. 9 has Liouville operator

[TABLE]

We introduce pairs of raising and lowering operators $\left(a^{\dagger},a\right)$, $\left(b^{\dagger},b\right)$, for the species ${\rm A}$ and ${\rm B}$ respectively. The convention we adopt for ordering the components of the (somewhat complicated) vectors $\psi\!\left(a,b\right)$ and ${\psi}^{\dagger}\!\left(a^{\dagger},b^{\dagger}\right)$ is given in Eq. (279), and the corresponding stoichiometric decomposition of the Liouville operator is given in Eq. (348).

VI.4.1 Mean-field solutions and scaling regimes

The truncated factorials $Y_{a}^{\underline{j_{a}}}$ and $Y_{b}^{\underline{j_{b}}}$ that govern the scaling regimes in the steady-state moment hierarchy are provided in Eq. (363). The lowest-order terms in the descaled form (83) of the moment equation are

[TABLE]

The two vanishing conditions for $Y_{a}\hat{\mathbb{A}}$ and $Y_{b}\hat{\mathbb{A}}$ (coming from the projectors for the two $s$-flows of this network) are

[TABLE]

These are solved at ${\xi}_{a}={\xi}_{b}=\xi$, with $\xi$ again satisfying the mean-field equation (214) for the single-species model of Sec. VI.3.

In the other asymptote, the highest-order terms in Eq. (83) that are not identically zero in the component-wise product of $Y_{a}^{\underline{j_{a}}}\cdot Y_{b}^{\underline{j_{b}}}$ are

[TABLE]

The vanishing conditions for these two projectors are

[TABLE]

They are again solved at ${\xi}_{a}={\xi}_{b}=\xi$ but now with $\xi=k_{1}/{\bar{k}}_{1}\equiv K_{1}$, reproducing the large-$k$ asymptotic condition from the single-species model of Sec. VI.3.

This two-species case may again be solved for mixed moments in the neighborhood of the diagonal $k_{a}=k_{b}$, writing coupled recursion relations for the ratios of the factorial moments, as shown in Krishnamurthy:CRN_moments:17 . In the following section, we illustrate an alternate solution method using the asymptotic expansions that we have developed in the earlier sections.

VI.4.2 Polynomial expansion of a solution for the moment

equation in a neighborhood of the diagonal $k_{a}=k_{b}$

We now illustrate how the representation $\Lambda$ of the generator for the stochastic process, acting similarly to a Laplacian on the two-dimensional lattice of moments ${\Phi}_{\left(k_{a},k_{b}\right)}$, can be approximately solved in a neighborhood of the diagonal $k_{a}=k_{b}$. The method of solution is to use the symmetry of the recursion equations under $k_{a}\leftrightarrow k_{b}$ to expand solutions in even powers of $\left(k_{a}-k_{b}\right)$, with coefficient functions of $\left(k_{a}+k_{b}\right)$ solved by asymptotic expansion in a manner similar to that used in the 1-species models of Sec. VI.2 and Sec. VI.3.

We do not have a proof that the radius of convergence of these solutions covers the entire lattice of $k$ values, but comparisons to Gillespie simulation show good agreement in neighborhoods of the diagonal, suggesting that the asymptotic boundary conditions we use are consistent with those of full solutions. Existence of approximate solutions with this form shows a strong breaking of factorability from the product-Poisson form that is a characteristic of the ACK solution for deficiency-zero networks.

In the following solutions, two combinations of the rate constants that will appear repeatedly are given a short-hand:181818The second term, $\eta$, has already appeared in Eq. (217).

[TABLE]

As in the asymptotic solutions for the 1-species models, we begin by recognizing that the large-$k$ asymptotic form is dominated by the scaling of the projection operator that is non-zero for the largest value of $j$ in the sum (83). This is the projector given in Eq. (228). Therefore we define the descaled moment operator by ${\Phi}_{\left(k_{a},k_{b}\right)}\equiv K_{1}^{k_{a}+k_{b}}{\hat{\Phi}}_{\left(k_{a},k_{b}\right)}$, and look for solutions in the form

[TABLE]

where ${\varphi}_{\left(k_{a},k_{b}\right)}\rightarrow 0$ at large $k_{a}$ or $k_{b}$.191919We justify this assumed scaling by reference to the large-$k$ limit (216) from the similar 1-species model, because the orders of catalysis are more similar to that case than to the model of Sec. VI.2 leading to the soft (logarithmic) divergence of Eq. (201).

We introduce diagonal and transverse variables, written as functions of the vector argument $k$:

[TABLE]

In matrix multiplications below, we will often use $\kappa$ and $q^{2}$ as function names, with the argument $k$ which is the index of summation suppressed as in usual matrix notation. We look for solutions to $\varphi$ in the form of power series,

[TABLE]

Each term is to be chosen so that ${\varphi}^{\left(\alpha\right)}_{\kappa}\rightarrow 0$ as $\kappa\rightarrow\infty$. The functions ${\varphi}^{\left(\alpha\right)}_{\kappa}$ obey recursion relations similar to those for an infinite sequence of 1-dimensional moment hierarchies labeled by $\alpha$, except that the vectors in the sequence are coupled across values of $\alpha$. Within the solution for each ${\varphi}^{\left(\alpha\right)}_{\kappa}$, we may treat $\kappa$ itself as the discrete index of the recursion. Here as in the 1-species models, the descaled recursion relation suggests leading-order asymptotics for ${\varphi}^{\left(\alpha\right)}_{\kappa}$ in powers of $1/\kappa$, which may be used to seed numerical solutions.

The conversion from the original lattice ${\hat{\Phi}}_{k}$ to the sequence of vectors ${\varphi}^{\left(\alpha\right)}_{\kappa}$ leads to the following approximation procedure to solve for steady states: The steady state condition from Eq. (83) is

[TABLE]

where $\underline{1}\equiv\left[1\right]{\left[1\right]}^{T}$ is the dyadic matrix of all 1s. We introduce a zeroth order source term $s^{\left(0\right)}_{\kappa}$ defined by

[TABLE]

so the steady-state condition is equivalent to the series solution of an inhomogeneous Laplacian equation

[TABLE]

in which $\hat{\Lambda}$ serves as Laplacian and $s^{\left(0\right)}_{k}$ is the source for the inhomogeneous solution.

The form of $\hat{\Lambda}$ can be described graphically in terms of difference operators acting across adjacent positions on the lattice of $k$ values, as shown in Fig. 10. $\hat{\Lambda}$ acts non-trivially on $q^{2\alpha}$ as well as on ${\varphi}^{\left(\alpha\right)}_{\kappa}$, so Eq. (236) induces connections across orders in $\alpha$, and we relegate the details of a solution by successive approximations to App. B.

VI.4.3 Properties of steady states in the 2-species model

The major features of the steady-state solution in this model, which we have verified against Gillespie simulations, are the following:

Order of terms versus $\alpha$: The naïve scaling dimensions implied for ${\varphi}^{\left(\alpha\right)}_{\kappa}$ by Eq. (83) suggest that these functions should decay at large $\kappa$ with increasing powers of $1/\kappa$. Numerically, this appears to be borne out, with indeed the entire series ${\varphi}^{\left(\alpha\right)}_{\kappa}$ decreasing in magnitude with increased $\alpha$. In addition to the on-diagonal terms (where $q^{2}\equiv 0$), which are defined entirely in terms of ${\varphi}^{\left(0\right)}_{\kappa}$, terms adjacent to the diagonal, which should be dominated by ${\varphi}^{\left(\alpha\right)}_{\kappa}$ at low orders in $\alpha$, are well approximated by the solution ${\varphi}^{\left(0\right)}_{\kappa}$ across the whole range of $k_{a}=k_{b}$.

Scaling of finite-order approximations along rays of $\left|q\right|/\kappa$: The measure of error – non-zero values of $\partial{\Phi}_{k}/\partial\tau$ – appears roughly constant along rays of fixed $\left|k_{a}-k_{b}\right|/\left(k_{a}-k_{b}\right)\equiv\left|q\right|/\kappa$ at finite orders of approximation in ${\varphi}^{\left(\alpha\right)}_{\kappa}$. Stabilizing the asymptotic expansion independently at each order of ${\varphi}^{\left(\alpha\right)}_{\kappa}$ becomes increasingly difficult as $\alpha$ increases, due to cross-level feedback and the successive-approximation algorithm we use for solution. Thus we obtain an approximate solution only through order ${\varphi}^{\left(5\right)}_{\kappa}$.

Comparison of the 2-species cross-catalytic and 1-species autocatalytic models: The ratios of adjacent moments in the value $k_{a}+k_{b}$ – which require at minimum comparing on-diagonal and first-off-diagonal moments – are shown in Fig. 11 and compared to the corresponding sequence of ratios derived from the moment solutions in Fig. 6. We find that the mean value $\left<{{\rm n}}_{a}\right>\equiv\left<{{\rm n}}_{b}\right>$ in the 2-species model is very close to the mean value $\left<{\rm n}\right>$ from the 1-species model, as suggested by the equivalence of their mean-field forms, even though both models differ significantly from the MFT-approximation, which is the solution $\xi$ to Eq. (214). Moreover, the second moments $\left<{{\rm n}}_{a}{{\rm n}}_{b}\right>$ remain close to the 1-species expectation $\left<{\rm n}\left({\rm n}-1\right)\right>$, and again different in both cases from the MFT prediction. At higher $k$, a different behavior is seen: the transition to scaling dominated by the term (228) is governed in the 2-species model by $k_{a}$ and $k_{b}$ comparable to 1-species $k$, and not by the sum $k_{a}+k_{b}$. This is expected by comparing the forms of the two Liouville operators (204) and (225).

VI.4.4 Breaking of the factorability of the ACK theorem

through cross-catalysis

The solution scheme defined in Sec. VI.4.2 and worked out in App. B.1 suggests that the moment hierarchy near the diagonal is well approximated (at least at large $k$) by a function of $k_{a}+k_{b}$, which is a strong deviation from the factorability that would be produced by the ACK theorem for a deficiency-0 network. To study the failure of factorability more directly than through the numerical approximation scheme of App. B.1, we may alternatively approximate the large-$k$ behavior of the moment hierarchy by a sum of products of powers of $1/k_{a}$ and $1/k_{b}$, solved fully for $\hat{\Lambda}\hat{\Phi}=0$ order-by-order in $1/k$. The two expansions do not have the same asymptotics along the boundaries $k_{a}=0$ and $k_{b}=0$, but they can be made to satisfy the same boundedness criteria at large $k$ in a neighborhood of $k_{a}=k_{b}$.

An expansion of the solution to $\hat{\Lambda}\hat{\Phi}=0$ to second order in $1/k$ is given by

[TABLE]

which we have checked numerically cancels the error term $\partial\hat{\Phi}/\partial\tau$ to the correct order of $1/k$.202020 Eq. (237) may be compared to Eq. (216) for the one-species model, and also to the leading-order scaling estimate for the term ${\varphi}^{\left(0\right)}_{\kappa}$ from Eq. (389). Along the diagonal $k_{a}=k_{b}$, Eq. (237) becomes

${\Phi}_{\left(k_{a},k_{b}\right)}\approx\mathcal{N}K_{1}^{k_{a}+k_{b}}\left[\frac{4\eta\left(k_{a}+k_{b}\right)}{{\left(k_{a}+k_{b}+K_{1}-\eta\right)}^{2}}.\right]$

The leading behavior differs from Eq. (389) by a factor $\left(4/3\right)$ multiplying $\left(K_{1}-\eta\right)$, which is consistent with the fact that the scaling solutions in App. B.1 only propagate the effects of $\hat{\Lambda}$ upward in a hierarchy of powers; feedbacks down the hierarchy are absorbed in higher-order correction terms that have the same large-$k$ order as corrections in the multiplier of $\left(K_{1}-\eta\right)$. The first departure from factorability occurs in the second-order term with numerator $-\eta K_{1}$, showing where the $\delta$-flow contributions create correlated fluctuations that would be ruled out in a ($\delta=0$)-network.

VI.4.5 Generalizing to a larger number of species

Our constructions apply to CRNs with arbitrary numbers of species, but the foregoing models show how the character of solution methods changes with increasing numbers. Because the generator (57) has finitely-many terms for any finite CRN, for one-species problems, the number of undetermined boundary data that must be sampled to search for stable asymptotic expansions is always finite. Moreover, the large-$k$ limit may be extended to improve the precision of approximations from coarser seed functions. In two or more dimensions, unknown boundary data can exist along all surfaces of co-dimension 1 or more, in which one or more $k_{p}=0$. If an asymptotic bounding surface is moved outward for the nonzero values of $k_{p}$, new unknown values are added to the set that must be sampled along small-$k_{p}$ boundaries. Therefore increasingly much of the information in a solution must come from boundary-condition data, compared to the constraint in the scalar condition $\Lambda\Phi=0$.

Other problems of convex analysis, analogous to the Feinberg deficiency-zero argument, are also left as questions for future work. Is there a systematic way to represent the number of distinct scaling regions controlled by terms in the sum (57) over $j$? Do large-$k$ asymptotic conditions on the vector $\xi$ of coherent-state parameters always possess unique solutions? When are they underdetermined, and in these cases do the solutions from $s$-flow conditions extend outward indefinitely?

Despite leaving several detailed questions to be addressed, we emphasize that the finite rank of the operator $\Lambda$ reduces the solving of a moment hierarchy to all orders, to a problem of equal complexity to solving a Laplacian diffusion equation, generally with beyond-nearest-neighbor couplings. This is a simpler and less-costly problem than direct simulation, especially for high-order moments in systems with large particle numbers.

VII Conclusions

The new observations and results of this work may be grouped under the following four main topics:

A shift in emphasis from topology to dynamics: In this article we have bypassed the use of deficiency as a topological index to categorize networks, and focused instead on the contrast between mean-regressing and non-mean-regressing flows on the complex graph, which is the dynamical property that causes deficiency to be important. The dynamical distinction, which we express in the stoichiometric representation (63) of the stochastic-process generator, continues to be definable in terms of the images and kernels of $\mathbb{A}$ and $Y\mathbb{A}$, even in ($\delta>0$)-networks where it cannot be associated with a deficiency-0 sub-network. As shown in Sec. V.3, the unique role of the mean-regressing flows as the determinants of the first-moment conditions, which persists to all orders in deficiency-0 networks, persists in a more limited form as the leading term in a $1/n$-expansion for the determination of moments in more general networks.

The Doi operator algebra, Laplace transforms, and an expanded role for Poisson basis distributions: The Doi operator algebra used to express the generators of the stochastic process is the tool that allows us to associate the $s$-flows and $\delta$-flows in the stoichiometric representation with corresponding product-Poisson distributions, the origin and meaning of which are the same as those of the unique steady-state distributions in the Anderson-Craciun-Kurtz theorem. In this way not only the linear algebra of first moments in the stoichiometric subspace, but the Poisson family of distributions as basis functions, extends directly from deficiency-0 to deficiency-nonzero cases. This simplification and clarification results from working with the Laplace transform and the Liouville operator: the elementary projection operators in terms of which $\mathcal{L}$ naturally decomposes, which annihilate particular Poisson distributions, describe collective motions that recursively relate all orders in the moment hierarchy.

The manifestation of fundamental symmetries in the generator of the stochastic process: The Doi operator algebra also exposes symmetries of the generators (27) for stochastic CRNs that are obscured in the combination of index shifts and number-dependence of rates in the master equation (26), and masked entirely in the asymmetric form of the Mean-Field mass-action rate equations (64,67). We use the simplifications this formalism affords to both derive the moment hierarchy as well as develop approximations to solve them.

Locality of scaling regions and connection to the Poisson approximations: The finiteness of the generator acting on the moment hierarchy is the feature that allows scaling regimes to be defined locally in different asymptotic ranges of moments, and that makes Poisson-form basis distributions good approximations to the moment recursion relations in such regions. Despite the fundamental underlying complexity of CRNs – remember that the search for non-locally defined topological properties such as shortest paths or feedback cycles can be NP-hard Andersen:NP_autocat:12 – the convergence of both high- and low-$k$ asymptotic expansions toward the matching region buffers the strength with which different regions are coupled, and allows the convergence toward the Poisson basis elements to be locally governed. n this way, Poisson distributions which are among the lowest-information distributions, serve as a basis for the evaluation of the moments in systems with the potential for very high information capacity. As a by-product, we also obtain a systematic approach to moment closure, which does not have any of the problems of ad-hocness and unphysical results Schnoerr:Survey:15 prevalent in existing schemes.

Acknowledgments

DES thanks the Physics Department at Stockholm University for support during visits in 2014 and 2016 when the bulk of this work was carried out. We also acknowledge Dan Rockmore, Scott Pauls, and Greg Leibon for hospitality and conversations on related topics in 2010, Artur Wachtel in 2015 and Nathaniel Virgo of ELSI in 2016.

Appendix A Supporting algebra for CRN examples

This appendix provides explicit forms for transfer matrices, Liouville operators, and truncated factorials of the stoichiometric matrix, for the CRN models in the main text.

A.1 Forms of the transfer matrices appearing in master

equations

Transfer matrices are given below for the indicated CRN models:

For the graph of Fig. 2 the evolution equation for the density that defines the transfer matrix has the form

[TABLE]

This equation illustrates in the simplest form how Poisson steady states come to be ruled out when the ACK theorem no longer applies. The shift operator acts by single units ${\rm n}\rightarrow{\rm n}\pm 1$, but the numerical factors in the rates for particle creation and annihilation differ by second-order terms in ${\rm n}$, which cannot be absorbed in any Poisson distribution.

For the weakly reversible graph of Fig. 3, the master equation adds a term affecting fluctuations though it preserves the first-moment rate equation:

[TABLE]

The master equation for the CRN of Fig. 5 is

[TABLE]

For the CRN of Fig. 9, ${\rm n}$ becomes a two-component index to $\rho$, and the master equation becomes

[TABLE]

A.2 Stoichiometric representation forms for Liouville

operators of the models

The stoichiometric decompositions of Liouville operators not explicitly written in the main text are provided here:

A.2.1 The 1-species, 4-complex model of

Sec. VI.3

The diagonalized form of the Liouville operator (204), corresponding to master equation (240) is

[TABLE]

in which the top line is the $s$-flow and the bottom line is the $\delta$-flow.

A.2.2 The 2-species, cross-catalytic model of

Sec. VI.4

For the 2-species model of Fig. 9, the vector of creation and annihilation operators with respect to which we will write $\mathcal{L}$ in matrix form is:

[TABLE]

Then $\mathcal{L}$ from Eq. (225), corresponding to the master equation (241), in matrix form and also diagonalized, becomes

[TABLE]

The top line is the $s$-flow projected out by $\left(a^{\dagger}a+b^{\dagger}b\right)$; the second line is the $s$-flow projected out by $\left(a^{\dagger}a-b^{\dagger}b\right)$. The next two lines are the $\delta$-flows with the same symmetry or antisymmetry.

The stoichiometric vectors appearing in Eq. (83), and their products with $\mathbb{A}$ are

[TABLE]

Appendix B Polynomial expansion for solutions to the two-species CRN

of Sec. VI.4

This appendix describes the combination of recursive solution, and successive approximation, used to solve the hierarchical expansion (233) for the lattice of moments in the two-species model of Sec. VI.4.2.

Here, to simplify notation and improve readability, we will regard matrices such as $\hat{\Lambda}$ as operators that shift the indices $\left(k_{a},k_{b}\right)$ in terms ${\varphi}_{k}^{\left(\alpha\right)}q^{2\alpha}$ by means of discrete index-shift operators $e^{\partial/\partial k_{a}}$, $e^{\partial/\partial k_{b}}$, as we did for transfer matrices in Eq. (26) et seq. Because of the exchange symmetry in the dynamical equations under $k_{a}\leftrightarrow k_{b}$, $\hat{\Lambda}$ acts on ${\varphi}_{k}^{\left(\alpha\right)}$ through a shift of the $\kappa$ value (by integers), and on $q^{2\alpha}$ through shifts in $\alpha$. This allows us to treat $\kappa$ as an index shifted by integers, analogous to $k$ in single-species models. Where we suppress subscript $\kappa$ indices, the whole vector is intended.

With these notational conventions, the action of the generator in Eq. (236) can be broken down into three terms:

[TABLE]

Here ${\hat{\Lambda}}^{0}$ is a diagonal operator (in $\alpha$) acting only on ${\varphi}^{\left(\alpha\right)}$, which takes the form (refer to Fig. 10):

[TABLE]

In the same way as the action of $\hat{\Lambda}$ on the constant background $\underline{1}$ produced the zero-th order source $-\eta s^{\left(0\right)}$ in Eq. (236), the action on each order ${\varphi}^{\left(\alpha\right)}$ generates a source term $s^{\left(\alpha+1\right)}$ in Eq. (364) that propagates cross-terms one order upward in $\alpha$, defined by

[TABLE]

In addition to an upward-propagating “source” term, Eq. (364) contains “feedback” terms, which propagate cross-terms downward from order $q^{2\alpha}$ to all lesser orders $q^{2\beta}$ with $\beta<\alpha$, given by

[TABLE]

In Eq. (387), referring to the graphical form of Fig. 10, all terms from $j_{a}+j_{b}=1$ are grouped in the first line, all terms from $j_{a}+j_{b}=2$ are grouped in the second line, and all terms from $j_{a}+j_{b}=3$ are grouped in the third line.

B.1 Solution by upward propagation in powers of $q^{2}$ and

perturbative correction in $1/\kappa$

The decomposition (364) of the steady-state condition can now be solved by alternating steps of exact cancellation of source terms at ascending orders in $\alpha$, and successive approximation to cancel feedback terms which takes the form of a perturbative expansion in $1/\kappa$.

Upward-propagation consists of solving a series of Laplacian equations in the operator ${\hat{\Lambda}}^{0}$ in terms of the sources (366), as

[TABLE]

for all $\alpha\geq 0$.

Examination of the scaling terms in $\kappa$ in Eq. (LABEL:eq:L0_indexed_form) suggests that a bounded large-$\kappa$ asymptotic approximation for each order is

[TABLE]

Solution of Eq. (388) at each order $\alpha$ by one-dimensional matched asymptotic expansion proceeds as for the 1-species models. At order $\alpha=0$, the solution can be stably extended down to $\kappa=0$, but for $\alpha>1$ the asymptotic approximation (389) is not sufficient to produce convergence below $\kappa\sim 20$, and a non-trivial matched expansion is needed to produce valid higher-order corrections in $q^{2}$ to the low-order moments.

One round of forward propagation (through all orders $\alpha$) will produce a solution that does not satisfy $\hat{\Lambda}\sum_{\alpha=0}^{\infty}{\varphi}^{\left(\alpha\right)}q^{2\alpha}-s^{\left(0\right)}=0$, but rather

[TABLE]

The full solution can be approached perturbatively by using the feedback terms on the right-hand side of Eq. (390) as sources for an iterative correction to the original $\varphi$.

A successive-approximation approach to full solutions

${\varphi}^{\left(\alpha\right)}$

A scaling analysis following from Eq. (389) suggests that the correction terms needed to cancel the residuals in Eq. (390) are suppressed by powers of $1/\kappa$, and thus that a method of successive approximations should converge.

No feedback term ${\sigma}_{0}^{\left(\alpha\right)}$ to $q^{0}$ order from $\alpha=0$ exists, because $\beta\leq\alpha-1$ in the sum (364). For all $\alpha\geq 1$, ${\sigma}_{0}^{\left(\alpha\right)}$ is given by

[TABLE]

Since the lowest-order contribution is from ${\varphi}^{\left(1\right)}_{\kappa}\sim 4/{\kappa}^{3}$, it follows that ${\sigma}_{0}^{\left(\alpha\right)}$ is no larger for any $\alpha$ than the leading term

[TABLE]

one order lower (in either $K_{1}$ or $\kappa$, according to the range of $\kappa$) than the scaling of $s^{\left(0\right)}$ in Eq. (235). The same argument extends to higher $\beta$ in Eq. (390); the first term that contributes at each order scales with two additional powers of $1/\kappa$, and so is smaller than the corresponding $s^{\left(\alpha\right)}$ term in the first iteration of Eq. (366).

We therefore introduce a second-order correction term ${\varphi}^{\left(\alpha\right)\prime}$, satisfying

[TABLE]

where ${\varphi}^{\left(\alpha\right)\prime}$ is solved by upward propagation in $q^{2\alpha}$, as in Eq. (388), but now with an entire tower of sources $-\sum_{\alpha=\beta+1}^{\infty}{\sigma}^{\left(\alpha\right)}_{\beta}$ at each order $q^{2\beta}$, rather than just the zeroth-order source $s^{\left(0\right)}$ that was used for ${\varphi}^{\left(\alpha\right)}$, and leaving its own residues ${\sigma}^{\left(\alpha\right)\prime}$. After infinitely many iterations, the residue terms go to zero as a sequence in powers of $1/\kappa$, and if the sequence converges, the sum over corrections will be a solution to the original steady-state condition (236).

An equivalent expression for the closed solution (summing over all orders of perturbative correction), expressed in terms of the homogeneous operator ${\hat{\Lambda}}^{0}$, is

[TABLE]

Eq. (394) is similar in form to a Schwinger-Dyson equation for Green’s function solution, in which $\left({\hat{\Lambda}}^{0}{\varphi}^{\left(\alpha\right)}-s^{\left(\alpha\right)}\right)$ serves as a “bare” Green’s function, which defines a basis for perturbative incorporation of an “interaction” term $\sum_{\alpha=\beta+1}^{\infty}{\sigma}_{\beta}^{\left(\alpha\right)}$.

B.2 Numerical evaluations

We have implemented the above solution method for a series (233) truncated at five successive orders of approximation ${\varphi}^{\left({\alpha}_{\rm max}\right)}$ for ${\alpha}_{\rm max}=0,\ldots,4$. Under perfect convergence of the perturbative recurrence (393), the error measure $\partial{\hat{\Phi}}_{\left(k_{a},k_{b}\right)}/\partial\tau$ would cancel to order $q^{2\left({\alpha}_{\rm max}+1\right)}$ around the diagonal (order $q^{10}$ for the highest-order approximation we compute). We obtain cancellation along the diagonal to machine precision for the lowest-order correction ${\varphi}^{\left(0\right)}$, but instability of the downward-going asymptotic expansion at higher orders degrades both the accuracy of ${\hat{\Phi}}_{\left(k_{a},k_{b}\right)}$ on the diagonal for $k_{a}+k_{b}\lesssim 20$, and convergence toward the $q^{2\left({\alpha}_{\rm max}+1\right)}$ residual.

Fig. 12 shows the successive approximations to the ratio ${\hat{\Phi}}_{\left(k_{a},k_{b}\right)}/{\hat{\Phi}}_{\left(k_{a}-1,k_{b}-1\right)}$ along the diagonal $k_{a}=k_{b}$, compared to values obtained from a Gillespie simulation. The approximation ${\alpha}_{\rm max}=0$ already shows good agreement with simulations. Divergence of the downward-going asymptotic expansion for higher-order terms begins around $k_{a}=k_{b}\approx 9$, which is the upper stable solution for the coherent-state mean number (214), as predicted in the scaling analysis of Sec. V.2.

Fig. 13 shows the error measure $\partial{\hat{\Phi}}_{\left(k_{a},k_{b}\right)}/\partial\tau$ (the deviation from a full steady-state condition) across the anti-diagonal contour $k_{a}+k_{b}=49$, testing the quality of the convergence of the ${\varphi}^{\left({\alpha}_{\rm max}\right)}$ approximation. This contour is in a range $k_{a},k_{b}\gg 9$ where the asymptotic expansions are still fairly well-controlled. The qualitative character of the convergence is well approximated along the diagonal, but exact cancellation to order $q^{2\left({\alpha}_{\rm max}+1\right)}$ degrades at higher ${\alpha}_{\rm max}$, as both the asymptotic expansion and the recursive solution (393) accumulate numerical errors. The figure also shows the crossing of error contours at $\left|q\right|/\kappa\sim 0.38$, suggesting a finite radius of convergence that does not cover the entire $\left(k_{a},k_{b}\right)$ lattice.

Fig. 14 is a contour plot of the errors $\partial{\hat{\Phi}}_{\left(k_{a},k_{b}\right)}/\partial\tau$ for ${\alpha}_{\rm max}=5$. The anti-diagonal cross section corresponds to the outermost curve from Fig. 13, raised to the 0.1 power appropriate if the total error scales as $q^{10}$. The region $k_{a}+k_{b}\lesssim 20$ shows the divergence of the asymptotic expansions already noted in Fig. 12. In the region where the asymptotic expansions converge, the errors are roughly constant along contours of fixed $\left|q\right|/\kappa$ for large $k$.