An algebraic method to calculate parameter regions for constrained steady-state distribution in stochastic reaction networks
AAn algebraic method to calculate parameter regions for constrainedsteady-state distribution in stochastic reaction networks
Tan Van Vu a) and Yoshihiko Hasegawa b) Department of Information and Communication Engineering, Graduate School of Information Science and Technology,The University of Tokyo, Tokyo 113-8656, Japan
Steady state is an essential concept in reaction networks. Its stability reflects fundamental characteristics ofseveral biological phenomena such as cellular signal transduction and gene expression. Because biochemicalreactions occur at the cellular level, they are affected by unavoidable fluctuations. Although several methodshave been proposed to detect and analyze the stability of steady states for deterministic models, these methodscannot be applied to stochastic reaction networks. In this paper, we propose an algorithm based on algebraiccomputations to calculate parameter regions for constrained steady-state distribution of stochastic reactionnetworks, in which the means and variances satisfy some given inequality constraints. To evaluate ourproposed method, we perform computer simulations for three typical chemical reactions and demonstratethat the results obtained with our method are consistent with the simulation results.
Many biological phenomena like cellular signaltransduction and gene expression can be de-scribed by stochastic reaction networks. It hasbeen known that these systems function robustlyat the steady state in the presence of noise.Therefore, it is pertinent to ask which conditionsof parameters define such stable operating regimeof the system. This kind of information providesinsights into understanding the underlying mech-anism of the system, and particularly into the de-sign processes of stochastic biocircuits. Here, wepropose an algebraic method to calculate the pa-rameter regions in which the means and variancesat the steady state satisfy some given constraints.These constraints can be, for example, an upperbound of variances, coefficients of variation, orFano factors to control fluctuations. We show inthe experiments that our approach gives resultscomparable with stochastic simulations. As allcomputations are symbolic, our method does notrequire the prior knowledge of the parameters,which are often unavailable in biological systems.The constraints can also be intentionally added toobtain the conditions of parameters under whichthe system is brought into a desired steady state.
I. INTRODUCTION
Biochemical reaction networks are mathematical mod-els used for describing biological processes, such as signaltransduction and gene expression, at the cellular level .Since most of the biological processes that constitutean organisms’ activity are either in or moving towarda steady state, elucidating the stability of the steadystate provides an understanding of the behavior of the a) Electronic mail: [email protected] b) Electronic mail: [email protected] biochemical processes . Although a great deal of ef-fort has been put into the formulation and analysis ofreaction networks, detecting and analyzing steady statesof reaction networks remain challenging . Two com-monly used models for reaction networks are continuousdeterministic models and discrete stochastic models. Indeterministic models, no randomness is involved and allstate variables are predictable. Deterministic models areoften described by a set of ordinary differential equationsand are most appropriate when the molecule numbers ofall reactant species are sufficiently large that underly-ing fluctuations can be ignored. However, in biologicalprocesses, noise is unavoidable and plays functional rolessuch as noise-induced bistability and oscillation . Lowmolecule numbers of only a few reactant species can leadto significant fluctuations. In such cases, deterministicmodels fail to accurately depict the dynamics of the sys-tem; therefore, stochastic models are necessary. Con-sequently, steady state in stochastic models becomes adistribution rather than a fixed point as in deterministicmodels. Previous studies have shown that the stochasticdynamics of a well-mixed chemically reacting system canbe accurately modeled by the chemical master equation(CME) . In most cases, the CME has not been ana-lytically solved. Hence numerical computations such asstochastic simulation algorithms , finite state projec-tion method , and quantized tensor trains are oftenconducted. Although these numerical methods can helpus obtain the distribution of the network with a par-ticular parameter value, it is intractable to apply thesemethods to calculate the parameters that yield a desiredsteady state, in which the means and variances satisfysome given constraints.Over the past few years, there have been many at-tempts to apply algebraic methods to the analysis ofreaction networks . The biggest advantage of thealgebraic method is that it does not require knowledgeof parameter values. In algebraic computations, the pa-rameters of reaction networks are treated as symbolicquantities rather than as numbers. Moreover, at thesteady state, rate equations of the deterministic model a r X i v : . [ q - b i o . M N ] J u l form an algebraic variety that can be studied using alge-braic geometry. Mart´ınez et al. proposed a procedureto locate the steady states of reaction networks from analgebraic geometry method-derived formula. By com-puting Gr¨obner basis of rate equations, one can derivethe bifurcation diagram of the reactant concentration atthe steady state in terms of a specific parameter. In ,by exploiting exact symbolic computation, an approachis presented for analyzing the stability of a large class ofbiological networks that modeled as autonomous systemsof differential equations. Siegal-Gaskins et al. appliedSturm’s theorem to the analysis of bistable biological cir-cuits. All of these works have been applied to determinis-tic models. One of the differences between stochastic anddeterministic models is the concept of the steady state.In deterministic models, the system may have multiplesteady states which are fixed points. On the contrary,in stochastic models, the system always holds a uniquesteady state depending on the initialization , whichis a probability distribution.In the present paper, we propose an algebraic methodto calculate parameter regions in which the steady-statedistribution of the reaction networks satisfies given con-straints of means, variances, or fluctuation characteristicslike coefficient of variation and Fano factor . Theprocedure of our method is as follows. First, we com-pute closed-moment equations of reaction networks viamoment closure approximations or linear noise approx-imation. At the steady state, moments are consideredto be time-invariant and moment equations form an al-gebraic variety whose solutions hold the information ofmoments at the steady state. Since moment equationsare obtained through approximation, there is the possi-bility that physically inappropriate values of momentswill be included in the solutions. To eliminate theseinadmissible solutions, inequality constraints, positivityof both means and variances, and the upper bound ofvariances are added . Eventually, a system of multi-variate polynomials that contains equations and inequa-tions is obtained. Finally, we apply an algebraic methodto compute the conditions of parameters such that thepolynomial system has exactly one solution. We demon-strate the validity of our method on three well-knownreaction network models: a gene regulatory system, two-component Michaelis–Menten enzyme reactions, and aBrusselator model. We perform stochastic simulationsto sample the desired regions of parameters. Althoughthe ranges of parameters obtained with our method areapproximate (for nonlinear systems), experiments showthat the results of our method agree with the simulation.The results of our proposed method provide an insightinto the dynamic behavior of the reaction network at thesteady state. The volume of the space of admissible pa-rameters can be considered as a quantity, which indicatesthe robustness of the system. Moreover, with the flexi-bility of adding constraint conditions, our method canbe used as a tool to explore the parameter configura-tions, which satisfy requirements in the design process of stochastic biocircuits . II. PRELIMINARIESA. Reaction networks
We consider a general reaction network with N reac-tant species X , . . . , X N interacting through M reactionchannels C , . . . , C M inside a cell with fixed volume Ω.The system is assumed to be well-mixed. The reactionchannel C j (1 ≤ j ≤ M ) is of the type: a j X + · · · + a Nj X N k j −→ b j X + · · · + b Nj X N , where a ij , b ij ∈ N ≥ are the stoichiometric coefficientsand k j ∈ R > is the macroscopic rate of reaction. Thestate of system is fully determined by the vector ofmolecule numbers of each species, n = ( n , . . . , n N ),where n i ∈ N ≥ is the molecule number of the species X i . The CME describing the time evolution of the sys-tem is given by ∂P ( n , t ) ∂t = M (cid:88) j =1 ( f j ( n − V j ) P ( n − V j , t ) − f j ( n ) P ( n , t )) , (1)where V = [ b ij − a ij ] ∈ Z N × M is a stoichiometric matrix, V j denotes the j th column of matrix V , P ( n , t ) is theprobability that the system will be in state n at time t ,and f j ( n ) represents the propensity function to accountfor the transition from a given state n to any other statein the reaction channel C j (1 ≤ j ≤ M ). Under the as-sumption of mass-action kinetics, the propensity function f j ( n ) has the following form: f j ( n ) = Ω − (cid:80) Ni =1 a ij k j N (cid:89) i =1 n i !( n i − a ij )! . (2)The solution of CME completely describes the stochasticdynamics of the system. However, in most cases, thisdifferential equation is extremely difficult to solve ex-plicitly. Consequently, stochastic simulation algorithms,such as the Gillespie algorithm or its modifications ,are often used to simulate the dynamics of the system.Although stochastic simulation can exactly describe thestochastic evolution of the system, these methods arevery computationally expensive when the number of re-actant species is large. In these cases, the precision ofthe stochastic simulation is often sacrificed for faster, yetmore approximate, methods. Various numerical and an-alytical methods have been proposed to approximatelysolve the CME, e.g., approximations of the CME so-lution by solving a truncated version of the Markovprocess , linear noise approximation , momentclosure methods , and chemical Langevin equationtreatments . B. Algebraic preliminaries
Let K be an arbitrary field and K [ x , . . . , x n ] be thering of multivariate polynomials over K with variables x ≺ x ≺ · · · ≺ x n , where ≺ denotes the ascending or-der of variables. x and x m ( m < n ) are used to denote x , . . . , x n and x , . . . , x m , respectively. For any polyno-mial P ∈ K [ x ] and variable x m , P can be viewed as aunivariate polynomial in x m over K [ x m − , x m +1 , . . . , x n ].deg( P, x m ) denotes the degree of P in x m and lc( P, x m )represents the leading coefficient of P with respect to(w.r.t.) x m . For convenience, we define deg(0 , x m ) (cid:44) − P is calledthe leading variable of P and is denoted by lv( P ). Iflv( P ) = x i then the initial and reductum of P will bedefined as follows:ini( P ) (cid:44) lc( P, x i ) , red( P ) (cid:44) P − lc( P, x i ) x deg( P,x i ) i . For any two nonzero polynomials
P, Q ∈ K [ x ] withdeg( P, x m ) = n and deg( Q, x m ) = n >
0, the pseudo-division algorithm computes two polynomials
S, R ∈ K [ x ] such that I r P = SQ + R , where I = lc( Q, x m ) , r =max( n − n + 1 , , deg( S, x m ) = max( n − n , − R, x m ) < n . The polynomials S and R are called the pseudo-quotient and pseudo-remainder of P w.r.t. Q in x m and denoted by pquo( P, Q, x m )and prem( P, Q, x m ), respectively. If lv( Q ) = x i then prem( P, Q ) (cid:44) prem( P, Q, x i ) , pquo( P, Q ) (cid:44) pquo( P, Q, x i ). gcd( P, Q, x m ) denotes the greatest com-mon divisor of P (cid:54) = 0 and Q (cid:54) = 0 w.r.t. x m .For any two polynomial sets P = { P , P , . . . , P N p } and Q = { Q , Q , . . . , Q N q } , the expressions P = 0 , Q (cid:54) =0 , Q > { P = 0 , P = 0 , . . . , P N p = 0 } , { Q (cid:54) = 0 , Q (cid:54) = 0 , . . . , Q N q (cid:54) = 0 } and { Q > , Q > , . . . , Q N q > } , respectively. Zero( P ) denotes the setof all common real zeros of the polynomials in P . Addi-tionally, we use the following notation:Zero( P \ Q ) (cid:44) { x | x ∈ Zero( P ) \ Zero( Q ) } , Zero( P , Q > (cid:44) { x | x ∈ Zero( P ) , Q ( x ) > ∀ Q ∈ Q} . For any set S , |S| denotes the cardinality of the set S ,i.e., the number of elements of S .
1. Resultant and subresultant
Suppose we are given two polynomials P ( x ) , Q ( x ) with n ≥ n > P ( x ) = p x n + p x n − + · · · + p n − x + p n ( p (cid:54) = 0) Q ( x ) = q x n + q x n − + · · · + q n − x + q n ( q (cid:54) = 0) . ( n + n ) × ( n + n ) Sylvester matrix S of P and Q hasthe following form: S = p p · · · p n · · · · · · · · · · · · p p · · · p n q q · · · q n · · · · · · · · · · · · q q · · · q n (cid:27) n (cid:27) n . Definition Resultant of two polynomials P and Q w.r.t. x , which is denoted by res( P, Q, x ), is the determi-nant of S .If P, Q ∈ K [ x ], then res( P, Q, x ) = 0 if and only if P and Q have common non-constant factors in K [ x ] .Moreover, if res( P, P (cid:48) , x ) (cid:54) = 0 then the polynomial P ( x )does not hold multiple roots.Let S i (cid:48) i be the submatrix of S obtained by deletingthe last i of the n rows of P coefficients, the last i ofthe n rows of Q coefficients and the last 2 i + 1 columns,excepting column n + n − i (cid:48) − i , for 0 ≤ i (cid:48) ≤ i ≤ n .The polynomial S i ( x ) = (cid:80) ii (cid:48) =0 det( S i (cid:48) i ) x i (cid:48) is then calledthe i th subresultant of P and Q w.r.t. x for 0 ≤ i ≤ n .If n > n + 1, the definition of i th subresultant S i ( x ) of P and Q w.r.t. x is extended as follows: S n ( x ) = q n − n − Q, S i ( x ) = 0 , n < i < n − .S i is said to be defective of degree r if deg( S i , x ) = r < i ,and regular otherwise. Definition Let
P, Q ∈ K [ x ] be two polynomials with n = deg( P, x ) ≥ deg( Q, x ) = n > n = (cid:40) n − , if n > n n , otherwise.Let S ¯ n +1 = P, S ¯ n = Q and S i be the i th subresultantof P and Q w.r.t. x for 0 ≤ i < ¯ n , then the sequenceof polynomials S ¯ n +1 , S ¯ n , . . . , S is called the subresultantchain of P and Q w.r.t. x .There is an effective method for constructing subresul-tant chains by means of pseudo-division . Definition Let S ¯ n +1 , S ¯ n , . . . , S be the subresultantchain of P and Q w.r.t. x . The sequence of regular sub-resultants S d , . . . , S d r is called the subresultant regularsubchain (s.r.s.) of P and Q w.r.t. x if1. ¯ n + 1 = d > d > · · · > d r ≥ S d i (cid:48) is regular for all 2 ≤ i (cid:48) ≤ r and S i is defectivefor all i ∈ { , . . . , ¯ n } \ { d , . . . , d r } .The s.r.s S d , . . . , S d r is renamed H , . . . , H r in Algo-rithm 5. For any two polynomials P, Q ⊂ Z [ x , . . . , x n ]with deg( P, x n ) ≥ deg( Q, x n ), the s.r.s { H i } ri =2 of P, Q provides an efficient way to calculate gcd(
P, Q, x n ) with-out computing multiple gcds. Moreover, this s.r.s can beexploited to calculate Zero( { P, Q } \ I ) as follows :Zero( { P, Q } \ I ) = r (cid:91) i =2 Zero( { H i , I i +1 , . . . , I r } \ { I, I i } ) , where I = lc( Q, x n ) , I i = lc( H i , x n ) for all i = 2 , . . . , r . Definition A polynomial set T = { T , T , . . . , T n } ⊂ K [ x ] is called a triangular set if1. T ∩ K = ∅ , i.e., T does not contain any constantpolynomial.2. lv( T i (cid:48) ) ≺ lv( T i ) for all 1 ≤ i (cid:48) < i ≤ n .We use the following notations: T ( m ) (cid:44) { T ∈ T | lv( T ) (cid:22) x m } , T (cid:104) m (cid:105) (cid:44) { T ∈ T | lv( T ) = x m } . The pseudo-remainder prem( P, T ) of any polynomial P ∈ K [ x ] w.r.t. T is defined recursively asprem( P, T ) (cid:44) prem (cid:16) prem( P, T n , lv( T n )) , T ( n − (cid:17) , where prem( P, ∅ ) (cid:44) P . Similarly, we define res( P, T ) (cid:44) res (cid:0) res( P, T n , lv( T n )) , T ( n − (cid:1) , where res( P, ∅ ) (cid:44) P . Ifres( P, T ) (cid:54) = 0 then Zero( { P } ) and Zero( T ) have no ele-ments in common, i.e., P has no common solution withthe system {T = 0 } . Definition A polynomial set T = { T , T , . . . , T n } ⊂ K [ x ] is called a regular set if1. T is a triangular set.2. res(ini( T i ) , { T , . . . , T i − } ) (cid:54) = 0 for all 1 ≤ i ≤ n ,i.e., the leading coefficient of T i is non-zero when x = ¯ x , . . . , x i − = ¯ x i − are substituted. Here, { ¯ x , . . . , ¯ x i − } is a solution of T ( i − = 0.For any two polynomial sets P and Q , if P is a regularset then we call {P , Q} a regular system. III. METHODSA. Moment equations
Because obtaining analytical probability distributionsthat satisfy the CME is very difficult in most cases, theirmoments (e.g., mean and variance) are often computedto capture the behaviors of stochastic reaction networks.For obtaining moment equations, two methods are com-monly used: moment closure approximation and linearnoise approximation (LNA). For linear systems where allpropensity functions of the system are linear, explicit ex-pressions for moments can be obtained. However, fornonlinear systems, this is not the case because each mo-ment depends on higher moments, resulting in an infi-nite hierarchy of moment equations. For that reason, several approximations have been proposed to obtainclosed moment equations. These approaches approxi-mate higher order moments by lower order ones assum-ing probability distributions (normal distribution ,lognormal distribution , Poisson distribution ) or to-tally ignore the higher order moments (central momentneglect ). If moment equations are truncated at thesecond order, i.e., only means and variances (covariances)are considered, approximations based on normal distri-bution and central moment neglect become equivalent. Ithas been reported that the moment closure approxima-tion based on normal distribution is advantageous overothers, which provides a larger range of parameter spacewhere it gives physically meaningful results . For thatreason, we adopt here the approximations based on nor-mal distribution and LNA to obtain moment equations.We stress that in our method, the approximation schemecan be replaced by other better ones if available withoutaffecting the method pipeline.
1. Approximation based on normal distribution
We define the first two moments, i.e., means µ andvariances σ , as follows: µ i = (cid:104) n i (cid:105) = (cid:88) n n i P ( n , t ) , (3) σ ii (cid:48) = (cid:104) ( n i − µ i )( n i (cid:48) − µ i (cid:48) ) (cid:105) = (cid:88) n ( n i − µ i )( n i (cid:48) − µ i (cid:48) ) P ( n , t ) , (4)and set all central moments above order two to zero. Af-ter some transformations (see Appendix A), the equa-tions of the first two moments are dµ i dt = M (cid:88) j =1 V ij (cid:32) f j ( µ ) + 12 (cid:88) h,l ∂ f j ( µ ) ∂n h ∂n l σ hl (cid:33) dσ ii (cid:48) dt = M (cid:88) j =1 (cid:32) V ij V i (cid:48) j (cid:32) f j ( µ ) + 12 (cid:88) h,l ∂ f j ( µ ) ∂n h ∂n l σ hl (cid:33) + V i (cid:48) j (cid:88) l ∂f j ( µ ) ∂n l σ il + V ij (cid:88) l ∂f j ( µ ) ∂n l σ i (cid:48) l (cid:33) . (5)At the steady state, all moments are time-invariant andmoment equations form an algebraic variety as below: M (cid:88) j =1 V ij (cid:32) f j ( µ ) + 12 (cid:88) h,l ∂ f j ( µ ) ∂n h ∂n l σ hl (cid:33) = 0 , ∀ i = 1 , . . . , N, M (cid:88) j =1 (cid:32) V ij V i (cid:48) j (cid:32) f j ( µ ) + 12 (cid:88) h,l ∂ f j ( µ ) ∂n h ∂n l σ hl (cid:33) (6)+ V i (cid:48) j (cid:88) l ∂f j ( µ ) ∂n l σ il + V ij (cid:88) l ∂f j ( µ ) ∂n l σ i (cid:48) l (cid:33) = 0 , ∀ i, i (cid:48) = 1 , . . . , N.
2. Linear noise approximation
The CME can be simplified in the LNA through Ω-expansion . To derive the LNA, we approximate theCME by a Taylor expansion of the state variables n around the mean values µ of the concentrations. Thefluctuations in n are assumed to be of the order of O (Ω / ). Then n can be expressed as n = Ω µ + √ Ω ξ , (7)where ξ are the fluctuating variables. The CME can berewritten as ∂P ( n , t ) ∂t = M (cid:88) j =1 (cid:34)(cid:32) N (cid:89) i =1 E − V ij i (cid:33) − (cid:35) f j ( n , Ω) P ( n , t ) , (8)where E vi is an operator that replaces n i by n i + v .The probability distribution P ( n , t ) can be replacedby the distribution of fluctuations Π( ξ , t ) as P ( n , t ) =Ω − N/ Π( ξ , t ). By using the Taylor expansion, we canapproximate the operator E vi as follows: E vi ≈ (cid:20) v √ Ω ∂∂ξ i + v ∂ ∂ξ i + . . . (cid:21) . (9)Substituting the approximation in Eq. (9) into Eq. (8)and collecting the terms of the order of O (Ω ) and O (Ω − / ), we obtain the rate equations for the meansand the Fokker-Planck equation for the fluctuations as d µ dt = V ν , (10) ∂ Π ∂t = − (cid:88) i,j Γ ij ∂ i ( ξ i Π) + 12 (cid:88) i,j D ij ∂ ij Π , (11)where ∂ i ≡ ∂/∂ξ i , ν j ( µ ) = lim Ω →∞ f j ( n ) / Ω andΓ ij = ∂ [ V ν ] i ∂µ j (cid:12)(cid:12)(cid:12)(cid:12) µ , D = V diag[ ν ] V (cid:62) . (12)Multiplying Eq. (11) with ξ i ξ j and integrating by partsgives the following equations of time evolution of vari-ances: d Σ dt = ΓΣ + ΣΓ (cid:62) + D , (13)where Σ ij = (cid:104) ξ i ξ j (cid:105) . We note that the means and vari-ances of molecule numbers can be calculated via µ , Σ as (cid:104) n (cid:105) = Ω µ , (cid:104) ( n − (cid:104) n (cid:105) )( n − (cid:104) n (cid:105) ) (cid:62) (cid:105) = Ω Σ . At the steadystate, the means and variances satisfy the following equa-tions: V ν = 0 , ΓΣ + ΣΓ (cid:62) + D = 0 . (14) Because the moment closure yields approximate equa-tions, the solutions of the variety defined by Eq. (6)may contain invalid steady states, including negativemeans or variances. Moreover, the variance should bebounded from above by a positive constant or a func-tion of mean to reflect the noise level. Consequently, forall 1 ≤ i ≤ N , we add into the variety the followinginequalities: µ i > , σ ii > c i − σ ii > d i µ i − σ ii > c i , d i are positive constants determined by the user. Forstochastic dynamics, an individual trajectory may un-dergo large fluctuations while its moment dynamics aretime-invariant, which cannot be regarded as a stablesteady state from a practical viewpoint. These inequal-ity constraints exclude such cases. The value d i can beinterpreted as an upper bound of the Fano factor, whichis a noise measure that correlates with the distributionwidth. One can also add the inequality d i µ i − σ ii > µ and variances σ , and parameters are macro-scopic rates of reactions k , k , . . . , k M . The problemis transformed to obtaining conditions for parameters k such that the semi-algebraic variety has exactly one realsolution. By renaming variables x = [ µ , σ ] and param-eters k = [ k , k , . . . , k M ], the semi-algebraic variety canbe viewed in the following form: (cid:40) P ( x , k ) = 0 , . . . , P N p ( x , k ) = 0 Q ( x , k ) > , . . . , Q N q ( x , k ) > . (15)Letting P = { P , . . . , P N p } , Q = { Q , . . . , Q N q } , we de-scribe semi-algebraic variety of Eq. (15) as P = 0 , Q > . (16)If a stochastic reaction network involves N reactantspecies, then the number of equations is N p = N ( N +3) /
2, and the number of inequations is N q ≥ N . B. Parameter analysis
We describe an algebraic method to analyze the condi-tions of parameters such that the semi-algebraic varietydefined by Eq. (16) has exactly one real solution.
1. Real solutions in a specific regular system
If the parameters of the regular system {P , Q} are as-signed specific values, then the distinct real solutions of {P = 0 , Q > } can be calculated precisely. Supposethat the polynomial set P has the form { P , P , . . . , P N p } ,here P i ∈ R [ x , . . . , x N p ] , ∀ i = 1 , . . . , N p . First, we com-pute distinct real roots r , . . . , r n of polynomial P ( x )and substitute each root x = r i into P to obtain a newregular set P ( i ) = { P ( i )2 , P ( i )3 , . . . , P ( i ) N p } ⊂ R [ x , . . . , x N p ].The regular set P ( i ) has N p − { x = α , . . . , x N p = α N p } is a solution of P ( i ) = 0 then { x = r i , x = α , . . . , x N p = α N p } is a solution of P = 0.By repeating this procedure, one can acquire all distinctreal solutions of P = 0. Finally, we only need to examinewhether a solution satisfies the conditions Q >
0. Thealgorithm
ExactSolve , which counts the number of dis-tinct real solutions of a regular system, is described inAlgorithm 3.
2. Regular system decomposition
A common way to analyze or solve a polynomial sys-tem is by computing the triangular decomposition of thesystem. It is known that a semi-algebraic system can bedecomposed into several regular systems . If we let P = { P , . . . , P N p } , Q = { Q , . . . , Q N q } be sets of multi-variate polynomials, then one can decompose the set P into a finite number of regular sets T , . . . , T J such thatZero( P \ Q ) = J (cid:91) i =1 Zero( T i \ Q ) . (17)In this paper, we adopt the algorithm RegSer proposedby Wang to decompose the polynomial system. Tospeed up this algorithm, we propose a parallel algorithmthat decomposes each polynomial system. A sequentialalgorithm based on RegSer is executed inside this par-allel algorithm. However, it can be observed that dur-ing the decomposition process the algorithm
RegSer produces many polynomial systems {P , Q} , which hasthe property Zero( P \ Q ) = ∅ . To reduce computa-tion and make the algorithm more efficient, we add aprobabilistic test into the sequential decomposition. Thistest probabilistically eliminates polynomial systems thathave no complex solution with some specific parametervalues. The details of each algorithm can be seen in Al-gorithm 1, 4, 5. Algorithm 1
ParallelDecomposition
Input:
A polynomial system {P , Q} in K [ x ]Output: List of regular sets Θ = [ T , . . . , T J ] set Φ ← [ {P , Q , N p } ] , Ψ ← ∅ while Φ (cid:54) = ∅ do parallel execute sequential decomposition for each1 ≤ i ≤ | Φ | [Φ ( i ) , Ψ ( i ) ] = SequentialDecomposition(Φ[ i ]) update Φ = (cid:83) i Φ ( i ) , Ψ = Ψ ∪ (cid:83) i Ψ ( i ) set Θ ← [ T , . . . , T J ], here Ψ =[ {T , U } , . . . , {T J , U J } ] return Θ To solve the parameter analysis problem, we computethe border polynomial B ( k ) , which has the followingproperty: the number of distinct real solutions of the va-riety defined by Eq. (16) is invariant over each connectedcomponent of the complement of B ( k ) = 0 in the param-eter space. First, we define the border polynomial of aregular system as follows. Definition Suppose a regular system {P , Q} is given,then the border polynomial B {P , Q} ( k ) of this system isdefined as follows: B {P , Q} ( k ) (cid:44) N p (cid:89) n =1 BP n ( k ) × N q (cid:89) m =1 BQ m ( k ) , where for all 1 ≤ n ≤ N p , ≤ m ≤ N q ,BP n ( k ) = res(res( P n , P (cid:48) n , x n ) , { P , . . . , P n − } ) ,BQ m ( k ) = res( Q m , P ) . For arbitrary semi-algebraic varieties, the border poly-nomial is defined as follows.
Definition Given a semi-algebraic variety as shownin Eq. (16), assume that the polynomial set P is decom-posed as Eq. (17), then the border polynomial B ( k ) ofthe variety will be B ( k ) (cid:44) J (cid:89) i =1 B {T i , Q} ( k ) . In the following lemma and theorem, we prove theabove-stated property of the border polynomial.
Lemma If a regular system {P , Q} satisfies that B {P , Q} ( k ) (cid:54) = 0 for all k ∈ R , where R is a continuousregion, then | Zero( P , Q > | is invariant over R . Proof. As B {P , Q} ( k ) = (cid:81) N p n =1 BP n ( k ) × (cid:81) N q m =1 BQ m ( k ) (cid:54) = 0, BP n ( k ) (cid:54) = 0 and BQ m ( k ) (cid:54) = 0 ∀ k ∈R . The conditions BP n ( k ) (cid:54) = 0 ( n = 1 , . . . , N p ) implythat the number of distinct real solutions of {P = 0 } isinvariant, while BQ m ( k ) (cid:54) = 0 ( m = 1 , . . . , N q ) indicatethat each polynomial Q m has no common solutionwith the system {P = 0 } . Let x (¯ k ) , . . . , x r (¯ k ) bethe distinct real solutions of {P = 0 } when k = ¯ k . Itcan be seen that x i (¯ k ) ( i = 1 , . . . , r ) are continuousfunctions of ¯ k over region R . Assume that indexes1 ≤ m ≤ N q , ≤ i ≤ r and parameters k (cid:54) = k ∈ R such that Q m ( x i ( k ) , k ) Q m ( x i ( k ) , k ) < Q m ( x i ( k ) , k ) is a continuous function of k , thisinequality implies that there exists ¯ k ∈ R such that Q m ( x i (¯ k ) , ¯ k ) = 0. It means that Q m has a commonsolution with the system {P = 0 } and is contradictory.Therefore, the sign of Q m ( x i ( k ) , k ) does not changeon region R for all m = 1 , . . . , N q and i = 1 , . . . , r .Consequently, | Zero( P , Q > | is invariant over region R . FIG. 1. (Color online) Outline of parameter analysis in our method. First, as in the left panel, parameter space is divided intoa finite number of subspaces such that the number of distinct real solutions of {P = 0 , Q > } is invariant over each subspace.Next, the middle panel shows that a specific value of the parameter is sampled from each subspace to calculate the number ofsolutions by exploiting the algorithm ExactSolve . Finally, all satisfied subspaces are gathered, and one obtains the conditionsof parameters as in the right panel.
Theorem Let B ( k ) be the border polynomial of thevariety defined by Eq. (16) and R be a continuous regionin the parameter space such that B ( k ) (cid:54) = 0 ∀ k ∈ R .Then the number of distinct real solutions of the varietyis invariant over R . Proof.
Since B ( k ) = (cid:81) Ji =1 B {T i , Q} ( k ) (cid:54) = 0, it is obviousthat B {T i , Q} ( k ) (cid:54) = 0 ∀ k ∈ R . According to Lemma 1,we obtain the result that | Zero( T i , Q > | is a constantover R for all i = 1 , . . . , J . Therefore, | Zero( P , Q > | = (cid:80) Ji =1 | Zero( T i , Q > | is also a constant. This meansthat the number of distinct real solutions of {P = 0 , Q > } is invariant over R .After obtaining the border polynomial B ( k ), one canuse cylindrical algebraic decomposition to decom-pose the complement of B ( k ) = 0 into finitely connectedregions such that the sign of B ( k ) does not change overeach region. The boundaries of these regions are the al-gebraic surfaces on which B ( k ) = 0 holds. Accordingto Theorem 1, the number of distinct real solutions ofthe variety is invariant over each region. By sampling anarbitrary value of the parameters from each region andapplying the algorithm ExactSolve to regular systemswith specific parameters assigned, we can easily calculatethe number of solutions of the variety in each region. Ul-timately, we gather all satisfied regions and obtain theconditions of the parameters such that the variety hasexactly one real solution. The schematic of our methodis shown in Fig. 1.Given all of these results, we propose the followingalgorithm to solve the problem of parameter analysis.
Algorithm 2
ParameterAnalysis
Input:
A polynomial system {P = 0 , Q > } Output:
The conditions of parameters k such that thesystem has exactly one real solution. set initial conditions of parameters A ← ∅ use algorithm ParallelDecomposition to decom-pose polynomial set P as Eq. (17) compute the border polynomial B ( k ) of the system {P , Q} decompose the complement of B ( k ) = 0 into finitelyconnected cells R , . . . , R n for i = 1 , . . . , n do sample an arbitrary point ¯ k i from cell R i substitute k = ¯ k i into the system {P = 0 , Q > } and use algorithm ExactSolve to count the num-ber of distinct real solutions s i (cid:44) | Zero( P , Q > | set A ← A ∪ R i if s i = 1 return A IV. RESULTS
Here, we illustrate how our proposed method works intypical stochastic reaction networks. For each case, wecompute the results with the approximations based onnormal distributions and the LNA individually. To verifythe validity of our proposed method, stochastic simula-tions are executed. For each specific parameter value, werun 10 realizations of stochastic trajectories to obtainthe means and variances at the steady state. First, weperform numerical simulations to find the boundaries ofthe regions of satisfied parameters. The boundary hereis considered to be the place at which a given constraintcondition becomes broken. After that, we uniformly sam-ple many points from the parameter space to determinethe interior of satisfied regions. By sequentially apply-ing these procedures, we obtain the simulation resultsfor each considered case. A. Gene regulatory system
We consider a simple single gene regulatory system .The system contains four reactions describing the tran-scription, translation, and degradation of the mRNA andprotein as follows (Fig. 2(a)): ∅ k (cid:29) k M, M k −→ M + P, P k −→ ∅ . (18)Here we denote the mRNA and protein by M and P ,respectively, and k i (0 ≤ i ≤
3) denote the reaction rates.The stoichiometric matrix V of the system is V = (cid:20) − − (cid:21) . Let n M ( t ) and n P ( t ) be the copy numbers of M and P ,respectively, at time t . Assuming mass-action kinetics,the propensity functions are given by f ( n M , n P ) = [Ω k , k n M , k n M , k n P ] (cid:62) . The corresponding CME is as follows: ∂P ( n M , n P , t ) ∂t = Ω k P ( n M − , n P , t )+ k ( n M + 1) P ( n M + 1 , n P , t )+ k n M P ( n M , n P − , t )+ k ( n P + 1) P ( n M , n P + 1 , t ) − (Ω k + k n M + k n M + k n P ) P ( n M , n P , t ) . (19)It has been known that depending on the parameters k i (0 ≤ i ≤ have revealed that the dy-namics of the mRNA is a Poisson process, i.e., the Fanofactor is 1, while the dynamic process of protein is super-Poissonian (the Fano factor is larger than 1). Therefore,we put here emphasis on inequation constraints of theFano factor of the protein.We denote by µ P and σ P the mean and variance, re-spectively, of the protein at the steady state. We add aconstraint σ P < µ P to set an upper bound of the Fanofactor of P . We fix parameters Ω = 100 , k = 1 and an-alyze three cases: the conditions of parameters ( k , k )when k = 10, of ( k , k ) when k = 10, and of ( k , k )when k = 10 such that the steady-state distributionsatisfies above constraint. The analytical and numericalresults of these three cases are shown in Figs. 2(b)-(d), re-spectively, in which colored regions represent conditionsobtained by the proposed method. Specifically, the bluelines and violet dashed lines express the boundaries ofthe regions obtained with approximations based on nor-mal distribution and the LNA, respectively. The circlesdenote those with numerical simulations.It can be seen that our results agree with simulation ( a ) ( b ) ( c ) NormalLNASimulation ( d ) FIG. 2. (Color online) Parameter analysis of gene regularotysystem with fixed parameters Ω = 100 , k = 1. Our analyt-ical results (blue regions) represent the conditions of param-eters such that the inequality σ P < µ P is satisfied at thesteady state. Blue lines express the boundaries of regions ob-tained with the approximation based on normal distribution.Violet dashed lines represent the counterpart obtained withthe LNA. The simulation results (orange circles) indicate theboundaries of regions of satisfied parameters. Figures fromleft to right correspond to (a) schematic diagram of a simplegene regulatory system, conditions of parameters when (b) k = 10, (c) k = 10, and (d) k = 10. results. Since the system is linear, the moment equa-tions obtained with the approximation based on normaldistribution and the LNA are same, and these equationsare not approximate but exact ones. Therefore, the re-sults obtained with both of approximation schemes areidentical with simulation. B. Michaelis–Menten enzyme reactions
The Michaelis–Menten enzyme reactions can be de-scribed by ∅ k −→ S, E + S k (cid:29) k ES k −→ E + P, (20)where E, S, ES , and P represent the free enzyme, inputsubstrate, enzyme-substrate complex, and product, re-spectively, and k i (0 ≤ i ≤
3) denotes the reaction rates.Let n E ( t ) , n S ( t ) , n ES ( t ), and n P ( t ) be the molecule num-bers of reactant species E, S, ES , and P , respectively, attime t . This reaction network has a conservation relation n E ( t ) + n ES ( t ) = n T , where n T is a constant positiveinteger. This implies that the sum of the molecule num-bers of E and ES is constant at all times. Therefore,the behavior of this system can be characterized by twovariables n E ( t ) and n S ( t ). The stoichiometric matrix V is V = (cid:20) − − (cid:21) . Assuming mass-action kinetics, the propensity functionsare given by f ( n E , n S ) = (cid:2) Ω k , Ω − k n E n S , k ( n T − n E ) , k ( n T − n E ) (cid:3) (cid:62) . (21)Substituting propensity functions in Eq. (21) into Eq. (1),we obtain the corresponding CME as follows: ∂P ( n E , n S , t ) ∂t = Ω k P ( n E , n S − , t )+ Ω − k ( n E + 1)( n S + 1) P ( n E + 1 , n S + 1 , t )+ ( n T − n E + 1) (cid:16) k P ( n E − , n S − , t )+ k P ( n E − , n S , t ) (cid:17) − (cid:16) Ω k + Ω − k n E n S + ( k + k )( n T − n E ) (cid:17) P ( n E , n S , t ) . (22)We fix parameters k = 5 and k = 4 and analyze theconditions of parameters k and k . Let µ E and σ E bethe mean and variance, respectively, of the species E atthe steady state ( µ S and σ S are defined analogously for S ). First, we add two constraints: σ E < c and σ S < c ,where c is a positive constant. The conditions of param-eters such that the means and variances satisfy the aboveconstraints are shown in Figs. 3(a)–(c), where we showthree c cases: (a) c = 200, (b) c = 300, and (c) c = 400.Next, we fix c = 500 and add two additional constraints: σ E < dµ E and σ S < dµ S , where d is a positive constant.The value of d can be considered as an upper bound of theFano factor, which is equal to one in the case of the Pois-son distribution. Here, we set several intermediate rangesfor the Fano factor and accept, to some extent, a largedispersion in the steady-state distribution. Figures 3(d)–(f) show results for three d cases with fixed c = 500: (d) d = 1 .
5, (e) d = 2 .
0, and (f) d = 2 .
5. In Fig. 3, themeanings of the colored regions, lines, and circles are thesame as in Fig. 2, and the other parameter values areshown in the caption of Fig. 3. When we relax the con-straints, the region of parameters also enlarges. Fromthe figures, for Michaelis–Menten enzyme reactions thatoften appear in biochemical reactions, we can concludethat our method gives results consistent with stochasticsimulations. Blue lines and violet dashed lines are almostidentical. This implies that the approximations based onnormal distribution and the LNA give the same accuracyin these cases. Interestingly, when we relax the constraintconditions σ E < dµ E and σ S < dµ S , the region of the ( a ) NormalLNASimulation ( b ) ( c ) ( d ) ( e ) ( f ) FIG. 3. (Color online) Parameter analysis of enzyme reac-tions with fixed parameters Ω = 100 , n T = 1000 , k = 5, and k = 4. Our analytical results (blue region) represent theregion of parameters such that inequalities σ E < c, σ S < c are satisfied at the steady state, where (a) c = 200, (b) c = 300, and (c) c = 400. The boundaries of regions ob-tained with the approximations based on normal distributionand the LNA are expressed by blue lines and violet dashedlines, respectively. Orange circles, which represent simulationresults, indicate the boundaries of regions of satisfied param-eters. When c is fixed to 500 and constraints σ E < dµ E and σ S < dµ S are added, the corresponding results are as in (d) d = 1 .
5, (e) d = 2 .
0, and (f) d = 2 . parameters undergoes a change to non-convex as shownin Figs. 3(e) and (f). The upper boundary curve indi-cates where the condition σ S < c becomes violated, i.e.,when σ S is equal to c . This yields the discontinuity ofthe parameter k if k is fixed near k = 15. C. Brusselator model
Next, we examine a nonlinear oscillating reaction net-work, the Brusselator model . This system is composed0of two reactant species X and X and the following fourreactions: ∅ k −→ X , X + X k −→ X , X k −→ X , X k −→ ∅ , (23)where k i (0 ≤ i ≤
3) denote the reaction rates. Wefix k = k = 1 and consider k and k parameters astarget of the analysis. Let n ( t ) and n ( t ) be the moleculenumbers of reactant species X and X , respectively, attime t . The stoichiometric matrix V is V = (cid:20) − − − (cid:21) . Assuming mass-action kinetics, the propensity functionsare given by f ( n , n ) = (cid:2) Ω k , Ω − k n ( n − n , k n , k n (cid:3) (cid:62) . (24)The master equation of the system is as follows: ∂P ( n , n , t ) ∂t = Ω k P ( n − , n , t )+ Ω − k ( n − n − n + 1) P ( n − , n + 1 , t )+ ( n + 1) (cid:16) k P ( n + 1 , n − , t ) + k P ( n + 1 , n , t ) (cid:17) − (cid:16) Ω k + Ω − k n ( n − n + ( k + k ) n (cid:17) × P ( n , n , t ) . (25)In the case of the deterministic model, depending on themagnitude relation of k and k + 1, the deterministicrate equations show sustained oscillations, damped oscil-lations, or overdamped oscillations. However, stochasticand deterministic models may behave qualitatively dif-ferently for some parameters. For instance, a stochasticmodel can exhibit a sustained oscillation where its corre-sponding deterministic model shows an overdamped os-cillation. Let µ X and σ X be the mean and variance,respectively, of the species X at the steady state ( µ X and σ X are defined analogously for X ). We add twoconstraints, σ X < c and σ X < c , where c is a positiveconstant, to control the noise level or amplitude of theoscillation at the steady state. We calculate the condi-tions of the parameters such that the steady-state distri-bution satisfies these constraints. The results are shownin Fig. 4 for three c cases: (a) c = 300, (b) c = 400,and (c) c = 500. Again, meanings of the colored regions,lines, and circles in Fig. 4 are the same as in Fig. 2, andthe other parameter values are shown in the caption ofFig. 4. When k is fixed, increasing k results in strongerfluctuations at the steady state. From Fig. 4, it can beseen that our method gives results comparable with thoseof stochastic simulations. When increasing the value of c , however, the results obtained with the LNA are not asgood as with the approximation based on normal distri-bution. The reason is that the LNA is derived under theassumption of small fluctuations, which are of the order ( a ) NormalLNASimulation ( b ) ( c ) FIG. 4. (Color online) Parameter analysis of the Brusselatormodel with fixed parameters Ω = 200 and k = k = 1. Blueregion represents the conditions of parameters such that in-equalities σ X < c, σ X < c are satisfied at the steady state.Blue lines and violet dashed lines express the boundaries ofregions obtained with the approximations based on normaldistribution and the LNA, respectively. The simulation re-sults (orange circles) indicate the boundaries of regions ofsatisfied parameters. Figures from left to right correspond tothe cases of (a) c = 300, (b) c = 400, and (c) c = 500. of O (Ω / ). Increasing c means that we allow more con-siderable fluctuations at the steady state, and leading topoor performance of the LNA. The LNA underestimatesthe values of variances at the steady state; therefore, theparameter regions are enlarged.In the case of the deterministic model, the system pos-sesses a stable steady state when k < k + 1. However,in the case of the stochastic model, the region of param-eters is much smaller. This can be explained by the factthat in the stochastic model, noise-induced oscillation oc-curs earlier, which means that the region of parametersis limited. V. CONCLUSIONS
Stochastic fluctuations are inevitable and ubiquitousin biological systems. Recent experimental studies haverevealed that noise plays a crucial role in the biochemicalreaction networks of living cells. For example, stochas-tic effects on gene expression lead to massive amountsof cell–cell variation observed in isogenic populations .Deterministic models, i.e., rate equations, give a macro-1scopic description of the dynamics of reaction networksand are incapable of capturing the features of the systemwhen the effects of stochastic fluctuations become signifi-cant. Thus, stochastic models, i.e., master equations, aretypically exploited to describe the dynamics of reactionnetworks as a stochastic process.Steady states play several important roles in many bi-ological functions and have been intensively studied inrecent years. However, stochastic models have not beenexplored. There is a possibility that certain steady statein the stochastic model cannot be observed with the de-terministic model . In this paper, we proposed an al-gebraic method to calculate parameter regions in whichsteady-state distribution of the reaction network satis-fies some given constraints. We examined our methodon three small reaction networks and performed numer-ical simulations to verify its validity. Through the ex-periments, it can be concluded that our method givesconsistent results with those of the simulations. Our ap-proach does not require prior knowledge of the parame-ters. This is a significant benefit since information aboutthe parameters is often unavailable in biological systems. One can also intentionally add constraints, which relatemeans and variances at steady state, to obtain the con-ditions of parameters under which the system is broughtinto a desired steady state.The precision of our method relies on moment closureapproximation. Approximations in our method (momentclosure based on normal distributions and LNA) givecomparable results in the cases of unimodal steady-statedistributions. For the networks characterized by a mul-timodal distribution, these approximations may provideunreliable results . In such cases, an approximationscheme that can handle multimodal distributions likeconditional LNA should be considered. We stress thatthe approximation of moment equations in our methodcan be flexibly replaced by other approximation schemes. ACKNOWLEDGMENTS
This work was supported by MEXT KAKENHI GrantNo.JP16K00325.
Appendix A: Derivation of moment equations
To obtain the equation of the first moment, we multiply n i by Eq. (1) and take the sum of all possible states n toget the following equation: (cid:88) n n i ∂P ( n , t ) ∂t = M (cid:88) j =1 (cid:88) n ( n i f j ( n − V j ) P ( n − V j , t ) − n i f j ( n ) P ( n , t )) . (A1)By applying the transformation n − V j → n in the first term of the right side of Eq. (A1), we obtain (cid:88) n n i ∂P ( n , t ) ∂t = M (cid:88) j =1 (cid:88) n (( n i + V ij ) f j ( n ) P ( n , t ) − n i f j ( n ) P ( n , t ))= M (cid:88) j =1 (cid:88) n V ij f j ( n ) P ( n , t ) . Thus, dµ i dt = M (cid:88) j =1 V ij (cid:104) f j ( n ) (cid:105) . (A2)Up till now, n is considered to be a vector of positive integers. To express the time derivative of the first moment byonly itself and the second central moment, we assume that n is a vector of continuous real numbers and apply theTaylor expansion for f j ( n ) around µ as follows: f j ( n ) = f j ( µ ) + ( n − µ ) (cid:62) ∂f j ( µ ) ∂ n + 12 ( n − µ ) (cid:62) ∂ f j ( µ ) ∂ n ( n − µ ) + O ( | n − µ | ) . (A3)2By approximating f j ( n ) to the order of | n − µ | and utilizing the fact that E [ n − µ ] = 0, we obtain the followingapproximation: (cid:104) f j ( n ) (cid:105) = f j ( µ ) + 12 (cid:88) h,l ∂ f j ( µ ) ∂n h ∂n l (cid:104) ( n h − µ h )( n l − µ l ) (cid:105) = f j ( µ ) + 12 (cid:88) h,l ∂ f j ( µ ) ∂n h ∂n l σ hl . Substituting the above result into Eq. (A2), we obtain a differential equation of the first moment as follows: dµ i dt = M (cid:88) j =1 V ij f j ( µ ) + 12 (cid:88) h,l ∂ f j ( µ ) ∂n h ∂n l σ hl . (A4)Similarly, to obtain the equation of the second central moment, we multiply ( n i − µ i )( n i (cid:48) − µ i (cid:48) ) by Eq. (1) and takethe sum of all possible states n to get the following equation: dσ ii (cid:48) dt = M (cid:88) j =1 (cid:88) n (( n i − µ i )( n i (cid:48) − µ i (cid:48) ) f j ( n − V j ) P ( n − V j , t ) − ( n i − µ i )( n i (cid:48) − µ i (cid:48) ) f j ( n ) P ( n , t )) . (A5)By applying the transformation n − V j → n in the first term of the right side of Eq. (A5), we obtain dσ ii (cid:48) dt = M (cid:88) j =1 (cid:88) n (( n i + V ij − µ i )( n i (cid:48) + V i (cid:48) j − µ i (cid:48) ) f j ( n ) P ( n , t ) − ( n i − µ i )( n i (cid:48) − µ i (cid:48) ) f j ( n ) P ( n , t ))= M (cid:88) j =1 (cid:88) n ( V ij V i (cid:48) j + V i (cid:48) j ( n i − µ i ) + V ij ( n i (cid:48) − µ i (cid:48) )) f j ( n ) P ( n , t ))= M (cid:88) j =1 ( V ij V i (cid:48) j (cid:104) f j ( n ) (cid:105) + V i (cid:48) j (cid:104) ( n i − µ i ) f j ( n ) (cid:105) + V ij (cid:104) ( n i (cid:48) − µ i (cid:48) ) f j ( n ) (cid:105) ) . Using the approximation of f j ( n ) in Eq. (A3) and truncating all central moments of order higher than two, we getthe following approximation: (cid:104) ( n i − µ i ) f j ( n ) (cid:105) = (cid:88) l ∂f j ( µ ) ∂n l (cid:104) ( n i − µ i )( n l − µ l ) (cid:105) = (cid:88) l ∂f j ( µ ) ∂n l σ il . Consequently, the equation of the second central moment is acquired as follows: dσ ii (cid:48) dt = M (cid:88) j =1 V ij V i (cid:48) j f j ( µ ) + 12 (cid:88) h,l ∂ f j ( µ ) ∂n h ∂n l σ hl + V i (cid:48) j (cid:88) l ∂f j ( µ ) ∂n l σ il + V ij (cid:88) l ∂f j ( µ ) ∂n l σ i (cid:48) l . Appendix B: Exact Solve algorithm
Algorithm 3
ExactSolve
Input:
A regular system {P , Q} ⊂ R [ x , . . . , x n ] Output:
The number of distinct real solutions of system {P = 0 , Q > } rs ← FindRoots( P ( x )) if rs = [ r , . . . , r h ] then pts ← [[ r ] , . . . , [ r h ]] else pts ← [ ] set i ← repeat if pts i = [ ] then pts ← [ ] and go to step 17 else let pts i = [[ r , . . . , r i ] , . . . , [ r h i , . . . , r h i i ]] set l ← substitute x = r l , . . . , x i = r li into the system {P = 0 , Q > } rs ( l ) i +1 ← FindRoots( P i +1 ( x i +1 )) if rs ( l ) i +1 = [ ] then pts ( l ) i +1 = [ ] else let rs ( l ) i +1 = [ r , . . . , r u ] pts ( l ) i +1 ← [[ r l , . . . , r li , r ] , . . . , [ r l , . . . , r li , r u ]]; l ← l + 1 if l ≤ h i then go to step 9 else pts i +1 ← (cid:83) h i i (cid:48) =1 pts ( i (cid:48) ) i +1 ; i ← i + 1 until i = n if pts = [ ] then return else return m if m members of pts make Q > Appendix C: Probabilistic Test algorithm
If a polynomial system always has at least one solution in C n , then the resultant of the system will be zero.Exploiting this property, one can calculate the resultant and check whether the system has solutions. However, ascomputational complexity of the resultant grows fast when the number of variables increases, it leads to a heavycomputation and poor time performance. Suppose that the resultant is R ( u ), where u is the parameter. Theninstead of calculating exact form of R ( u ), we compute R ( ¯ u ) (mod p ), where ¯ u is a rational value of the parameterand p is an arbitrary prime number. From a practical viewpoint, calculating R ( ¯ u ) (mod p ) is more efficient thansymbolic computation of R ( u ). If R ( u ) = 0 then we always obtain the result of zero, since R ( ¯ u ) (mod p ) = 0. Alarge value of p may lead to the high probability of the elimination; however, it also reduces the time performance.4 Algorithm 4
ProbabilisticTest
Input:
A polynomial system
P ⊂ K [ x , . . . , x n ] and a prime number p Output:
Return true if {P = 0 } probably has at least one solution in C n , else return false substitute a random rational value of the parameter u = ¯ u into P if P = { P , P , . . . , P m } , m > n then set i n ← n, i m ← m while i n > do ps ← { P , P , . . . , P i m } if ps contains nonzero constant then return false sort ps according to ascending degree of x i n if deg( ps [ i m ] , x i n ) > then h ← minimum index such that deg( ps [ h ] , x i n ) > update P i = ps [ i ] , ∀ ≤ i < h, P i m = ps [ i m ] ,P i = res( ps [ i m ] , ps [ i ] , x i n ) (mod p ) , ∀ h ≤ i < i m set i n ← i n − , i m ← i m − if P contains nonzero constant then return false return true Appendix D: Sequential Decomposition algorithm
This sequential decomposition algorithm is based on
RegSer in . Unlike RegSer which produces a list of regularsystems, this algorithm returns a list of decomposed systems and a regular system which may be empty. Severalprocesses are added to eliminate the systems that have no solutions. The prime number p is set p = 3 in Probabilis-ticTest . The pseudocode of the algorithm is shown as follows.5
Algorithm 5
SequentialDecomposition
Input:
A polynomial system {T , U} ⊂ K [ x ] and a positive integer number n Output:
Return [Φ , Ψ], where Φ is a list of decomposing systems and Ψ is a regular system set Φ ← ∅ , Ψ ← ∅ for m = n, . . . , do set T ← T \ { } , U ← U \ ( K \ { } ) if T ∩ K (cid:54) = ∅ or 0 ∈ U then go to 25 if ∃ u [ x ] ∈ U , t [ x ] ∈ T such that prem( u, t ) = 0 then go to 25 if ProbabilisticTest( T ,
3) is false then go to step 25 if T (cid:104) m (cid:105) = ∅ then go to step 21 make all polynomials in T and U to be square-free while true do let P be an element of T (cid:104) m (cid:105) with minimal degree in x m and setΦ ← Φ ∪ [ {T \ { P } ∪ { ini( P ) , red( P ) } , U , m } ] U ← U ∪ { ini( P ) } if |T (cid:104) m (cid:105) | = 1 then go to step 16 else take a polynomial P from T (cid:104) m (cid:105) \ { P } compute the s.r.s H , . . . , H r of P and P w.r.t x m set I i ← lc( H i , x m ) for 2 ≤ i ≤ r if lv( H r ) ≺ x m then set r ← r − else set r ← r set Φ ← Φ ∪ [ {T \ { P , P } ∪ { H i , I i +1 , . . . , I r } , U ∪ { I i } , m } | ≤ i ≤ r − T ← T \ { P , P } ∪ { H r , H r } , U ← U ∪ { I r } while U (cid:104) m (cid:105) (cid:54) = ∅ and lv( P ) = x m do let P be a polynomial in U (cid:104) m (cid:105) , compute s.r.s H , . . . , H r of P and P w.r.t x m set I i ← lc( H i , x m ) for 2 ≤ i ≤ r set Φ ← Φ ∪ [ {T \ { P } ∪ { pquo( P , H i , x m ) , I i +1 , . . . , I r } , U ∪ { I i } , m } | ≤ i ≤ r − T ← T \ { P } ∪ { pquo( P , H r , x m ) } , P ← pquo( P , H i , x m ) if lv( H r ) ≺ x m then set U ← U \ { P } ∪ { I r } else set U ← U ∪ { I r } if U (cid:104) m (cid:105) (cid:54) = ∅ then for all P ∈ U (cid:104) m (cid:105) do set Φ ← Φ ∪ [ {T ∪ { ini( P ) } , U \ { P } ∪ { red( P ) } , m } ] U ← U ∪ { ini( P ) } set Ψ ← {T , U} return [Φ , Ψ] M. Thattai and A. van Oudenaarden, “Intrinsic noise in generegulatory networks,” Proc. Natl. Acad. Sci. U.S.A , 8614–8619 (2001). T. Shibata and K. Fujimoto, “Noisy signal amplification in ultra-sensitive signal transduction,” Proc. Natl. Acad. Sci. U.S.A ,331–336 (2005). G. Fichera, M. A. Sneider, and J. Wyman, “On the existenceof a steady state in a biological system,” Proc. Natl. Acad. Sci.U.S.A , 4182–4184 (1977). D. Angeli, J. E. Ferrell, and E. D. Sontag, “Detection of multi- stability, bifurcations, and hysteresis in a large class of biologicalpositive-feedback systems,” Proc. Natl. Acad. Sci. U.S.A ,1822–1827 (2004). K. Gatermann, M. Eiswirth, and A. Sensse, “Toric ideals andgraph theory to analyze hopf bifurcations in mass action sys-tems,” J. Symb. Comput. , 1361 – 1382 (2005). G. Craciun, A. Dickenstein, A. Shiu, and B. Sturmfels, “Toricdynamical systems,” J. Symb. Comput. , 1551–1565 (2009). M. Mincheva, “Oscillations in non-mass action kinetics modelsof biochemical reaction networks arising from pairs of subnet- works,” J. Math. Chem. , 1111–1125 (2012). M. Domijan and M. Kirkilionis, “Bistability and oscillations inchemical reaction networks,” J. Math. Biol. , 467–501 (2009). L. S. Tsimring, “Noise in biology,” Rep. Prog. Phys. , 026601(2014). Y. Hasegawa and M. Arita, “Optimal implementations for reli-able circadian clocks,” Phys. Rev. Lett. , 108101 (2014). D. A. McQuarrie, “Stochastic approach to chemical kinetics,” J.Appl. Probab. , 413–478 (1967). D. T. Gillespie, “A rigorous derivation of the chemical masterequation,” Physica A , 404–425 (1992). D. T. Gillespie, “Exact stochastic simulation of coupled chemicalreactions,” J. Phys. Chem. , 2340–2361 (1977). M. A. Gibson and J. Bruck, “Efficient exact stochastic simulationof chemical systems with many species and many channels,” J.Phys. Chem. A , 1876–1889 (2000). D. T. Gillespie, “Approximate accelerated stochastic simulationof chemically reacting systems,” J. Chem. Phys. , 1716–1733(2001). B. Munsky and M. Khammash, “The finite state projection algo-rithm for the solution of the chemical master equation,” J. Chem.Phys. , 044104 (2006). A. Gupta, J. Mikelson, and M. Khammash, “A finite state pro-jection algorithm for the stationary solution of the chemical mas-ter equation,” J. Chem. Phys. , 154101 (2017). V. Kazeev, M. Khammash, M. Nip, and C. Schwab, “Directsolution of the chemical master equation using quantized tensortrains,” PLoS Comput. Biol. , 1–19 (2014). D. Wang and B. Xia, “Stability analysis of biological systemswith real solution classification,” in
Proceedings of the 2005 In-ternational Symposium on Symbolic and Algebraic Computation ,ISSAC ’05 (ACM, New York, 2005) pp. 354–361. A. K. Manrai and J. Gunawardena, “The geometry of multisitephosphorylation,” Biophys. J. , 5533–5543 (2008). M. Thomson and J. Gunawardena, “The rational parameteri-sation theorem for multisite post-translational modification sys-tems,” J. Theor. Biol. , 626–636 (2009). I. Mart´ınez-Forero, A. Pel´aez-L´opez, and P. Villoslada, “Steadystate detection of chemical reaction networks using a simplifiedanalytical method,” PLoS ONE , 1–6 (2010). P. M. Loriaux, G. Tesler, and A. Hoffmann, “Characterizing therelationship between steady state and response using analyticalexpressions for the steady states of mass action models,” PLoSComput. Biol. , 1–19 (2013). D. Siegal-Gaskins, E. Franco, T. Zhou, and R. M. Murray, “Ananalytical approach to bistable biological circuit discriminationusing real algebraic geometry,” J. R. Soc. Interface , 20150288(2015). I. Otero-Muras, P. Yordanov, and J. Stelling, “Chemical reactionnetwork theory elucidates sources of multistability in interferonsignaling,” PLoS Comput. Biol. , 1–28 (2017). J. M. M. Gonz´alez, “Revealing regions of multiple steady statesin heterogeneous catalytic chemical reaction networks usingGr¨obner basis,” J. Symb. Comput. , 521 – 537 (2017). D. A. Cox, J. Little, and D. O’Shea,
Ideals, Varieties, and Al-gorithms (Springer-Verlag, New York, 2007). A. Gupta, C. Briat, and M. Khammash, “A scalable computa-tional framework for establishing long-term behavior of stochas-tic reaction networks,” PLoS Comput. Biol. , 1–16 (2014). A. Gupta and M. Khammash, “Computational identification ofirreducible state-spaces for stochastic reaction networks,” SIAMJ. Appl. Dyn. Syst. , 1213–1266 (2018). E. M. Ozbudak, M. Thattai, I. Kurtser, A. D. Grossman, andA. van Oudenaarden, “Regulation of noise in the expression of asingle gene,” Nat. Genet. , 69 (2002). M. Kæn, W. J. Blake, and J. Collins, “The engineering ofgene regulatory networks,” Annu. Rev. Biomed. Eng. , 179–206(2003). J. Kuntz, P. Thomas, G.-B. Stan, and M. Barahona, “Rigor-ous bounds on the stationary distributions of the chemical mas- ter equation via mathematical programming,” arXiv:1702.05468(2017). Y. Sakurai and Y. Hori, “Optimization-based synthesis ofstochastic biocircuits with statistical specifications,” J. R. Soc.Interface (2018). J. Elf and E. M˚ans, “Fast evaluation of fluctuations in biochemi-cal networks with the linear noise approximation,” Genome Res. , 2475–2484 (2003). M. Scott, B. Ingalls, and M. Kærn, “Estimations of intrinsic andextrinsic noise in models of nonlinear genetic networks,” Chaos , 026107 (2006). N. V. Kampen,
Stochastic Processes in Physics and Chemistry (Elsevier, Amsterdam, 2007). E. W. J. Wallace, D. T. Gillespie, K. R. Sanft, and L. R. Petzold,“Linear noise approximation is valid over limited times for anychemical system that is sufficiently large,” IET Syst. Biol. ,102–115 (2012). C. A. G´omez-Uribe and G. C. Verghese, “Mass fluctuation kinet-ics: Capturing stochastic effects in systems of chemical reactionsthrough coupled mean-variance computations,” J. Chem. Phys. , 024109 (2007). L. Ferm, P. L¨otstedt, and A. Hellander, “A hierarchy of approx-imations of the master equation scaled by a size parameter,” J.Sci. Comput. , 127–151 (2008). C. H. Lee, K. H. Kim, and P. Kim, “A moment closure methodfor stochastic reaction networks,” J. Chem. Phys. , 134107(2009). C. S. Gillespie, “Moment-closure approximations for mass-actionmodels,” IET Syst. Biol. , 52–58 (2009). R. Grima, “A study of the accuracy of moment-closure approx-imations for stochastic chemical kinetics,” J. Chem. Phys. ,154105 (2012). A. Ale, P. Kirk, and M. P. H. Stumpf, “A general momentexpansion method for stochastic kinetic models,” J. Chem. Phys. , 174101 (2013). D. T. Gillespie, “The chemical Langevin equation,” J. Chem.Phys. , 297–306 (2000). B. Mishra,
Algorithmic Algebra (Springer-Verlag, New York,1993). D. Wang, “Decomposing polynomial systems into simple sys-tems,” J. Symb. Comput. , 295–314 (1998). M. J. Keeling, “Multiplicative moments and measures of persis-tence in ecology,” J. Theor. Biol. , 269 – 281 (2000). I. N˚asell, “An extension of the moment closure method,” Theor.Popul. Biol. , 233–239 (2003). D. Schnoerr, G. Sanguinetti, and R. Grima, “Comparison ofdifferent moment-closure approximations for stochastic chemicalkinetics,” J. Chem. Phys. , 185101 (2015). D. Wang, “Computing triangular systems and regular systems,”J. Symb. Comput. , 221–236 (2000). L. Yang, X. Hou, and B. Xia, “A complete algorithm for au-tomated discovering of a class of inequality-type theorems,” Sci.China Ser. F , 33–49 (2001). L. Yang and B. Xia, “Real solution classification for paramet-ric semi-algebraic systems,” in
Algorithmic Algebra and Logic (2005). D. S. Arnon, G. E. Collins, and S. McCallum, “Cylindrical alge-braic decomposition I: The basic algorithm,” SIAM J. Comput. , 865–877 (1984). D. S. Arnon, G. E. Collins, and S. McCallum, “Cylindrical alge-braic decomposition II: An adjacency algorithm for the plane,”SIAM J. Comput. , 878–889 (1984). S. K. Hortsch and A. Kremling, “Characterization of noise inmultistable genetic circuits reveals ways to modulate heterogene-ity,” PLoS ONE , 1–20 (2018). I. Prigogine and R. Lefever, “Symmetry breaking instabilities indissipative systems. II,” J. Chem. Phys. , 1695–1700 (1968). M. B. Elowitz, A. J. Levine, E. D. Siggia, and P. S. Swain,“Stochastic gene expression in a single cell,” Science , 1183–1186 (2002). P. Smadbeck and Y. N. Kaznessis, “On a theory of stability fornonlinear stochastic chemical reaction networks,” J. Chem. Phys. , 184101 (2015). D. Schnoerr, G. Sanguinetti, and R. Grima, “Validity condi- tions for moment closure approximations in stochastic chemicalkinetics,” J. Chem. Phys. , 084103 (2014). P. Thomas, N. Popovi´c, and R. Grima, “Phenotypic switchingin gene regulatory networks,” Proc. Natl. Acad. Sci. U.S.A111