An algebraic approach to product-form stationary distributions for some reaction networks
aa r X i v : . [ m a t h . P R ] D ec AN ALGEBRAIC APPROACH TO PRODUCT-FORMSTATIONARY DISTRIBUTIONS FOR SOME REACTIONNETWORKS
LINARD HOESSLY AND BEATRIZ PASCUAL-ESCUDERO
Abstract.
Exact results for product-form stationary distributions of Markovchains are of interest in different fields. In stochastic reaction networks (CRNs),stationary distributions are mostly known in special cases where they are ofproduct-form. However, there is no full characterization of the classes of net-works whose stationary distributions have product-form.We develop an algebraic approach to product-form stationary distributionsin the framework of CRNs. Under certain hypotheses on linearity and de-composition of the state space for conservative ergodic CRNs, this gives suf-ficient and necessary algebraic conditions for product-form stationary distri-butions. Correspondingly we obtain a semialgebraic subset of the parameterspace that captures rates where, under the corresponding hypotheses, CRNshave product-form. We employ the developed theory to CRNs and some mod-els of statistical mechanics, besides sketching the pertinence in other modelsfrom applied probability. Introduction
Many stochastic processes of interest count abundances of different entities, andserve as basic models in applied or theoretical probability. Often, such models canbe expressed as continuous-time Markov chains (CTMCs) on Z n ≥ . In terms ofexamples, consider, e.g., reaction networks [30], interacting particle systems [27],ecology [29], networks in the sciences [19] or Queues [35]. In this work, we study thestationary distributions of closed systems of this type, focussing on conditions forproduct-form. In particular, our results exhibit that in different cases of interestthe subset of the parameter space with product-form stationary distribution issemi-algebraic. While we focus on stochastic CRNs, the developed theory appliesbroadly to similar models, i.e. as long as the parametrisation of the CTMC modelis linear and has a certain state space decomposition (see § Reaction networks.
CRNs are an effective way to describe interactions ofdifferent species through mathematical models in the sciences. They are of interestin biochemistry, systems biology and cellular biology, besides other applications[19, 18, 15].A CRN consists of a set of reactions, that come together with associated reactionrates that govern their speed. They are usually expressed via their reaction graph.
Mathematics Subject Classification.
Key words and phrases.
Chemical reaction network, mass-action system, product-form sta-tionary distribution, Markov process, particle systems . I.e. for the parameters involved in the stochastic model in question.
As an example, consider the reversible Michaelis-Menten mechanism:(1.1) S + E ⇋ ES ⇋ P + E. Two approaches are adopted when it comes to studying the dynamics of CRNs.Deterministic models mostly consist of ordinary differential equations (ODE) whichdescribe the changes in species’ concentrations in time. If the molecular counts inthe system are low, stochastic models are usually better suited to describe thedynamics. Such systems can be analysed both in terms of the transient behavior,i.e. when it is changing in time, or the stationary behavior when the system hasreached an equilibrium, if it does. In the following we will mostly focus on ergodicstochastic CRNs and their stationary behavior.1.2.
Stochastic reaction networks and analytical results on stationary dis-tributions.
Traditionally, deterministic models have been the preferred modellingchoice. However, the analysis of living systems with small molecular counts andstochasticity has become essential. Correspondingly there is an increased interestin the dynamics of the CTMC for CRNs, albeit due to the difficulty in analyticaltractability most studies are based on approximations or simulations [15, 33, 31, 17].Distinctive differences separate the stochastic and deterministic models, both interms of their transient and stationary behavior. Besides the obvious differences,studied phenomena include, e.g., the small number effect [33], absorption in abso-lutely robust CRNs [2], condensation [23], or discreteness induced transitions [31].Characterising transient and stationary behaviour of stochastic CRNs are formi-dable tasks in general, and they are often examined via simulations [17]. Analyticalsolutions for the stationary distribution (in case of existence) are unknown for mostsystems, except for some special cases. Some stationary distributions of CRNs arewell-understood. Both complex balanced [1] and autocatalytic CRNs [23] haveproduct-form stationary distributions. Furthermore [26] offers sufficient conditionsfor product-form stationary distributions through decompositions of CRNs as wellas different examples. Note that solutions of product-form stationary distributionscan be used to prove ergodicity, see, e.g. [4, II, Corollary 4.4] for the generalcase, or [12] for the application to complex balance. Beyond these results, littleis known concerning explicit stationary distributions. Moreover, analytical resultson product-form stationary distributions allow convenient expressions for the long-term dynamics, and enable the rigorous study of scaling limits, condensation, orother features of the dynamics [23, 6, 25].While previous results were mostly focussed on particular classes of stochasticCRNs for which the stationary distributions can be given, i.e. [1, 13, 23, 26], weaim to give tools to systematically find regions of the reaction rate space where thecorresponding stochastic CRN has or has not product-form stationary distribution.The underlying idea is to notice that the condition for product-form stationary dis-tributions for closed ergodic CTMC dynamics on Z n ≥ with linear parametrisationsis canonically related to conditions expressible through algebraic geometry. Thisenables to study the semi-algebraic set of product-form for corresponding CRNscomputationally via a bottom-up approach.1.3. Overview and main results.
We recall in § § N ALGEBRAIC APPROACH TO STOCHASTIC CRN 3 (1) We show that stationary distributions have polynomial parametrisations aslong as the reaction rates are linear in the kinetics (Lemma 3.1).(2) We give necessary, and sometimes equivalent, conditions for product-formstationary distributions for sequences of probability distributions along ir-reducible components of specific forms (Theorem 3.5).(3) In § § § Notation.
We will write [ k ] for the set { , · · · , k } , and π, Z for the stationarydistributions and the normalising constants, when we do not further specify thedomain. Γ will be used for an irreducible component, and if we want to specify thedomain we write π Γ for the stationary distribution on Γ and Z Γ for the normalisingconstant.For a finite set Γ, let ∆ Γ represent the probability simplex, i.e.∆ Γ := { x ∈ R Γ | X i ∈ Γ x i = 1 , x i ≥ i ∈ Γ } , and we denote the ordinary interior consisting of probability distributions which arenonzero on all coordinates by ˚∆ Γ . Let us denote by CP Γ , RP Γ , RP Γ ≥ and RP Γ > thecomplex, real, nonnegative real and positive real projective space on Γ respectively.2. Stochastic reaction networks
Basic terminology. A reaction network consists of three finite sets G =( S , C , R ) where S is the set of species S = { S , · · · , S n } , C is the set of complexes and R ⊆ C × C is the set of reactions .Complexes correspond to nonnegative linear combinations of species over Z ≥ ,which we write in the form ν = P ni =1 ν i S i , and identify with vectors ν ∈ Z n ≥ .Reactions consist of ordered tuples ( ν, ν ′ ) ∈ R . The sets C and R are such that( ν, ν )
6∈ R and that if ν ∈ C , there exists ν ′ ∈ C such that either ( ν, ν ′ ) ∈ R , or( ν ′ , ν ) ∈ R , where we will write reactions ( ν, ν ′ ) ∈ R as ν → ν ′ . We say that areaction consumes the reactant ν and creates the product ν ′ . Each reaction ν → ν ′ has a positive (reaction) rate constant κ ν → ν ′ ; the vector of reaction rates is definedby κ ∈ R R > and the CRN with rates κ is denoted by ( G , κ ).We usually describe a reaction network by its reaction graph which is the directedgraph with vertices C and edge set R . A connected component of the reaction graphof G is termed a linkage class . We say ν ∈ C reacts to ν ′ ∈ C if ν → ν ′ is a reaction. ACRN is reversible if ν → ν ′ ∈ R whenever ν ′ → ν ∈ R . It is weakly reversible if forany reaction ν → ν ′ ∈ R , there is a sequence of directed reactions beginning with ν ′ as a reactant and ending with ν as a product complex. If it is not weakly reversiblewe say it is non-weakly reversible . The molecularity of a reaction ν → ν ′ ∈ R isequal to the number of molecules in the reactant | ν | = P i ν i . The stochiometricsubspace is defined as T = span ν ′ → ν ∈R { ν − ν ′ } ⊂ R n . The deficiency of a reaction LINARD HOESSLY AND BEATRIZ PASCUAL-ESCUDERO network G is defined as δ = |C| − ℓ − dim( T ) , where ℓ is the number of linkageclasses. A CRN G = ( S , C , R ) is conservative (resp. subconservative) if there isa vector of positive weights w ∈ R S > such that for any reaction ν → ν ′ we havethat P i ∈S ν i w i = P i ∈S ν ′ i w i (resp. P i ∈S ν i w i ≥ P i ∈S ν ′ i w i ), where ν, ν ′ ∈ N S are reactant and product.2.2. Stochastic model.
The evolution of species counts follows a Markov process,where the vector X ( t ) = x ∈ Z n ≥ changes according to transitions from the reac-tions ν → ν ′ by jumping from x to x + ν ′ − ν with transition intensity λ ν → ν ′ ( x ).The transition intensity functions λ ν → ν ′ : Z n ≥ → R ≥ then give the Q -matrix ofthe form Q ( x, x + ξ ) := X ν → ν ′ ∈R | − ν + ν ′ = ξ λ ν → ν ′ ( x )such that the Markov process satisfies P ( X ( t + ∆ t ) = x + ξ | X ( t ) = x ) = X ν → ν ′ ∈R | − ν + ν ′ = ξ λ ν → ν ′ ( x )∆ t + o (∆ t ) . The transition intensity under Mass-action kinetics (more general kinetics arepossible as well [1, 13]) associated to the reaction ν → ν ′ is(2.1) λ ν → ν ′ ( x ) = κ ν → ν ′ ( x )!( x − ν )! x ≥ ν (where z ! := n Y i =1 z i ! for z ∈ Z n ≥ ) . Stochastic CRNs are analysed via inspection of the underlying CTMC, where thestate space is decomposed into different types of states (i.e. transient, recurrent,positive recurrent, see, e.g., [32]). On irreducible components, positive recurrence(also called ergodicity) is equivalent to non-explositivity together with existence ofan invariant distribution [32]. The classification of stochastic CRNs is challengingin general. While different results exist like e.g. for positive recurrence [3], non-explositivity of complex balanced CRN [12], extinction/absorption events [24, 28],quasi-stationary distributions [28] or 1-d stochastic CRNs [37], we are far from acomplete characterization for most.We next introduce some terminology for stochastic CRNs. A reaction ν → ν ′ is active on x ∈ Z n ≥ if x ≥ ν , with ≥ the usual order structure for R n . A state u ∈ Z n ≥ is accessible from x ∈ Z n ≥ if u is reachable from x via its CTMC dynamics.We will denote this by x → G u . Definition 2.1. [Decomposition of state space] A non-empty set Γ ⊂ Z n ≥ is an irreducible component of G if for all x ∈ Γ and all u ∈ Z n ≥ , u is accessible from xif and only if u ∈ Γ.An irreducible component Γ is positive if for all reactions there is some x ∈ Γsuch that the reaction is active on that state. We say G is essential if the statespace is a union of irreducible components, and G is almost essential if the statespace is a union of irreducible components except for a finite number of states.2.3. Stationary distributions of reaction networks.
Let X ( t ) denote the un-derlying stochastic process from a CRN on an irreducible component Γ. Then,given that the stochastic process X ( t ) is ergodic and starts in Γ, we have that thelimiting distribution is the stationary distribution, i.e.,lim t →∞ P ( X ( t ) ∈ A ) = π Γ ( A ) , for any A ⊂ Γ . N ALGEBRAIC APPROACH TO STOCHASTIC CRN 5
Hence for ergodic X ( t ), the stationary distribution π Γ on an irreducible componentΓ is unique and describes the long-term behavior [32]. The stationary distributionis determined by the master equation of the Markov chain:(2.2) X ν → ν ′ ∈R π ( x + ν − ν ′ ) λ ν → ν ′ ( x + ν − ν ′ ) = π ( x ) X ν → ν ′ ∈R λ ν → ν ′ ( x ) , for all x ∈ Γ. We will denote by
M E ( G , Γ) the set of equations of the Masterequation on the irreducible component Γ.Under Mass-action kinetics, the master equation then takes the following form:(2.3) X ν → ν ′ ∈R π ( x + ν − ν ′ ) κ ν → ν ′ ( x − ν ′ + ν )!( x − ν ′ )! x ≥ ν ′ = π ( x ) X ν → ν ′ ∈R κ ν → ν ′ ( x )!( x − ν )! x ≥ ν . Remark . Given a reaction network G and assuming Mass-action kinetics, theMaster equation on an irreducible component M E ( G , Γ) is defined by (1 , κ i and stationary distributions indexed byelements of the irreducible component π Γ ( x ) , x ∈ Γ. This then defines a complexprojective algebraic set in product space X C Γ := V ( M E ( G , Γ)) ⊂ CP R × CP Γ .Solving equation (2.2) is in general a challenging task, even when restricting tothe Mass-action case. Remark . Sub-conservative CRNs have finite irreducible components, hence theirCTMC dynamics are non-explosive [32]. Note that for infinite CTMCs, existenceof stationary distribution does not imply positive recurrence, cf., e.g. [32, Ex 3.5.4].Consider a reaction network ( G , κ ) with corresponding CTMC dynamics thathas at least one nontrivial and positive recurrent irreducible component. Definition 2.4.
We say:(1) ( G , κ ) has product-form stationary distribution if there exist functions in-dexed by species f i : N → R > ; i ∈ S such that for any positive recurrentirreducible component Γ the stationary distribution has the form π Γ ( x ) = Z − Y i ∈S f i ( x i )for any x ∈ Γ, where Z Γ is a normalizing constant defined by Z Γ = P x ∈ Γ Q i ∈S f i ( x i ).(2) ( G , κ ) has I -product-form stationary distribution if there is a nontrivial covering I = ( I j ) j ∈ [ k ] of the species set S (i.e. S = ∪ j ∈ [ k ] I j ) with functions f I j : Z I j ≥ → R > such that for any positive recurrent irreducible componentΓ the corresponding stationary distribution has the form π Γ ( x ) = Z − k Y j =1 f I j ( x I j ) I.e. with more than one element. i.e S 6∈ I , and w.l.o.g. assume that I j = I l for j = l , j, l ∈ [ k ] LINARD HOESSLY AND BEATRIZ PASCUAL-ESCUDERO where Z Γ is the normalizing constant.We note that such considerations only matter on nontrivial ergodic irreduciblecomponents (consider e.g. A → B ). Note that for CRNs with infinite irreduciblecomponents, positive recurrence can be hard to check, cf., e.g., [37]. Remark . Observe that Definition 2.4 is more restrictive than, e.g., only askingto have a factorisation for one positive recurrent irreducible component.Notice that product-form corresponds to I -product-form with I all 1-elementsets as partition. For I -product-form, the following holds concerning their compar-ison. Lemma 2.6.
For two covers of S , if I ≤ I , i.e. I refines I as a cover , then I product-form implies I product-form for ( G , κ ) . Note that the notion of product-form is used in numerous works in appliedprobability [27, 19, 35, 25] and is the strongest notion, i.e. it implies any other I -product-form by Lemma 2.6. The notion of pair-product-form stationary distri-bution from statistical physics (see, e.g. [14, Section 6.1.1] or [34, § I -product-form, where I consists of pairs of overlapping 2-sets. Remark . While it is rare not to have I -product-form in studied examples, theyexist. For a mass transport model with explicit stationary distributions that issimilar, but not expressible as I -product-form we refer to [22].In the following we mostly focus on almost-essential conservative CRNs. Byfiniteness, CTMC dynamics are positive recurrent on irreducible components. Definition 2.8.
We distinguish three cases for a CRN G with given kinetics:(I) G has product-form stationary distribution independently of the rate if forany rate κ ∈ R R > , ( G , κ ) has product-form stationary distribution.(N) G has non-product-form stationary distribution independently of the rate iffor any rate κ ∈ R R > , ( G , κ ) has no product-form stationary distribution.(E) G can attain both product- and non-product-form stationary distribution ifthat holds for different ( G , κ ) with rates in κ ∈ R R > .2.4. Results on stationary distributions.
The classification of the stationarybehaviour of reaction networks is challenging, and often studied via simulations[17]. In particular, analytical solutions for the stationary distribution (in case ofexistence) are not known for most systems. Some stationary distributions of weaklyreversible reaction networks are well-understood. Complex balanced CRNs have aPoisson product-form stationary distribution [1] and can be characterized by that[7]: For a complex balanced CRN ( G , κ ) and an irreducible component Γ, thecorresponding CTMC has product-form stationary distribution of the form π Γ ( x ) = Z − c x x ! , where c ∈ R n> is a point of complex balance and Z Γ is a normalizing constant.As deficiency zero weakly reversible CRN are complex balanced independently ofthe rate, they have product-form stationary distribution independently of the rate.While it is currently not known which weakly reversible CRN admit a product-form stationary distribution beyond complex balance, at least some more have this say I ≤ I if ∀ I j ∈ I ∃ I ′ j ∈ I such that I j ⊆ I ′ j N ALGEBRAIC APPROACH TO STOCHASTIC CRN 7 property. In particular, many weakly reversible non deficiency zero CRN haveproduct-form stationary distribution independently of the rate [26].Furthermore so-called autocatalytic CRNs, a class of non-weakly-reversible CRNs,also have product-form stationary distributions, with product form functions froman infinite family of functions, where the first one specializes to the Poisson formas above. So for an autocatalytic CRN in the sense of [23, § π Γ ( x ) = Z − Y S i ∈S f i ( x i ) , with product-form functions f i ( x i ) = λ x i i x i ! x i Y l =1 (1 + n i X k =2 β ki k − Y r =1 ( l − r ))on its irreducible components ( λ i and β ki are determined by the autocatalytic CRN,cf. [23, § Z Γ the normalising constant. Results on systematical deriva-tion of sufficient conditions for product-form stationary distributions through de-compositions of CRNs with examples can be found in [26]. For a class of CRNswith so-called discreteness-induced transitions there are parameter values for whichthe form of the stationary distribution can be given [5], however they are not ofproduct form.2.5. Examples.
We end this section with some illustrative examples to introducethe reader further into the setting. Note that both reversible, weakly reversible andautocatalytic CRNs are essential, cf., e.g. [26].
Example 2.9.
Consider the following non-weakly-reversible CRN A β −→ B, A α ←− B, which is almost essential, with stochiometric compatibility classes dividing the statespace into parts Γ j = { x ∈ N | x + x = j } . The Γ j are positive irreduciblecomponents for j ≥ Example 2.10.
The following weakly-reversible deficiency one CRN2 A BA + B αγ β is almost essential with state space decomposition Γ j = { x ∈ N | x + x = j } ,where the Γ j are positive for j ≥
2. The CRN is complex balanced if and only if γ = αβ .3. Geometric view on product-form stationary distributions
In this section we first study geometric and algebraic aspects of the solutions tothe Master equation
M E ( G , Γ) on an irreducible component for a given network G . Then we give equivalent formulations for product-form stationary distributionswhich we utilize to define an ideal and a semi-algebraic set that, under certainhypotheses, characterise the reaction rates leading to product-form stationary dis-tributions. LINARD HOESSLY AND BEATRIZ PASCUAL-ESCUDERO
Expressing the solution set as a kernel.
For ergodic CTMC dynamics of a CRN on a finite irreducible component, thesolution set of the Master equation on Γ,
M E ( G , Γ) with inserted κ ∈ R R > is a1-dimensional subspace of R Γ (see § Γ ⊂ R Γ . We consider this map in a constructive way as follows, R R > h Γ −−→ R Γ > p −→ RP Γ > ≃ ∆ Γ (3.1) κ h Γ ( κ ) h Γ ( κ ) Z h Γ ( κ ) := ( π Γ ( x ) : x ∈ Γ) , and show that the map h Γ is polynomial and unique (modulo normalisation, cf.Corollary 3.2). Note that p is simply the projection into real projective space, andthe last map is the normalisation to a probability distribution (i.e. the well-knownisomorphism RP Γ ≥ ≃ ∆ Γ ). For more on basic projective algebraic geometry werefer to, e.g. [10, Section 8].Suppose that the following holds for the kinetics: Assumption 1.
The parametrisation of the kinetics λ ν → ν ′ ( · ) are linear polynomi-als in the reaction rate parameters κ ν → ν ′ for all reactions ν → ν ′ ∈ R . As an example, consider Mass-action kinetics as given in equation (2.1). Considerthe CTMC on a non-trivial irreducible component Γ. Then the system of equationsdefining
M E ( G , Γ) from equation (2.2) is given as X ν i → ν ′ i ∈R π Γ ( x + ν i − ν ′ i ) λ i ( x + ν i − ν ′ i ) = π Γ ( x ) X ν i → ν ′ i ∈R λ i ( x ) . Note that the left hand side is not zero for any x in Γ. Taking everything to theleft hand side we can write M E ( G , Γ) as X ν i → ν ′ i ∈R π Γ ( x + ν i − ν ′ i ) λ i ( x + ν i − ν ′ i ) − π Γ ( x ) X ν i → ν ′ i ∈R λ i ( x ) = 0We want to express M E ( G , Γ) in matrix form. Let Γ be the set { c , c , · · · , c m } .Then the system of equations M E ( G , Γ) can be written as:(3.2) − P ν i → ν ′ i ∈R λ i ( c ) − P ν i → ν ′ i ∈R λ i ( c m ) π Γ ( c ) π Γ ( c )... π Γ ( c m ) = 0where the entries outside of the diagonal are of the form P c k + ν i − ν ′ i = c j λ i ( c k + ν i − ν ′ i )in the k-row and j-column for j = k . Note that all λ i ( c k + ν i − ν ′ i )-entries arehomogeneous polynomials of degree 1 in κ i -variables under Assumption 1 on therate functions.We denote the corresponding matrix by A ( κ ). The following holds for A ( κ ). Lemma 3.1.
Under Assumption 1, A ( κ ) has rank m − as a matrix over the fieldof fractions R ( κ i | ν i → ν ′ i ∈ R ) .Proof. For any a ∈ R R > , the kernel of A ( a ) is 1-dimensional (by Markov chaintheory). The determinant det ( A ( κ )) is an element of R [ κ i | ν i → ν ′ i ∈ R ], which N ALGEBRAIC APPROACH TO STOCHASTIC CRN 9 vanishes for all a ∈ R R > . So it is zero in R [ κ i | ν i → ν ′ i ∈ R ].For any a ∈ R R > there is some first minor of A ( a ) which does not vanish. Hencesome first minors of A ( κ ) are nonzero, and we conclude that A ( κ ) has rank m − R ( κ i | ν i → ν ′ i ∈ R ). (cid:3) Corollary 3.2.
Under Assumption 1, the map h Γ can be given as a unique vector ( F , · · · , F m ) , up to multiplication by a constant in R ∗ , whose entries are homoge-neous polynomials in R [ κ i | ν i → ν ′ i ∈ R ] of the same degree, with positive coefficientsand no non-constant common divisor. Hence the F i never vanish in RP R > Proof.
For the ease of the reader we complete the computation of h Γ by recallingthe computation of the kernel of A ( κ ) (which is one-dimensional by Lemma 3.1)via Gaussian elimination over the field of fractions R ( κ i | ν i → ν ′ i ∈ R ). This givesa constructive proof for h Γ and ( F , · · · , F m ). We can proceed as follows:(0) (cid:20) A ( κ ) I m (cid:21) (1)&(2) ===== ⇒ (cid:20) BC (cid:21) (3) , (0) We start by creating the row-augmented matrix, where I m is the m × m identity matrix.(1) We divide each j-column by − a j,j = P ν i → ν ′ i ∈R λ i ( x j ). Hence we get − A ( κ ).(2) We apply elementary column operations to get the upper matrix in columnechelon form. Note that we only add columns multiplied by elements in R ( κ i | ν i → ν ′ i ∈ R ) with all coefficients positive.(3) We get the matrices B, C , where B is in column echelon form of degree m −
1. Furthermore matrix C has the form C = c , · · · c ,m − v ... . . . ... ... c m, · · · c m,m − v m such that the last column of C , i.e. ( v , · · · , v m ), is a non-trivial vector ofthe kernel.Hence ( v , · · · , v m ) is a basis of the kernel over R ( κ i | ν i → ν ′ i ∈ R ), and we canchoose the corresponding column vector as homogeneous polynomials. If the great-est common divisor(GCD) in R [ κ i | ν i → ν ′ i ∈ R ] is nontrivial we divide all polyno-mials of the column vector by the GCD. Denote the corresponding vector by w withcoefficients in R [ κ i | ν i → ν ′ i ∈ R ], which is unique up to multiplication by a constantin R ∗ . Then the vector w has the form w = ( F , · · · , F m ), and by construction the F , · · · , F m are homogeneous polynomials of the same degree with positive coeffi-cients by (2). As all the coefficients of the polynomials are positive, it is clear that RP R > ⊆ Z ( F , · · · , F m ) c . (cid:3) From h Γ as defined above we obtain an algebraic counterpart: the map thatsends the variables π ( x ), x ∈ Γ to h ∗ Γ ( π ( x )) = F x ( κ ) ∈ R [ κ i | ν i → ν ′ i ∈ R ] for each π ( x ), x ∈ Γ. This then defines a map between two polynomial rings: the ringof polynomials whose variables are π ( x ), x ∈ Γ and the variables κ i encoding thereaction rates, and the ring containing polynomials whose variables are only thelatter: R [ κ i | ν i → ν ′ i ∈ R ][ π ( x ) : x ∈ Γ] h ∗ Γ −→ R [ κ i | ν i → ν ′ i ∈ R ] . We further note that the map sending κ ∈ R R ( κ, h Γ ( κ )) = ( κ, F ( κ ) , · · · , F ( κ )) ∈ R R × R Γ is trivially a morphism of algebraic varieties. By Corollary 3.2, the image restrictedto R R > is contained in R R > × R Γ > , which maps into the solution of V ( M E ( G , Γ)).In particular, h ∗ Γ maps the Master equation to 0. We will use this algebraic maplater in this section.From now on we allow the following abuse of notations. We will write π forthe stationary distribution without specifying the irreducible component, π Γ whenspecifying it to Γ, π x for the value of the stationary distribution on Γ with coordinate x ∈ Γ, h Γ ( κ ) for the vector of polynomials from Corollary 3.2 on Γ and h x ( κ ) forthe corresponding polynomial of h Γ ( κ ) in coordinate x ∈ Γ. Furthermore, if thecontext is clear and we consider a sequence of irreducible components Γ i for i ∈ I along I an index set, we might write π i and h i ( κ ) for π Γ i resp. h Γ i ( κ ). Example 3.3.
Consider example 2.9 with Mass-action kinetics. The Master equa-tion is given as απ ( x − , x + 2)( x + 2)( x + 1) { x ≥ } + βπ ( x + 1 , x − x + 1) { x ≥ } = π ( x , x )( αx ( x − { x ≥ } + βx { x ≥ } )We restrict our treatment to the first three positive irreducible components Γ , Γ , Γ ,where we write the defining equations in matrix form: • M E ( G , Γ ) is defined by: − β α β − β β − α π π π = 0 • M E ( G , Γ ) is defined by: − β α β − β α β − (2 α + β ) 00 0 β − α π π π π = 0 • M E ( G , Γ ) is defined by: − β α β − β α
00 3 β − (2 α + 2 β ) 0 12 α β − (6 α + β ) 00 0 0 β − β π π π π π = 0By Lemma 3.1 the kernels A ( κ ) over R ( κ i | ν i → ν ′ i ∈ R ) are one-dimensional. Asoutlined they can be computed (e.g. using Maple, cf. § h Γ of Proposition 3.2, with positive entries in R [ κ i | ν i → ν ′ i ∈ R ]. The vector ofpolynomials from Corollary 3.2 in the irreducible components of the CTMC are thefollowing: • h Γ ( α, β ) = ( α, α, β ) • h Γ ( α, β ) = (4 α , (2 α + β )3 α, αβ, β ) • h Γ ( α, β ) = (3 α (6 α + β ) , α + 28 α β, αβ (6 α + β ) , αβ , β ) N ALGEBRAIC APPROACH TO STOCHASTIC CRN 11
Remark . Under Assumption 1, we can express (the points of) the projectivealgebraic set from Remark 2.2 set as X C Γ = V ( M E ( G , Γ)) = { ( z, w ) ∈ CP R × CP Γ | w ∈ ker( A ( z )) } and ker ( A ( z )) for z ∈ R R > corresponds to a 1-dimensional subspace of R Γ .3.2. Algebraic relations for product-form distributions.
Consider an in-dexed set of pairwise disjoint subsets Γ l of Z n ≥ with probability distributions π l in the interior of the probability simplices on these sets, π l ∈ ˚∆ Γ l . According toDefinition 2.4, (Γ l , π l ) l ∈I has product-form distribution if there are product-formfunctions f i : N → R > ; i ∈ S such that for any subset Γ l the distribution has theform π Γ l ( x ) = Z − l Y i ∈S f i ( x i ) , Z Γ l = X x ∈ Γ l Y i ∈S f i ( x i ) . The next result follows from combining Lemma 6.2 and Proposition 6.5, whichare given in § Theorem 3.5.
Consider (Γ l , π l ) l ∈I with π l ∈ ˚∆ Γ l , where Γ l = { x ∈ Z S≥ | P S i ∈S x i = l } with index set I = Z ≥ . Then the following are equivalent: (1) (Γ l , π l ) l ∈I has product-form (2) The following non-linear relations are satisfied by the distributions: (a)
For all x, x + e j − e k ∈ Γ i and y, y + e j − e k ∈ Γ i +1 , i ∈ I with x j = y j , x k = y k we have π i ( x + e j − e k ) π i +1 ( y ) = π i +1 ( y + e j − e k ) π i ( x )(b) For all x, z ∈ Γ i ; x + e j , y, z + e k , w ∈ Γ i +1 and y + e j , w + e k ∈ Γ i +2 , i ∈ I , with x j = y j and z k = w k we have π i +1 ( x + e j ) π i +1 ( y ) π i +2 ( w + e k ) π i ( z ) = π i +1 ( z + e k ) π i +1 ( w ) π i +2 ( y + e j ) π i ( x ) Moreover, if I = Z ≥ m for some m ≥ , then (1) ⇒ (2) still holds. Theorem 3.5 gives equivalent conditions for the sequence of probability distri-butions to have product-form. These conditions are expressed as relations betweenthe values of the distribution on different elements of the disjoint subsets of Z n ≥ .However, only under strict assumptions on the decomposition of Z n ≥ the conditionsfrom (2) are equivalent to product-form. For more general decompositions of thestate space (i.e. into irreducible components), analogous algebraic relations oftenexist, but they give only necessary conditions for product-form stationary distri-butions. Likewise, algebraic conditions for I -product-form stationary distributionsmight be worked out.3.3. The semi-algebraic set of product-form stationary distributions.
Inthis part we use the compatibility conditions for product-form stationary distri-butions from § h Γ i ( κ ) representing the stationary distributions from § R R > and the real algebraic set defined by the real zeros of the mentioned ideals. For basics on algebraic/semi-algebraic sets we refer to, e.g., [10].Consider a CRN G with given kinetics under Assumption 1, such that the indexedset (Γ l , h l ( κ )) l ∈I is the decomposition of the state space into irreducible componentsΓ l with unnormalised stationary distributions h l ( κ ) (parametrised by polynomials,see § l = { x ∈ Z S≥ | P S i ∈S x i = l } , with I = Z ≥ m . Further recall that the h l ( κ ) are givenby polynomials in variables κ i as proven in Proposition 3.2 of § Definition 3.6.
Consider the following ideals J j ⊂ R [ π Γ i ( x ) : x ∈ Γ i , m ≤ i ≤ j ]: J j = h homogeneous polynomials from Theorem 3.5 (a), (b) for Γ i : m ≤ i ≤ j, i , for j ∈ I ,J = h all homogeneous polynomials from Theorem 3.5 (a), (b) for all i ≥ m i We denote their images in R [ κ i | ν i → ν ′ i ∈ R ] via h ∗ Γ by I j , for j ∈ I , and I respectively. Correspondingly we have the following zero sets:- The affine algebraic set: V G := V ( I ) = { κ ∈ R R | f ( κ ) = 0 ∀ f ∈ I } - The affine semi-algebraic set: V G ,> := V ( I ) > = V ( I ) ∩ R R > Remark . Note that by construction the corresponding ideals I j , I are then sim-ply generated by the relations from Theorem 3.5 (a), (b) but with π Γ i ( x ) replacedby the polynomial in the x -coordinate of the vector from the map h Γ i ( κ ). Remark . - We have I j ⊆ I j +1 ⊆ I , hence V ( I ) ⊆ V ( I j +1 ) ⊆ V ( I j ) for all j ≥ m - R [ κ i | ν i → ν ′ i ∈ R ] is Noetherian, hence I has a finite set of generators. Equiva-lently, there exists some N ≥ m such that I m ⊆ . . . ⊆ I N = I N +1 = . . . = I .Combining Definition 3.6 with Theorem 3.5 gives the following result: Theorem 3.9.
For a conservative ergodic CRN G with irreducible components ofthe form Γ l = { x ∈ Z S≥ | P S i ∈S x i = l } for the index set I = Z ≥ under Assumption1, the following statements are equivalent for a ∈ R R > (1) ( G , a ) has product form stationary distribution (2) a ∈ V ( I ) Moreover, if I = Z ≥ m for some m ≥ , then (1) ⇒ (2) still holds.Proof. (1) = ⇒ (2): If a ∈ R R > is such that ( G , a ) has product-form stationarydistribution, any relation from Theorem 6.2 (2) is satisfied for the images of theparametrised stationary distributions via h Γ by definition. Hence a ∈ V ( I ). Notethat if I = Z ≥ m for some m ≥
2, then the implication remains valid as the corre-sponding implication from Theorem 6.2 (1) = ⇒ (2) still holds.(2) = ⇒ (1) follows similarly from the corresponding implication of Theorem6.2 (cid:3) N ALGEBRAIC APPROACH TO STOCHASTIC CRN 13
Example 3.10.
We continue examples 2.9 and 3.3. To check whether V ( I ) > isnon-empty amounts to check whether for all j ≥ V ( I j ) > is nonempty (by Remark3.8). Let us point out first that there are no relations coming from Remark 6.3 (a),so we only have to check the relations from (b). Consider the following elementof I ⊆ I from Definition 3.6 (where we are using the polynomials indexed bycoordinates of irreducible components, see example 3.3). h (0 , ( κ ) h (1 , ( κ ) h (1 , ( κ ) h (2 , ( κ ) − h (2 , ( κ ) h (1 , ( κ ) h (0 , ( κ ) h (1 , ( κ ) . Here we write h x ( κ ) for the polynomial of x -coordinate of the unnormalised sta-tionary distribution on the irreducible component (of the CRN) from § h Γ ( κ ) of example 3.3 and setting it to zerogives 6 α β [(6 α + β ) − α + β )] = 6 α β ( − β ) = 0It is easy to see that there are no positive solutions to this equation. We concludethat V ( I ) > is empty and this CRN has non-product-form stationary distributionindependently of the rate. Remark . For more general conservative or subconservative CRNs, correspond-ing algebraic relations as in Theorem 3.5 usually still exist. However, they are onlynecessary for product-form stationary distribution.4.
Examples
We consider here several classes of examples, most of which have a particu-lar state space decomposition into irreducible components. We start with CRNs,where we first study weakly reversible CRNs, and then move beyond the case ofweak reversibility. Finally we illustrate some other Markov processes from appliedprobability where the algebraic relations from Theorem 3.5 apply as well.4.1.
Reaction networks.
We consider CRNs with stochastic Mass-action kinetics that have a state spacedecomposition along the same irreducible components as in Theorem 3.5. Henceon each such irreducible component its CTMC dynamics are ergodic, i.e. positiverecurrent. We mainly aim to distinguish the three possibilities for the rate spaceconcerning product-form stationary distributions (see Definition 2.8), which weabbreviate by N = nonexistence of product-form stationary distribution ,I = existence of product-form stationary distribution independently of the rate ,E = existence of both product and non-product form stationary distribution . Remark . As a consequence of Theorem 3.5, for a CRN G with I = Z ≥ we havethe following equivalences: N ⇐⇒ V G ,> = ∅ , I ⇐⇒ V G ,> = R R > , E ⇐⇒ ∅ 6 = V G ,> ( R R >
04 LINARD HOESSLY AND BEATRIZ PASCUAL-ESCUDERO
Consider the following CRN with reactions of molecularity up to two for futurereference.(4.1) S S S S S + S αβ λ λ λ λ λ λ Examples of weakly reversible CRNs.
Let G be a weakly reversible CRN. Then we know that there exist complex balancedreaction rates for which the system has product-form stationary distribution ofPoisson form [1]. In particular, we recall that the reaction rates where the CRNis complex balanced is given by the toric variety from the moduli ideal, which wedenote here by M G as in [11]. Hence we have the following: Lemma 4.2.
Let G be weakly reversible with state space decomposition as in The-orem 3.5 and with I = Z ≥ m . Then the following relation holds between the positivepart of the toric variety V > ( M G ) from [11] and V G ,> from § ∅ 6 = V > ( M G ) ⊆ V G ,> Furthermore if the deficiency of G is zero, we have V > ( M G ) = V > ( I G ) = R R > . As there are always reaction rate values where weakly reversible CRNs are com-plex balanced, they are at least of type E or even of type I . Indeed, it is well-knownthat many weakly reversible CRNs are complex balanced independently of the rate(type I), for example if they have zero deficiency. Furthermore, it was shown in[26] that many weakly reversible CRNs of higher deficiencies are also of type I , i.e.they satisfy V > ( M G ) ( V G ,> = R R > .We next collect some examples and refer to § S −− ⇀↽ −− S S −− ⇀↽ −− S S + S −− ⇀↽ −− S W3 S −− ⇀↽ −− S S −− ⇀↽ −− S W4 2 S −− ⇀↽ −− S + S S −− ⇀↽ −− S W5 S −− ⇀↽ −− S −− ⇀↽ −− S S −− ⇀↽ −− S + S , S −− ⇀↽ −− S + S W6 S −− ⇀↽ −− S −− ⇀↽ −− S , S −− ⇀↽ −− S S −− ⇀↽ −− S + S , S −− ⇀↽ −− S + S The CRNs W1,W2, W5 and W6 all have product form stationary distributionsindependently of the rate. While the stationary distribution of a weakly reversibleCRN is projectively equivalent to a Poisson form if and only if the reaction ratevalues are complex balanced [7], different kinds of weakly reversible CRNs haveproduct-form stationary distributions beyond complex balance. On the other hand,from our investigation into the above examples we observe the following (proofs canbe found in § N ALGEBRAIC APPROACH TO STOCHASTIC CRN 15
Lemma 4.3.
W3 and W4 have product-form stationary distribution if and only ifthey are complex balanced, i.e. for W3 and W4 V > ( M G ) = V G ,> . This leads us to conjecture the following to be true.
Statement 1.
Consider a conservative and weakly reversible CRN G under Mass-action kinetics. Assume that for each r ∈ R there is no ˜ r ∈ R , r = ˜ r such that ν ′ r − ν r = ν ′ ˜ r − ν ˜ r . Then ( G , κ ) has product-form stationary distribution if and onlyif ( G , κ ) is complex balanced. In particular, assuming the positive recurrence conjecture holds (see, e.g., [3])we believe the following might be true.
Statement 2.
Consider a weakly reversible CRN G under Mass-action kineticssuch that for each r ∈ R there is no ˜ r ∈ R , r = ˜ r with ν ′ r − ν r = ν ′ ˜ r − ν ˜ r . Then ( G , κ ) has product-form stationary distribution if and only if ( G , κ ) is complex balanced,i.e. it only exhibits Poissonian product-form stationary distributions. Another natural question, assuming the positive recurrence conjecture holds,is whether the region of the reaction rate space where the CRN has product-formstationary distribution is semi-algebraic for any weakly reversible CRN under Mass-action kinetics.4.1.2.
Examples of CRNs beyond weak reversibility.
Let G be a non-weakly reversible CRN whose state space is as in Theorem 3.5. Thenwe cannot exclude any of the types I, N or E. We know that so-called autocatalyticCRNs are at least of type E or of type I [23]. Besides this, not much is currentlyknown for non-weakly reversible CRNs. We collect some examples and refer to § S −−→ S , S −−→ S S −− ⇀↽ −− S S −−→ S NW3 S −−→ S S −− ⇀↽ −− S NW4 2 S −−→ S + S S −− ⇀↽ −− S NW5 2 S ←−− S + S S −− ⇀↽ −− S NW6 S −− ⇀↽ −− S S −−→ S + S NW7 S −− ⇀↽ −− S S ←−− S + S NW8 S −− ⇀↽ −− S S ←−− S + S −−→ S NW9 S −− ⇀↽ −− S S ←−− S + S ←−− S NW10 S −− ⇀↽ −− S S −−→ S + S ←−− S
26 LINARD HOESSLY AND BEATRIZ PASCUAL-ESCUDERO
NW11 S −− ⇀↽ −− S −− ⇀↽ −− S S ←−− S + S −−→ S S ←−− S + S −−→ S NW12 S −−→ S −−→ S −−→ S S + S −−→ S , S + S −−→ S S + S −−→ S Note that NW4, NW5 have slightly different state space decomposition, i.e. theconditions from Theorem 3.5 are only necessary for product form.The so-called small number effect appears in NW7 and NW6, which correspondto [31, Motif F] and [31, Motif G]. NW1 was already briefly studied concerningits stationary distribution on some irreducible components in [7, Example 2], andit follows from our study that this is the smallest CRN with non-product-formstationary distribution independently of the rate. NW9 is an example of a simpleautocatalytic CRN on 2 species [23], while NW11 is an instance of an autocatalyticCRN on 3 species, and an inclusion process. NW12 was studied via simulationsconcerning stochastic oscillations in [36], besides being a particular instance of theinclusion process where dynamics can only move in one direction (see § § § Statement 3.
Consider a CRN G that is conservative, almost essential and non-weakly reversible under Mass-action kinetics. Then there is a connected component R ′ ⊆ R that is not weakly reversible. Assume there is r ∈ R ′ such that there is no ˜ r ∈ R \ R ′ , ˜ r = r in a weakly reversible component such that ν ′ r − ν r = ν ′ ˜ r − ν ˜ r .Then G has non-product form stationary distribution independently of the rate, i.e.it is in N. Other CTMC models.
We consider next some classes of examples of inter-est, first treating CTMCs in particle systems from statistical mechanics and thenQueuing networks. Note that we do not aim to be comprehensive in scope, insteadwe sketch cases of interest where states x ∈ Z n ≥ transition to x − e i + e j . Recall from § § Inclusion process.
The inclusion process [16, 23] is a particle system whereparticles can move from one site of a lattice to another. The CTMC dynamicsevolves in Z S≥ , where S is a finite number of sites (or set of species) which, as aCTMC, is defined by a generator L of the form L h ( x ) = X i = j p ij x i ( m x j )( h ( x + e j − e i ) − h ( x )) . N ALGEBRAIC APPROACH TO STOCHASTIC CRN 17
As observed in [23] the inclusion process can be expressed as a stochastic CRN,which is represented with sets of reactions R ij given by S i S j S i S i + S j S jα i,j α j,i α j,i α i,j The homogeneous case is a special case of the Misanthrope process on a finitelattice [9], and for p ij = p ji (where α i,j = p ij m and α i,j = p ij ) it corresponds tothe symmetric inclusion process, which is generalised by autocatalytic CRNs [23].These systems admit product-form stationary distributions for different parts ofparameter space (see, i.e., [21] or [23]), however, sufficient and necessary conditionson the rates are not known for product-form. After reparametrising inclusion pro-cesses, by Theorem 3.9 the set of rate space where the processes have product-formstationary distribution is semi-algebraic.While this is nice, in practice it is nonetheless challenging to use the approachwe derived, see NW8, NW11, NW12 for examples.4.2.2. More general particle systems.
Consider CTMC dynamics on a finite digraphof sites (i.e. G = (Λ , E ) strongly connected with | Λ | < ∞ ), where particles canjump from one site to another if the two vertices are connected in that directionunder closed but inhomogeneous dynamics. Consider the dynamics given by thegenerator L h ( x ) = X i → j b ij ( x i , x j )( h ( x + e j − e i ) − h ( x )) , where b ij are functions from Z ≥ × Z ≥ → R ≥ which satisfy(4.2) b ij ( n, m ) = 0 ⇐⇒ n = 0 , which we consider with linear parametrisation. This allows general types of func-tions b ij , allowing for inhomogeneities which can be spatial or coming from differenttypes of transition of one site to another (see, e.g. [23] or [20]). The setting aboveincorporates different well-known examples and generalisations, among which arethe following: • Zero-range process: b ij ( n, m ) = r ij f ( n ) where r ij is a parameter. • Target process: b ij ( n, m ) = n> r ij g ( n ) where r ij is a parameter. • Inclusion process: b ij ( n, m ) = n ( d + m ), where d is a parameter. Note thate.g. symmetric autocatalytic CRNs [23] with molecularity two representingthe inclusion process have functions b ij ( n, m ) = α ij n + α ij nm .Corresponding CTMCs are popular models for particle systems, and often admitproduct-form stationary distributions. We refer to [8, 21, 20] for more on the settingwe consider and to [25, 14, 34, 27] for more on context. In examples where the b ij satisfy equation (4.2), after choosing the graph structure and the constants (i.e. thefunctions f, g, v ) with linear parametrisation, the subset of parameter space where the CTMC has product-form stationary distribution is semi-algebraic (by Theorem3.9).4.2.3. Whittle networks.
In queuing networks, the particles represent customerswhich are in a location (or other entities). A fairly broad class of such CTMCsare Whittle networks [35], where we again focus on the case of closed dynamicson a finite digraph of nodes (i.e. G = ( M, E ) a strongly connected digraph with | M | < ∞ ). The generator is given by L h ( x ) = X i → j λ ij φ i ( x )( h ( x + e j − e i ) − h ( x )) , where λ ij are linear coefficients, and again reasonable conditions are usually im-posed on the functions φ i (i.e. φ i ( x ) = 0 ⇐⇒ x i = 0). Again in such a setting,it follows that, after choosing the graph and the constants (i.e. the functions φ i )with linear parametrisation, the subset of rate space with product-form stationarydistributions is semi-algebraic (by Theorem 3.9). Furthermore, more general dy-namics could be considered as well, i.e., e.g., the functions could be φ ij instead of φ i . We refer to [35, Theorem 1.15] or [4, IV § §
5] for different known sufficientconditions for product-form stationary distributions.5.
Discussion
We developed, analysed and applied algebraic approaches for product-form sta-tionary distributions to CRNs and other models from applied probability. Whilethe algebraic characterisation is only valid for CTMC dynamics with special statespace decompositions and linear parameters, it applies to many examples and givesa procedure to find the region of the reaction rate (parameter) space where thereis product-form. The approach can potentially be extended to more complicatedstate space decompositions, as long as all irreducible components are finite, or tomore general I product-form stationary distributions. Besides the applicability inCRN theory, we outlined other models of applied probability where the set of rateswith product form stationary distribution is semi-algebraic.We analysed some examples of CRNs under stochastic Mass-action kinetics, andfound the smallest CRN with no product form stationary distribution. Weaklyreversible CRNs always have some reaction rates that are complex balance, leadingto product-form stationary distribution. Therefore the set of rates giving productform is always non-empty for them. Among the examples considered here, this setis either equal to the toric variety described in [11], or equals the whole positiveorthant (i.e. of type I and with δ ≥ N ALGEBRAIC APPROACH TO STOCHASTIC CRN 19 to extend the approach developed here to more general CRNs (i.e. more generalstate space decompositions), or to systematically analyse examples for CRNs (fordifferent types of kinetics), models of statistical mechanics or stochastic networksof interest (i.e. for different configurations). In terms of CRNs it would be inter-esting to know which weakly reversible CRNs are complex balanced if and onlyif they have product form stationary distribution besides the ones with deficiencyzero. More generally, it would be interesting to know structural conditions whichimply either non-product-form stationary distributions independently of the rate,or product-form stationary distribution independently of the rate.
Acknowledgement
The authors would like to thank Jan Draisma, Christian Mazza, Manoj Gopalkr-ishnan and Elisenda Feliu for helpful discussions. LH is supported by the Swiss Na-tional Science Foundations Early Postdoc.Mobility grant (P2FRP2 188023). BPEacknowledges funding from the European Union’s Horizon 2020 research and in-novation programme under the Marie Sklodowska-Curie IF grant agreement No794627. 6.
Appendix
Some properties of product-form distribution.
In this part we take a slightly more abstract view. Consider an indexed set ofpairwise disjoint subsets Γ l ⊆ Z n ≥ with probability distributions on these sets π l ∈ ∆ Γ l , and l ∈ I which overall we denote by (Γ l , π l ) l ∈I . We say this set hasproduct-form distribution if there are product-form functions f i : N → R > ; i ∈ S such that for any subset Γ l the distribution has the form(6.1) π Γ l ( x ) = Z − l Y i ∈S f i ( x i )where Z Γ l = X x ∈ Γ l Y i ∈S f i ( x i )is the finite normalising constant. The derivations to come are consequences of thefollowing observation. Remark . Assume (Γ l , π l ) l ∈I has product form distribution. Then, if both x ∈ Γ l and x + e k ∈ Γ l ′ , it follows from (6.1) that(6.2) π Γ l ( x ) π Γ l ′ ( x + e k ) = Z Γ l ′ Z Γ l f k ( x k ) f k ( x k + 1)Our next result is the key to Theorem 6.5, and follows from the definition ofproduct form. Lemma 6.2.
Assume (Γ l , π l ) l ∈I with π l ∈ ˚∆ Γ l has product form distribution. Thenthe following two conditions are satisfied: (a) Assume that x, x + e j − e k ∈ Γ i and y, y + e j − e k ∈ Γ i +1 are such that x j = y j , x k = y k . Then π i ( x + e j − e k ) π i ( x ) = π i +1 ( y + e j − e k ) π i +1 ( y ) (b) Assume that x ∈ Γ i ; x + e j , y ∈ Γ i +1 and y + e j ∈ Γ i +2 are such that x j = y j . Then π i +1 ( x + e j ) π i ( x ) = a i π i +2 ( y + e j ) π i +1 ( y ) where a i is a constant that only depends on the index i ∈ I .In particular, these conditions are necessary for product form distribution. It is possible to express the previous conditions which were defined via fractionsas products. Together with the fact that the constants a i only depend on theirreducible component, this gives the following observation: Remark . Equivalent conditions to the ones of Lemma 6.2 that have to hold forproduct form stationary distributions are the following:(a) Assume that x, x + e j − e k ∈ Γ i and y, y + e j − e k ∈ Γ i +1 are such that x j = y j , x k = y k . Then π i ( x + e j − e k ) π i +1 ( y ) = π i +1 ( y + e j − e k ) π i ( x )(b) Assume that x, z ∈ Γ i ; x + e j , y, z + e k , w ∈ Γ i +1 and y + e j , w + e k ∈ Γ i +2 are such that x j = y j and z k = w k . Then π i +1 ( x + e j ) π i +1 ( y ) π i +2 ( w + e k ) π i ( z ) = π i +1 ( z + e k ) π i +1 ( w ) π i +2 ( y + e j ) π i ( x ) Remark . More general conditions can be given by replacing i + 1 , i + 2 byarbitrary indices i ′ , i ′′ in the formulation of Lemma 6.2 and Remark 6.3, but thenthey become less transparent. Furthermore we will only use them here in the formof Lemma 6.2. Note that, depending on the subsets Γ l of Z n ≥ , the conditions canbe empty. Proposition 6.5.
Consider (Γ l , π l ) l ∈I , where Γ l = { x ∈ Z S≥ | X S i ∈S x i = l } and π l ∈ ˚∆ Γ l . Assume that I = Z ≥ . Then the two conditions of Lemma 6.2 aresufficient for product-form.Proof. We assume n ≥
2, and instead of Lemma 6.2, we use the equivalent con-ditions in Remark 6.3. Let l = 1 and let us denote α i = π ( e i ) for i = 1 , . . . , n ,defining completely the stationary distribution on Γ . Then P ni =1 α i = 1. For thestationary distribution to be of product form we require, at this level, the existenceof functions f i : N −→ R > for i = 1 , . . . , n such that(6.3) π ( e i ) = (cid:16)Q j = i f j (0) (cid:17) f i (1) Z Γ ,where Z Γ = P x ∈ Γ Q ni =1 f i ( x i ). We can choose f i for i = 1 , . . . , n mapping 0 to1 and 1 to α i , obtaining Z Γ = P ni =1 α i = 1 such that the product-form functionssatisfy equation (6.3) for i = 1 , . . . , n .For l > f i ( l ) = Z Γ l π l ( le i ) N ALGEBRAIC APPROACH TO STOCHASTIC CRN 21 where(6.5) Z Γ l = Z Γ l − π l − ( x ) π l ( x + e k ) f k ( x k + 1) f k ( x k )for any k ∈ [ n ] and x ∈ Γ l − such that x = ( l − e k , in such a way that: (1) thedefinition of Z Γ l does not depend on the choice of x (as long as x = ( l − e k ),and (2) the f i give product form stationary distribution for Γ l (i.e. as in equation(6.1)).We proceed by 2-step induction on l : Let l = 2. Using the special structureassumed for the irreducible components, the stationary distribution is determinedby π (2 e i ) for i = 1 , . . . , n , and π ( e i + e k ) for i = k . Requiring that it has productform imposes the following conditions: π (2 e i ) = (cid:16)Q j = i f j (0) (cid:17) f i (2) Z Γ , π ( e i + e k ) = (cid:16)Q j = i,k f j (0) (cid:17) f i (1) f k (1) Z Γ .Note that for any k and for any x = e i where i = k in equation (6.5), we obtain Z Γ = α i α k π ( e i + e k ) according to the proposed definition. Let z = e j for any j = k, i .Then x, z ∈ Γ , x k , z k < x and z can both be used for (6.5)), and Remark 6.3(a), when applied to x ∈ Γ and y = x + e k ∈ Γ (where x s = y s for any s = k ),ensures π ( e j ) π ( e i + e k ) = π ( e j + e j ) π ( e i ), which implies α i α k π ( e i + e k ) = α j α k π ( e j + e k ) ,proving that Z Γ is well defined. Therefore we only need to define f i ( l ) = Z Γ π (2 e i )as in (6.4) for each i ∈ [ n ], which gives product-form stationary distribution.We assume now that we have defined f i ( l ) for all l ≤ m + 1 as in (6.4) by meansof Z Γ l and π l ( le i ), while Z Γ l itself is defined as in (6.5) in terms of Z Γ l − , π l , π l − and the f k (0) , . . . , f k ( l −
1) for k ∈ [ n ], and that (6.1) holds for these. Let us showthat then Z Γ m +2 can be defined as in (6.5) and does not depend on the choice of k and x , and that choosing then f i ( m + 2) = Z Γ m +2 π m +2 (( m + 2) e i ) for i ∈ [ n ] givesthe desired functions f , . . . , f n for the product form of π m +2 as in (6.1).Choose any w ∈ Γ m +1 , and let k be such that w k < m , so that w + e k ∈ Γ m +2 can be chosen for (6.5). Let y ∈ Γ m +1 be different than w and let j be such that y j < m , so y + e j ∈ Γ m +2 can be chosen for (6.5) too. Then there exist some x, z ∈ Γ m such that x j = y j and z k = w k . We claim that(6.6) Z Γ m +1 π m +1 ( y ) π m +2 ( y + e j ) f j ( y j + 1) f j ( y j ) = Z Γ m +1 π m +1 ( w ) π m +2 ( w + e k ) f k ( w k + 1) f k ( w k ) .The induction hypothesis ensures that Z Γ m +1 = Z Γ m π m ( x ) π m +1 ( x + e j ) f j ( x j + 1) f j ( x j ) = Z Γ m π m ( z ) π m +1 ( z + e k ) f k ( z k + 1) f k ( z k ) ,because x and z can also be chosen in (6.5), so π m ( x ) π m +1 ( z + e k ) π m +1 ( x + e j ) π m ( z ) = f k ( z k + 1) f j ( x j ) f j ( x j + 1) f k ( z k ) .But Remark 6.3 (b) makes the left hand side equal to π m +1 ( y ) π m +2 ( w + e k ) π m +2 ( y + e j ) π m +1 ( w ) and, by theway in which x and z have been chosen, the right hand side equals f k ( w k +1) f j ( y j ) f j ( y j +1) f k ( w k ) ,proving (6.6). To finish, note that the corresponding product-form functions solve the MEon Γ m +2 , because the constructed functions f i give the values of π ( x ) also for x = ( m + 2) e : this follows from the relation in Remark 6.3 b), which is guaranteedto hold between the π ( x ) for different points x of Γ m +2 . (cid:3) Remark . If |S| = 2, then the two conditions of Lemma 6.2 are sufficient forproduct-form even if I = Z ≥ . To see this, it is enough to show that the Ansatzfor product-form stationary distribution has a positive real solution if consideredon Γ , Γ , e.g. by using log ( · ) and the Rouch´e-Capelli Theorem, from which wededuce that there are infinitely many solutions. Note that this does not hold forarbitrary index sets I = Z ≥ j . Furthermore for |S| = 2 only relations of Remark6.3 (b) have to be considered.6.2. Derivations of conclusions for the examples of § For the conve-nience of the reader, we give some more detail on the analysis of the examplesconsidered in § Analysis of the semi-algebraic sets for CRNs.
Let us partially present herethe computations involving the examples in section 4.1. For this, we considerthe network in (4.1), such that all the examples on two species from § x ≥ (1 , [ π ( x , x ) αx − π ( x − , x + 1) β ( x + 1)] + x ≥ (0 , [ π ( x , x ) βx − π ( x + 1 , x − α ( x + 1)] + x ≥ (1 , [ π ( x , x ) λ x x − π ( x + 1 , x − λ ( x + 1) x ] + x ≥ (2 , [ π ( x , x )( λ + λ ) x ( x − − π ( x − , x + 2) λ ( x + 2)( x + 1) − π ( x − , x + 1) λ ( x − x + 1)] + x ≥ (0 , [ π ( x , x ) λ x ( x − − π ( x + 2 , x − λ ( x + 2)( x + 1)] = 0We will next focus on the computations for W3, W4, NW2, NW3, NW4, NW5,for which it will be sufficient to consider the corresponding ideals. For each ofthem the corresponding master equation can be obtained by setting the rates ofthe reactions that are not present to zero. Then we determine their irreduciblecomponents, and compute the kernel of the matrices A ( κ ) as in section 3.1 usingthe algorithm given there, and obtain expressions for the ideals I l defined in section3.2. For the computations we used Maple, and even though we did not performan analysis on the computational cost, for the examples considered here this wasnot an important issue, since they took less than 3 seconds on the biggest example(including kernel computation, definition of J i ideals and computation of I i for i ≤ l , l ∈ I are as in Theorem 3.5): • I = Z ≥ :W3, NW2, NW3. • I = Z ≥ :W4. • I = Z ≥ :NW4, NW5. N ALGEBRAIC APPROACH TO STOCHASTIC CRN 23
We note here that while Theorem 3.5 does not apply to NW4, NW5 as the indexsets are different, the ideals give necessary conditions for product form stationarydistribution by Lemma 6.2. Furthermore for W3, NW2, NW3 and W4 we knowby Theorem 3.5 and Remark 6.6 that the conditions in I N for N ∈ N such that I = I N , are also sufficient conditions.As a consequence of the irreducible decompositions, I is trivial for W4, NW4and NW5, and the ideal I is still trivial for NW4 and NW5, forcing us to computefurther ideals in order to obtain nontrivial necessary conditions for product form,according to Lemma 6.2.We next indicate for each of the mentioned examples the features that allows usto classify them as E/N form § J j canbe written in terms of the π ( x ), and this expression depends only on the irreduciblecomponents. Since the structure of the irreducible components is shared for allnetworks we are considering (for most of them the index set starts in 2 or 3), thereare common expressions: J = h π (0 , π (2 , π (1 , − π (1 , π (0 , π (2 , i J = J + h π (1 , π (3 , π (2 , − π (2 , π (1 , π (3 , ,π (0 , π (2 , π (1 , − π (1 , π (0 , π , ,π (0 , π (2 , π (3 , π (1 , − π (2 , π (1 , π (0 , π (3 , i J = J + h π (2 , π (4 , π (3 , − π (3 , π (2 , π (4 , ,π (1 , π (3 , π (4 , π (2 , − π (3 , π (2 , π (1 , π (4 , ,π (0 , π (3 , π (4 , π (1 , − π (3 , π (1 , π (0 , π (4 , ,π (1 , π (3 , π (2 , − π (2 , π (1 , π (3 , ,π (0 , π (2 , π (3 , π (1 , − π (2 , π (1 , π (0 , π (3 , ,π (0 , π (2 , π (1 , − π (1 , π (0 , π (2 , i After substituting each π l ( x ) by elements in the kernel of A ( κ ) for each exam-ple, i.e. h x ( κ ), we obtain the corresponding ideals I , I , I , and the followingobservations can be made: • In example NW2, the ideal I contains the polynomial β λ ( α + λ ), whichhas no solutions in the positive orthant of the rate space R > . Hence NW2is of type N. • For example NW3, I contains a polynomial of the form λ ( α +2 λ ), makingNW3 of type N. • For NW5, as we mentioned already, I = I = (0), so we need I for thefirst necessary conditions. It turns out again that I contains a polynomialwith only positive coefficients, which has also no positive solutions in therate space R . • For W3, we obtain the condition α λ − β λ = 0 already in I . Further-more we note that the above condition is the condition for complex balance,so M G ⊂ I (product form implies complex balanced) and ∅ 6 = V > ( M G ) = V G ,> by Lemma 4.2. • For W4, I is empty, but I contains the element λ λ − λ λ , which mustvanish at any set of rates potentially having product form stationary distri-bution. We note that the previous condition ensures complex balance forthe CRN. So again by Lemma 4.2, ∅ 6 = V > ( M G ) = V G ,> . • Finally, NW4 has trivial ideals I and I , but I yields some necessaryconditions in this case, which include 6 λ − λ + 7 λ = 0, λ + (2 λ − λ ) λ + λ ( λ − λ ) = 0. In the following we prove that this implies thatNW4 is of type N.W.l.o.g., let λ = 1 (the polynomials are homogeneous). Then the firstcondition gives that λ = λ +76 , the second gives 1 + (2 λ − λ ) + λ ( λ − λ ) = 0. Hence we can insert the expression for λ in the other, and weget the polynomial 24 λ + 59 λ + 36. This has no positive solution.To complete our analysis, we note that NW11, NW12 are known as particlesystems to have values with product-form stationary distribution (i.e. they are I orE), see [21, Theorem 2.1] or [23]. They both have irreducible components Γ l , l ∈ Z ≥ as in Theorem 3.5. After a check of I it is easy to find nontrivial elements of theideal showing that not all values give product-form stationary distributions (i.e. itis enough to look at, e.g. π (0 , , π (1 , , − π (1 , , π (0 , , CRNs with product-form stationary distributions independently of the rate.
We next consider the examples of type I. As solutions for stationary distributionsfor such examples were already covered in depth in [26], our presentation will berelatively brief. • Weakly reversible examples:W1 is clear, W2 is from [26, Example 5.1], and similarly W6 follows as itis a union of CRNs with compatible product-forms (see [26, Remark 5.3]). • Non-weakly reversible examples:NW6-NW10 follow as in [26, Example 5.1].We next give the product-form functions for the stationary distributions, g ( m ) = a m m ! , g ( m ) = 1 m ! m Y l =1 a + a ( l − a + a ( l − ,g ( m ) = 1 m ! m Y l =1 a a + a ( l − , g ( m ) = 1 m ! m Y l =1 a + a ( l − a , where the a i are positive constants. Then we classify the CRNs with type I accord-ing to the product-form functions in their product-form stationary distributions inthe following table. Index in S in S in S in S W1 g g − − W2 g g − − W5 g g g − W6 g g g g NW6 g g − − NW7 g g − − NW8 g g − − NW9 g g − − NW10 g g − − References [1] D. Anderson, G. Craciun, and T. Kurtz. Product-form stationary distributions for deficiencyzero chemical reaction networks.
Bul. Math. Biol. , 72:1947–1970, 2010.
N ALGEBRAIC APPROACH TO STOCHASTIC CRN 25 [2] D. Anderson, G. Enciso, and M. Johnston. Stochastic analysis of biochemical reaction net-works with absolute concentration robustness.
J. Royal Soc. Interface , 11, 2014.[3] D. Anderson and J. Kim. Some network conditions for positive recurrence of stochasticallymodeled reaction networks.
SIAM J. Appl. Math. , 78(5):2692–2713, 2018.[4] S. Asmussen and S. Asmussen.
Applied Probability and Queues . Applications of mathematics: stochastic modelling and applied probability. Springer, 2003.[5] E. Bibbona, J. Kim, and C. Wiuf. Stationary distributions of systems with discreteness-induced transitions.
Journal of The Royal Society Interface , 17(168):20200243, 2020.[6] J. Cao, P. Chleboun, and S. Grosskinsky. Dynamics of condensation in the totally asymmetricinclusion process.
Journal of Statistical Physics , 155(3):523–543, May 2014.[7] D. Cappelletti and C. Wiuf. Product-form Poisson-like distributions and complex balancedreaction systems.
SIAM J. Appl. Math. , 76(1):411–432, 2016.[8] P. Chleboun and S. Grosskinsky. Condensation in stochastic particle systems with stationaryproduct measures.
Journal of Statistical Physics , 154(1):432–465, 2014.[9] Ch. Cocozza-Thivent. Processus des misanthropes.
Zeitschrift f¨ur Wahrscheinlichkeitstheorieund Verwandte Gebiete , 70(4):509–523, Dec 1985.[10] D. A. Cox, J. Little, and D. O’Shea.
Ideals, Varieties, and Algorithms: An Introduction toComputational Algebraic Geometry and Commutative Algebra, 3/e (Undergraduate Texts inMathematics) . Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2007.[11] G. Craciun, A. Dickenstein, A. Shiu, and B. Sturmfels. Toric dynamical systems.
J. SymbolicComput. , 44(11):1551–1565, 2009.[12] M. Koyama D. F. Anderson, D. Cappelletti and T. G. Kurtz. Non-explosivity of stochasticallymodeled reaction networks that are complex balanced.
Bull. Math. Biol. , 80(10):2561–2579,2018.[13] T. D. Nguyen D. F. Anderson. Results on stochastic reaction networks with non-mass actionkinetics.
Mathematical Biosciences and Engineering , 16(mbe-16-04-103):2118, 2019.[14] M. R. Evans and T. Hanney. Nonequilibrium statistical mechanics of the zero-range processand related models.
J. Phys. A , 38(19):R195–R240, 2005.[15] C. W. Gardiner.
Handbook of stochastic methods for physics, chemistry and the naturalsciences , volume 13 of
Springer Series in Synergetics . Springer-Verlag, Berlin, third edition,2004.[16] C. Giardina, J. Kurchan, F. Redig, and K. Vafayi. Duality and hidden symmetries in inter-acting particle systems.
J. Stat. Phys. , 135:25–55, 2009.[17] D. Gillespie. Exact stochastic simulation of coupled chemical reactions.
J. Chem. Phys. ,81:2340–2361, 1977.[18] A Gorban and Gregory Y. Three waves of chemical dynamics.
Math. Model. Nat. Phenom. ,10:1–5, 08 2015.[19] J. Goutsias and G. Jenkinson. Markovian dynamics on complex reaction networks.
PhysicsReports , 529(2):199 – 264, 2013.[20] S. Grosskinsky. Stochastic and kinetic models of condensation, 2019.[21] S. Grosskinsky, F. Redig, and K. Vafayi. Condensation in the inclusion process and relatedmodels.
Journal of Statistical Physics , 142(5):952–974, Mar 2011.[22] J. Guioth and E. Bertin. A mass transport model with a simple non-factorized steady-statedistribution.
Journal of Statistical Mechanics: Theory and Experiment , 2017, 02 2017.[23] L. Hoessly and C. Mazza. Stationary distributions and condensation in autocatalytic reactionnetworks.
SIAM J. Appl. Math. , 79(4):1173–1196, 2019.[24] M. D. Johnston, D. F. Anderson, G. Craciun, and R. Brijder. Conditions for extinction eventsin chemical reaction networks with discrete state spaces.
J. Math. Biol. , 76(6):1535–1558,2018.[25] C. Kipnis and C. Landim.
Scaling Limits of Interacting Particle Systems . Springer, 1999.[26] L.Hoessly. Stationary distributions via decomposition of stochastic reaction networks. arXiv:1910.02871 , 2019.[27] T.M. Liggett.
Interacting Particle Systems . Grundlehren der mathematischen Wis-senschaften. Springer New York, 2012.[28] C. Wiuf M. C. Hansen. Existence of a unique quasi-stationary distribution for stochasticreaction networks. arxiv , 2018.[29] R. M. May. Qualitative stability in model ecosystems.
Ecology , 54(3):638–641, 1973. [30] D. A. McQuarrie. Stochastic approach to chemical kinetics.
Journal of Applied Probability ,4(3):413?478, 1967.[31] Y. Sughiyama N. Saito and K. Kaneko. Motif analysis for small-number effects in chemicalreaction dynamics.
The Journal of Chemical Physics , 145:094111, 2016.[32] J. R. Norris.
Markov Chains . Cambridge University Press, Cambridge., 1997.[33] N. Saito and K. Kaneko. Theoretical analysis of discreteness-induced transition in autocat-alytic reaction dynamics.
Phys. Rev. E , 91:022707, Feb 2015.[34] A. Schadschneider, D. Chowdhury, and K. Nishinari.
Stochastic Transport in Complex Sys-tems: From Molecules to Vehicles . 10 2010.[35] R. Serfozo.
Introduction to Stochastic Networks . Springer, New York, 1999.[36] D Spieler. Characterizing oscillatory and noisy periodic behavior in markov population mod-els. In
Quantitative Evaluation of Systems , pages 106–122, Berlin, Heidelberg, 2013. SpringerBerlin Heidelberg.[37] C. Xu, M. C. Hansen, and C. Wiuf. The asymptotic tails of limit distributions of continuoustime markov chains.
ArXiv:2007.11390 , 2020.
Department of Mathematical Sciences, University of Copenhagen, Denmark
Email address : [email protected] Department of Mathematics, Universidad Carlos III de Madrid, Spain
Email address ::