Groebner bases of reaction networks with intermediate species
GGR ¨OBNER BASES OF REACTION NETWORKS WITHINTERMEDIATE SPECIES
AMIRHOSEIN SADEGHIMANESH , ELISENDA FELIU , Abstract.
In this work we consider the computation of Gr¨obner bases of the steadystate ideal of reaction networks equipped with mass-action kinetics. Specifically, we focuson the role of intermediate species and the relation between the extended network (withintermediate species) and the core network (without intermediate species).We show that a Gr¨obner basis of the steady state ideal of the core network alwayslifts to a Gr¨obner basis of the steady state ideal of the extended network by means oflinear algebra, with a suitable choice of monomial order. As illustrated with examples,this contributes to a substantial reduction of the computation time, due mainly to thereduction in the number of variables and polynomials. We further show that if the steadystate ideal of the core network is binomial, then so is the case for the extended network,as long as an extra condition is fulfilled. For standard networks, this extra condition canbe visually explored from the network structure alone.
Keywords: binomial ideals, mass-action kinetics, steady state ideal, invariant, Gr¨obnerbasis
Introduction
Parametric polynomial systems of equations arise in the natural sciences when modelingecosystems, cell behavior, the spread of an illness, and molecular interactions within thecell, to name a few examples. In these scenarios questions of interest often boil down todescribing the solutions to these systems for varying values of the parameters. Althoughonly non-negative solutions are typically meaningful, the standard tool in computationalalgebraic geometry to study algebraic varieties, namely Gr¨obner bases, has proven useful.However, due to the parametric coefficients, the computation of a reduced Gr¨obner basiscan be time consuming for realistic examples, which typically involve many variables andparameters. The computation time depends mainly on the degree of the polynomials, thenumber of variables and coefficients, the choice of the monomial order and the used method[1, 2, 5, 9, 29]. These universal considerations target generic polynomial systems, but, inapplications, the structure of the particular system of interest might favor one method orone monomial order over another.We focus on a specific type of polynomial systems that arise when modeling chemicalreaction networks with mass-action kinetics [10, 12]. Specifically, the evolution of theconcentrations of the species of a chemical reaction network in time is described undermass-action by a system of ordinary differential equations in R ndx i dt = f κ,i ( x ) , i = 1 , . . . , n with f κ,i ( x ) polynomial. The monomials of each f κ,i ( x ) depend on the reaction networkstructure alone, and the coefficients depend on the reaction rate constants κ , which areoften unknown and thus treated as parameters. The steady states, or equilibrium points,of the system are the non-negative points of the variety defined by the steady state ideal I κ = (cid:104) f κ, ( x ) , . . . , f κ,n ( x ) (cid:105) . Department of Mathematical Sciences, University of Copenhagen, Universitetsparken 5, 2100 Copen-hagen, Denmark. Corresponding author: [email protected]
Date : April 5, 2018. a r X i v : . [ c s . S C ] M a r R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 2
The question of restricting to non-negative steady states remains challenging and nostraightforward solutions have been proposed. Despite of this, Gr¨obner bases have beenfor example used for model discrimination [13, 14, 18–20]. They are also used to decidewhether the steady state ideal is binomial, that is, whether any reduced Gr¨obner basisconsists of polynomials with at most two terms. If this is the case, then methods to detectthe existence of multiple steady states can be applied [22, 25].In this work we exploit the specific structure of the steady state ideal, which reflectsthe structure of the reaction network, to guide the selection of good monomial orders andto compute reduced Gr¨obner bases faster. Specifically, we consider a frequent and nicely-behaved class of species introduced in [11] called intermediate species (or intermediates, forshort). Intermediates give rise to linear terms in the steady state polynomials, and theycan be removed from a reaction network resulting in a smaller core network with only thenon-intermediates. A key property is that steady states of the core network can be liftedto steady states of the extended network .The first main result of this work is Theorem 3.4, where we show how to obtain a Gr¨obnerbasis of the extended network from one of the core network using linear algebra. The resultimplicitly gives good monomial orders, namely, those for which the concentration of the in-termediates are larger than for the non-intermediates, and are lexicographic in the variablescorresponding to the intermediates. Example 3.5 illustrates the computational advantageof using our approach compared with other methods. Additionally, we conclude that theanalysis of the steady state ideal of the core network is sufficient for model discrimination.The second main result, Theorem 3.10, addresses how to decide whether the steady stateideal is binomial. We show that if the steady state ideal of the core network is binomial,then this is also the case for the steady state ideal of the extended network provided anextra condition is fulfilled. In typical networks, this extra condition can be readily checkedfrom the network structure alone. When the core network has a homogeneous steady stateideal (which happens frequently for realistic reaction networks), then one can employ thelinear algebra-based method introduced in [4] to detect whether the steady state ideal ofthe core network is binomial. Then, combined with our result, we obtain a faster methodto address whether the steady state ideal of the original network is binomial, which doesnot rely on the computation of a Gr¨obner basis.The key property behind our results is that intermediates define a square linear subsys-tem of full rank among the steady state polynomials. Its solution and posterior substitutioninto the remaining polynomials gives rise to a smaller ideal in the non-intermediates. AGr¨obner basis of the small ideal can then be lifted to a Gr¨obner basis of the original ideal.Our approach can be theoretically applied to arbitrary parametric ideals, after detectionof linear subsystems among a set of generators. However, technical conditions that arenecessary for our results to hold might not be straightforward to check, since we overcomethis difficulty by exploiting the network structure.The structure of the paper is as follows. We start by introducing reaction networks andbasic concepts such as the steady state ideal. Intermediates are introduced in Section 2.In Section 3 we address Gr¨obner bases of networks with intermediates, discuss binomialsteady state ideals and relate our work to [4]. In Section 4 a technical condition of algebraicindependence of a set of rational functions, which is assumed in the former sections, isdiscussed. Finally, in the last section, we discuss another class of special species, namelyenzymes, that might lead to similar results concerning the computation of Gr¨obner bases.1.
The steady state ideal of a reaction network
We follow the formalism of [11]. See also [10, 12] for an introduction to reaction networks.Subscripts ≥ , > R (resp. Z ) refer to the non-negative and positive real numbers(resp. integer numbers). R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 3 A reaction network is an ordered triple N = ( S , C , R ) where S , C and R are three setscalled the set of species , complexes and reactions , respectively. Here S is a finite set and C is a finite set of linear combinations of elements of S with coefficients in Z ≥ . A reaction is an ordered pair of complexes ( c, c (cid:48) ) in C , usually denoted as c → c (cid:48) . For the reaction c → c (cid:48) , the complex c is called the reactant and c (cid:48) is called the product .A digraph is associated with a reaction network as follows. The vertex set is C and thereis a directed edge from the reactant to the product of every reaction. If both reactions c → c (cid:48) and c (cid:48) → c for two complexes c and c (cid:48) exist, then the notation c (cid:10) c (cid:48) is used andthe reaction is said to be reversible.Complexes that are not part of any reaction or species that are not part of any complex donot appear in the digraph. Therefore, the reaction network cannot uniquely be determinedfrom the digraph alone. For simplicity, however, we often introduce a reaction network byits digraph and implicitly assume that the set of complexes equals the set of vertices andthe set of species consists of the species that appear in at least one complex.Write S = { X , . . . , X n } , such that the set of species is implicitly ordered. Then acomplex c is of the form c X + · · · + c n X n , which we also write in vector form as c =( c , c , . . . , c n ) ∈ Z n ≥ . With this representation, c i is called the stoichiometric coefficient of X i in c . Example 1.1.
Let S = { X , X , X , X } , C = { X + X , X , X + X } , R = { X + X → X , X → X + X , X → X + X } . The network N = ( S , C , R ) is represented with thefollowing digraph X + X −− (cid:42)(cid:41) −− X −−→ X + X . The complexes X + X and X appear both as reactants and products while X + X appears only as a product.We next construct a system of Ordinary Differential Equations (ODEs) that models thevariation of the concentration of each species in time and introduce the relevant polynomials F i ( x ) that are the focus of this work. We denote the concentration of each species X i inlower-case x i . For each reaction c → c (cid:48) , we introduce a parameter k c → c (cid:48) , and a polynomial F i ( x ) is associated with every species X i as follows: F i ( x ) = (cid:88) c → c (cid:48) ∈R ( c (cid:48) i − c i ) k c → c (cid:48) x c ∈ R ( k )[ x ] , (1)where x c = x c . . . x c n n . Here x = ( x , . . . , x n ) and R ( k ) is the field of rational functionswith variables k c → c (cid:48) and real coefficients. The symbol k stands for the parameter vector k = ( k c → c (cid:48) | c → c (cid:48) ∈ R ) . For a chosen positive value k (cid:63) ∈ R R > of the parameter vector, we let F k (cid:63) ,i ( x ) ∈ R [ x ]denote the image of F i ( x ) under the evaluation map R ( k ) → R , k c → c (cid:48) (cid:55)→ k (cid:63)c → c (cid:48) . With this choice of k (cid:63) , the ODE system of the reaction network under mass-actionkinetics is ˙ x i = F k (cid:63) ,i ( x ) , i = 1 , . . . , n, x ∈ R n ≥ . (2)The value k (cid:63)c → c (cid:48) > reaction rate constant of c → c (cid:48) and is usually depictedas a label of the reaction in the associated digraph. By [27], if the starting condition of(2) belongs to R n> (resp. R n ≥ ), then so does the trajectory for all positive times in theinterval of definition.The steady states of the network are the common zeros of F k (cid:63) ,i ( x ), i = 1 , . . . , n . Inapplications, only non-negative real solutions have meaning and mostly, positive steadystates are interesting, meaning all concentrations are positive. Since the values of thereaction rate constants are in general unknown, they are treated as parameters of the R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 4 system. Thus we aim at studying the zeros of the system of polynomials F i ( x ) = 0, for i = 1 , . . . , n in R ( k ) and specially the positive zeros after specifying values for k . Definition 1.2.
Let N = ( S , C , R ) be a reaction network with S = { X , . . . , X n } .(a) F i ( x ) ∈ R ( k )[ x ] is called the steady state polynomial of X i .(b) The ideal generated by the steady state polynomials of all the species in the networkin the ring R ( k )[ x ] is called the steady state ideal of the network: I N = (cid:10) F i ( x ) | i = 1 , . . . , n (cid:11) ⊆ R ( k )[ x ] . The set of steady states for a vector of reaction rate constants k (cid:63) is thus the solution setto any basis (set of generators) of I N specialized to k (cid:63) .It follows from (1) and (2) that for all x ∈ R n , the vector F k ( x ) = ( F k, ( x ) , . . . , F k,n ( x ))lies in the vector subspace S = span( { c − c (cid:48) | c → c (cid:48) ∈ R} ) ⊆ R n . If s = dim( S ), then n − s of the steady state polynomials can be written as linear combinations of the remaining s polynomials. We conclude that it is always possible to find a basis of I N with cardinalitydim( S ). Example 1.3. (continued from Example 1.1) The ODE system of the reaction networkwith digraph X + X k −− (cid:42)(cid:41) −− k X k −−→ X + X is ˙ x = − k x x + k x ˙ x = k x ˙ x = − k x x + k x + k x ˙ x = k x x − k x − k x . In this case dim( S ) = 2, k = ( k , k , k ) and the steady state ideal is I N = (cid:10) − k x x + k x , k x (cid:11) ⊆ R ( k )[ x ] . Intermediates and steady states
In this subsection we introduce a special type of species of interest: intermediates.
Definition 2.1.
We say that
Y ⊆ S is a subset of intermediates if each Y ∈ Y fulfills: • Y ∈ C and the stoichiometric coefficient of Y in all other complexes is zero, and • there exists at least one reaction having Y as reactant and at least one reactionhaving Y as product.Each Y ∈ Y is called an intermediate .Whenever a set of intermediates Y is given, we partition the set of species into twodisjoint subsets Y = { Y , . . . , Y m } and X = { X , . . . , X n } of non-intermediates. We as-sume further that the set of species is ordered such that the species Y , . . . , Y m are first.With this convention, we let ( y, x ) denote the concentration vector of all species: x is theconcentration vector of the species in X and y of the species in Y . A complex is either anintermediate in Y or it contains only non-intermediates. In the latter case we say that c isa non-intermediate complex .Note that given Y , we refer to the intermediates of the network as the species in Y , eventhough there might be other species in X , regarded as non-intermediates, that fulfill thetwo items in Definition 2.1. Example 2.2.
The most common mechanism involving intermediates is of the followingform: X + E −−→ Y −−→ X (cid:48) + E R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 5 or variations of it by letting one or both reactions being reversible.
Isomerism mechanisms among intermediates are also common: Y −−→ Y (cid:48) Y −− (cid:42)(cid:41) −− Y (cid:48) . Combination of these mechanisms yields to more elaborate networks involving intermedi-ates, as in Examples 2.6 and 3.5 below.
Definition 2.3.
Let Y be a set of intermediates and Y ∈ Y . • A non-intermediate complex c is called an input for Y if there is a directed pathfrom c to Y in the digraph associated with the network, such that all vertices otherthan c belong to Y . • Y is called an (cid:96) -input intermediate if there are (cid:96) inputs for Y . Example 2.4.
Consider the following network with Y = { Y , Y , Y } : X + X −−→ Y −− (cid:42)(cid:41) −− Y −− (cid:42)(cid:41) −− Y −−→ X + X . There are two non-intermediate complexes, X + X and X + X . The species Y , Y , Y are all 1-input intermediates. Note that Y is however the product of two reactions.Consider now the following network with Y = { Y } : X + X −− (cid:42)(cid:41) −− Y −− (cid:42)(cid:41) −− X + X . The species Y is a 2-input intermediate and X + X and X + X are both inputs for Y .2.1. Intermediates and steady states.
Let (cid:101) N be a reaction network with a set of in-termediates Y = { Y , . . . , Y m } . Consider the steady state polynomials of the intermediatesand denote the parameter vector of reaction rate constants by κ (the reason why will bemade clear below). By definition, for every intermediate Y i , the variable y i is only part ofthe monomial y i in (1). Thus, the system with m equations F ( y, x ) = · · · = F m ( y, x ) = 0is linear in y , . . . , y m . It is shown in [11] that this system has a unique solution for fixedpositive values of κ and x , which is further positive. The solution is of the form y i = (cid:88) c ∈C µ i,c x c , where µ i,c ∈ R ≥ ( κ ) , i = 1 , . . . , m. The explicit dependence of µ i,c on κ is omitted from the notation for simplicity. Anexplicit description of µ i,c can be found using the Matrix-Tree theorem on a suitable labeleddigraph, see [11]. Example 2.5.
Consider the following reaction network with X = { X , X , X } and Y = { Y , Y , Y } : X + X κ (cid:47)(cid:111) κ Y κ (cid:47) (cid:47) κ (cid:34) (cid:34) Y κ (cid:47) (cid:47) X Y κ (cid:47)(cid:111) κ κ (cid:59) (cid:59) X κ (cid:79) (cid:79) The linear system in y , y , y that the steady state polynomials of Y , Y , Y define is: κ x x − ( κ + κ + κ ) y = 0 ,κ y − κ y = 0 ,κ y − ( κ + κ ) y + κ x = 0 , and its solution is y = κ κ + κ + κ x x , y = κ κ κ ( κ + κ + κ ) x x ,y = κ κ ( κ + κ )( κ + κ + κ ) x x + κ κ + κ x . R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 6
This gives µ , X X = κ κ + κ + κ , µ , X = 0 , µ , X = 0 ,µ , X X = κ κ κ ( κ + κ + κ ) , µ , X = 0 , µ , X = 0 ,µ , X X = κ κ ( κ + κ )( κ + κ + κ ) , µ , X = κ κ + κ , µ , X = 0 . Example 2.6.
The following digraph corresponds to the Mitogen-Activated Protein Kinasecascade (MAPK) given in [3]: X + E κ −− (cid:42)(cid:41) −− κ Y κ −−→ X + E κ −− (cid:42)(cid:41) −− κ Y κ −−→ X + EX + F κ −− (cid:42)(cid:41) −− κ Y κ −−→ Y κ −− (cid:42)(cid:41) −− κ X + F κ −− (cid:42)(cid:41) −− κ Y κ −−→ Y κ −− (cid:42)(cid:41) −− κ X + F. Species Y , . . . , Y are intermediates. The non-zero coefficients µ i,c are: µ , X E = κ κ + κ , µ , X E = κ κ + κ , µ , X F = κ κ + κ ,µ , X F = κ κ ( κ + κ ) κ , µ , X F = κ κ , µ , X F = κ κ + κ ,µ , X F = κ κ ( κ + κ ) κ , µ , X F = κ κ . Extended and core networks.Definition 2.7.
Let N = ( S , C , R ) and (cid:101) N = ( (cid:101) S , (cid:101) C , (cid:101) R ) be two reaction networks. We saythat (cid:101) N is an extension of N via the addition of intermediates Y , . . . , Y m if(i) Y = { Y , . . . , Y m } is a set of intermediates of (cid:101) N .(ii) S ∪ Y ⊆ (cid:101) S and C ∪ Y ⊆ (cid:101) C .(iii) c → c (cid:48) ∈ R if and only if there is a directed path from c to c (cid:48) in the digraph associatedwith (cid:101) N , such that all vertices other than c and c (cid:48) belong to Y (there might be none).In this case N is called the core network of (cid:101) N . Example 2.8.
The core network associated with the network in Example 2.5 is: X + X k (cid:47) (cid:47) k (cid:38) (cid:38) X X k (cid:59) (cid:59) Example 2.9.
The core network of the network in Example 2.6 has digraph X + E k −−→ X + E k −−→ X + E X + F k −−→ X + F k −−→ X + F. Notations κ, (cid:101) I, (cid:101) F are used to address reaction rate constants, steady state ideal andsteady state polynomials of the extended network respectively. This notation is fixed fromnow on whenever we study extensions via the addition of intermediates.Given (cid:101) N an extension of N via the addition of intermediates Y , . . . , Y m , we define amap φ : R ( k ) −→ R ( κ ) k c → c (cid:48) (cid:55)−→ φ c → c (cid:48) ( κ ) , such that for every reaction c → c (cid:48) ∈ R , φ c → c (cid:48) ( κ ) is the rational function φ c → c (cid:48) ( κ ) = κ c → c (cid:48) + m (cid:88) i =1 κ Y i → c (cid:48) µ i,c , (3)where it is understood that κ c → c (cid:48) = 0, κ Y i → c (cid:48) = 0 if respectively c → c (cid:48) , Y i → c (cid:48) do notbelong to (cid:101) R . Note that φ c → c (cid:48) ( κ ) (cid:54) = 0 for all c → c (cid:48) by Definition 2.7(iii) and that φ c → c (cid:48) ( κ )is a rational function with positive coefficients. R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 7
The map φ extends to a map Φ : R [ k ][ x ] → R ( κ )[ y, x ] . For example, if F i is a steady state polynomial of N , Φ( F i ) is the polynomial obtained byreplacing k c → c (cid:48) by the rational function φ c → c (cid:48) ( κ ). If the rational functions φ c → c (cid:48) ( κ ) arealgebraically independent over R , then φ extends to a map of polynomial ringsΦ : R ( k )[ x ] → R ( κ )[ y, x ] . We explore in Section 4 ways to check whether the algebraic independence condition holds,and provide types of intermediates for which it holds and no extra check is required.We introduce the following polynomials H i ( y, x ) = y i − (cid:88) c ∈C µ i,c x c ∈ R ( κ )[ y, x ] , i = 1 , . . . , m. (4) Theorem 2.10. ([11, Theorems 3.1 and 3.2])
Let (cid:101) N be an extension of N via the additionof intermediates Y , . . . , Y m .(i) The coefficient µ i,c is nonzero if and only if the non-intermediate complex c is aninput for Y i in (cid:101) N .(ii) The set of steady state polynomials of non-intermediate species and the polynomials H , . . . , H m in (4) form a basis of (cid:101) I .(iii) (cid:101) F i (cid:16) (cid:80) c ∈C µ ,c x c , . . . , (cid:80) c ∈C µ m,c x c , x , . . . , x n (cid:17) = Φ( F i ( x )) for i = 1 , . . . , n . Statements (ii) and (iii) of the previous theorem constitute the proof of the followingcorollary.
Corollary 2.11.
Let B be the set of steady state polynomials of N . Then (cid:101) I = (cid:68) Φ( B ) ∪ { H ( y, x ) , . . . , H m ( y, x ) } (cid:69) . We conclude this section with basic properties of Φ.
Lemma 2.12.
With the notation above, assume φ c → c (cid:48) ( κ ) for all c → c (cid:48) ∈ R are alge-braically independent over R . Let B = { f , . . . , f (cid:96) } and B (cid:48) = { f (cid:48) , . . . , f (cid:48) (cid:96) (cid:48) } be two sets in R ( k )[ x ] .(i) If f ∈ (cid:104) B (cid:105) , then Φ( f ) ∈ (cid:104) Φ( B ) (cid:105) .(ii) If (cid:104) B (cid:105) = (cid:104) B (cid:48) (cid:105) , then (cid:104) Φ( B ) (cid:105) = (cid:104) Φ( B (cid:48) ) (cid:105) . Thus Φ( (cid:104) B (cid:105) ) is well defined.Proof. (i) Write f = (cid:80) (cid:96)j =1 α j f j with α j ∈ R ( k )[ x ]. ThenΦ( f ) = (cid:96) (cid:88) j =1 Φ( α j )Φ( f j ) ∈ (cid:104) Φ( B ) (cid:105) . (ii) It is enough to show inclusion ⊆ , since the other inclusion is analogous. If g ∈ (cid:104) Φ( B ) (cid:105) ,we have g = (cid:96) (cid:88) i =1 λ i Φ( f i ) , λ i ∈ R ( κ )[ y, x ] . Since f i ∈ (cid:104) B (cid:48) (cid:105) , we have by (i) that Φ( f i ) ∈ (cid:104) Φ( B (cid:48) ) (cid:105) . In particular, g is an algebraiccombination of the polynomials Φ( f (cid:48) ) , . . . , Φ( f (cid:48) (cid:96) (cid:48) ) with coefficients in R ( κ )[ y, x ]. Thus g ∈(cid:104) Φ( B (cid:48) ) (cid:105) . (cid:3) R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 8 Gr¨obner bases and intermediates
Typically, the values of the reaction rate constants are unknown and reaction networksof interest involve a considerable number of variables. As a consequence, finding a Gr¨obnerbasis of the steady state ideal over the field R ( κ ) can be a demanding task, and some-times even impossible with standard computers. However, the presence of intermediates,a common feature of reaction networks, can reduce the computation time substantially, byexploiting the structure of the steady state polynomials associated with intermediates givenin Theorem 2.10. The main result of this section is Theorem 3.4. Example 3.5 illustrateshow the computation time can be reduced by applying our results.We start with some concepts from computational algebraic geometry.3.1. Monomial orders and Gr¨obner bases.
We follow the notation on Gr¨obner basesfrom [5]. We give here a brief overview of the results required in this text.Given a monomial order on R = K [ x , . . . , x n ], let LM( f ) and LT( f ) denote respectivelythe leading monomial and leading term of f . That is, LT( f ) = α LM( f ) if α is the coefficientof the greatest monomial of f . Then, for a subset A ⊆ R , one defines LT( A ) = (cid:8) LT( f ) | f ∈ A (cid:9) and LM( A ) = (cid:8) LM( f ) | f ∈ A (cid:9) . Clearly, (cid:104)
LT( A ) (cid:105) = (cid:104) LM( A ) (cid:105) . (5)For an ideal I , the initial ideal is the ideal generated by the leading terms of the elementsof I , (cid:104) LT( I ) (cid:105) . A subset G ⊆ I is called a Gr¨obner basis for I if (cid:10) LT( I ) (cid:11) = (cid:10) LT( G ) (cid:11) , (equiv. (cid:10) LM( I ) (cid:11) = (cid:10) LM( G ) (cid:11) ) . A Gr¨obner basis is a basis of I as well. Further, G is a reduced Gr¨obner basis if additionallyfor every element g ∈ G none of its terms can be divided by the leading monomial of anelement in G − { g } , and the coefficient of LM( g ) is 1.Whether a basis of an ideal is a Gr¨obner basis depends on the chosen monomial order.Given an ideal and a monomial order, the Gr¨obner basis is not unique but there is a uniquereduced Gr¨obner basis (see [5]).We will use the following lemma, which follows from Lemma 2.3.1 and Theorem 2.3.2 of[15]. Lemma 3.1.
Let B be a basis of I . If the leading monomials of every pair f, g ∈ B arerelatively prime, then B is a Gr¨obner basis. All monomial orders are defined via a matrix in the following way (though not allmatrices M define a monomial order in this way, [5, 26]). For M ∈ R n × n with full rank,the associated order fulfills x c > x c if the first non-zero entry of the vector M ( c − c ) ispositiveA typical order is the lexicographic monomial order, lex . After choosing a variable order x a > · · · > x a n , lex( x a , . . . , x a n ) is the order defined by the matrix with 1 in positions( i, a i ) for all i = 1 , . . . , n and zero otherwise.Another monomial order of interest is the graded reverse-lexicographic order, abbrevi-ated grevlex . With this order, x c > x c if the total degree of the first monomial is largerthan the second. If they are equal, then the monomial with the smallest variable with leastexponent is the greatest one. Grevlex with order of variables x > · · · > x n is defined bythe matrix . . . . . . −
10 0 . . . − − . . . . The choice of order plays an important role in the computation time for Gr¨obner bases,performing lex typically worse than grevlex. However, lex, as any other elimination type
R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 9 order, has a crucial property on elimination of variables. Given a partitioning of the set ofvariables, { x , . . . , x n } = { x j , . . . , x j n − s }∪{ x i , . . . , x i s } , a monomial order is of eliminationtype if x j (cid:96) , for (cid:96) = 1 , . . . , n − s , is larger than any monomial in K [ x i , . . . , x i s ] [6, § x j , . . . , x j n − s , x i , . . . , x i s ) is of elimination type. If G is a Gr¨obnerbasis of I with respect to an elimination type order as above, then G ∩ K [ x i , . . . , x i s ] is aGr¨obner basis of I ∩ K [ x i , . . . , x i s ] with respect to the induced monomial on K [ x i , . . . , x i s ],which for lex is lex( x i , . . . , x i s ).3.2. Gr¨obner bases and intermediates.
In this subsection we fix a reaction network N and an extension (cid:101) N via the addition of intermediates Y , . . . , Y m . We show that anyGr¨obner basis of the steady state ideal of N can be extended to one of (cid:101) N by simply addingthe polynomials H , . . . , H m given in Equation (4). By default, we order the variables y > · · · > y m > x > · · · > x n . We start with some general lemmas. Lemma 3.2.
Let I = (cid:104) f , f , . . . , f s (cid:105) ⊆ K [ y, x , . . . , x n ] be an ideal such that f i ∈ K [ x , . . . , x n ] for i = 1 , . . . , s and f = y + f (cid:48) , with f (cid:48) ∈ K [ x , . . . , x n ] . Consider a monomial order de-fined by a matrix M whose first row is (cid:0) . . . (cid:1) . Then (cid:10) LT( I ) (cid:11) = (cid:10) y (cid:11) + (cid:10) LT( (cid:104) f , . . . , f s (cid:105) ) (cid:11) . Further given G ⊆ K [ x , . . . , x n ] , G is a Gr¨obner basis of (cid:104) f , . . . , f s (cid:105) if and only if { f } ∪ G is a Gr¨obner basis of I .Proof. By the choice of monomial order, the monomial y is larger than any monomial notinvolving y . Consider a reduced Gr¨obner basis G (cid:48) of (cid:104) f , . . . , f s (cid:105) . Then the leading termsof the elements in G (cid:48) are relatively prime with each other and with the leading term of f .Since { f } ∪ G (cid:48) is a basis of I , then by Lemma 3.1 { f } ∪ G (cid:48) is a Gr¨obner basis of I . Now,the initial ideal of I is generated by the leading terms of { f } ∪ G (cid:48) . So: (cid:104) LT( I ) (cid:105) = (cid:104) LT( { f } ∪ G (cid:48) ) (cid:105) = (cid:104){ y } ∪ LT( G (cid:48) ) (cid:105) = (cid:104) y (cid:105) + (cid:104) LT( G (cid:48) ) (cid:105) = (cid:104) y (cid:105) + (cid:104) LT( (cid:104) f , . . . , f s (cid:105) ) (cid:105) . This proves the first part of the lemma.For the second part, note that (cid:104) y (cid:105) + (cid:104) LT( G ) (cid:105) = (cid:104){ LT( f ) } ∪ LT( G ) (cid:105) = (cid:104) LT( { f } ∪ G ) (cid:105) . Using this equality and the first part of the lemma, we have { f } ∪ G is a Gr¨obner basisof I if and only if (cid:104) y (cid:105) + (cid:104) LT( G ) (cid:105) = (cid:104) y (cid:105) + (cid:104) LT( (cid:104) f , . . . , f s (cid:105) ) (cid:105) . Since y is not part of anypolynomial in G , this equality holds if and only if (cid:104) LT( G ) (cid:105) = (cid:104) LT( (cid:104) f , . . . , f s (cid:105) ) (cid:105) , i.e. G isa Gr¨obner basis of (cid:104) f , . . . , f s (cid:105) . (cid:3) Recall that we write I ⊆ R ( k )[ x , . . . , x n ] and (cid:101) I ⊆ R ( κ )[ y , . . . , y m , x , . . . , x n ] for thesteady state ideals of N and (cid:101) N respectively. For the rest of the section, we assume thatthe rational functions φ c → c (cid:48) ( κ ) are algebraically independent over R , such that Φ( A )is defined for all subsets A of R ( k )[ x ].For an arbitrary basis B of I , define (cid:101) B = Φ( B ) ∪ (cid:8) H ( y, x ) , . . . , H m ( y, x ) (cid:9) ⊆ R ( κ )[ y, x ] . (6) Lemma 3.3. If B is a basis of I ⊆ R ( k )[ x ] , then (cid:101) B is a basis of (cid:101) I ⊆ R ( κ )[ y, x ] .Proof. Let B (cid:48) be the set of steady state polynomials of N . By Corollary 2.11 (cid:101) I = (cid:68) Φ( B (cid:48) ) ∪ { H ( y, x ) , . . . , H m ( y, x ) } (cid:69) . Let now B be an arbitrary basis of I . Then (cid:104) B (cid:105) = I = (cid:104) B (cid:48) (cid:105) and thus by Lemma 2.12(ii), (cid:104) Φ( B ) (cid:105) = (cid:104) Φ( B (cid:48) ) (cid:105) . Therefore (cid:68) Φ( B ) ∪ { H ( y, x ) , . . . , H m ( y, x ) } (cid:69) = (cid:68) Φ( B (cid:48) ) ∪ { H ( y, x ) , . . . , H m ( y, x ) } (cid:69) = (cid:101) I. R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 10
This completes the proof. (cid:3)
Let rem( p, q ) be the remainder of the division of the polynomial p by q . Theorem 3.4.
Fix a monomial order on R ( k )[ x ] associated with an n × n matrix Q , andlet G be a Gr¨obner basis of I with this order. Then, (cid:101) G is a Gr¨obner basis of (cid:101) I with themonomial order on R ( κ )[ y, x ] associated with the matrix (cid:101) Q = (cid:18) Id m Q (cid:19) , (7) where Id m is the identity matrix of size m .If G is reduced, then Φ( G ) ∪ (cid:110) y i − rem (cid:0) (cid:80) c ∈C µ i,c x c , Φ( G ) (cid:1)(cid:111) is the reduced Gr¨obner basisof (cid:101) I .Proof. First note that by the monomial order given by (cid:101) Q , we have y > · · · > y m > x i forall i = 1 , . . . , n . Also, the polynomial H i has degree one in y i and none of the elements ofΦ( G ) ∪ { H j | j (cid:54) = i } involves y i .Let us assume we have shown that Φ( G ) is a Gr¨obner basis of (cid:104) Φ( G ) (cid:105) with the givenorder, that is (cid:10) LT( (cid:104) Φ( G ) (cid:105) ) (cid:11) = (cid:10) LT(Φ( G )) (cid:11) . (8)Then by Lemmas 3.2 and 3.3, Φ( G ) ∪ { H ( y, x ) , . . . , H m ( y, x ) } is a Gr¨obner basis of (cid:101) I .Therefore the first part of the statement holds provided (8) holds.Let us show (8). We start by noting that for a subset J in R ( k )[ x ], the set LM( J )consists only of monomials in x , . . . , x n , and thus is naturally included in R ( κ )[ y, x ] aswell. Further LM( J ) = LM(Φ( J )) . (9)Let G (cid:48) be a reduced Gr¨obner basis of I . Since G (cid:48) is reduced, pairs of monomials inLM( G (cid:48) ) = LM(Φ( G (cid:48) )) are relatively prime. Since Φ( G (cid:48) ) is a basis of (cid:104) Φ( G (cid:48) ) (cid:105) , then byLemma 3.1 and Equation (5), it is actually a Gr¨obner basis and (8) holds for G (cid:48) . Now,consider an arbitrary Gr¨obner basis G of I . In R ( k )[ x ] it holds (cid:104) LM( G ) (cid:105) = (cid:104) LM( G (cid:48) ) (cid:105) . (10)This means that every monomial in (cid:104) LM( G (cid:48) ) (cid:105) is divisible by a monomial in (cid:104) LM( G ) (cid:105) andviceversa [5, § R ( κ )[ y, x ], (10) holds also in R ( κ )[ y, x ]. Combined with (9) this gives (cid:10) LM(Φ( G )) (cid:11) = (cid:10) LM(Φ( G (cid:48) )) (cid:11) . By Lemma 2.12(ii), (cid:104) G (cid:105) = (cid:104) G (cid:48) (cid:105) in R ( k )[ x ] implies (cid:104) Φ( G ) (cid:105) = (cid:104) Φ( G (cid:48) ) (cid:105) . Thus in R ( κ )[ y, x ]we have (cid:10) LM(Φ( G )) (cid:11) = (cid:10) LM(Φ( G (cid:48) )) (cid:11) = (cid:10) LM( (cid:104) Φ( G (cid:48) ) (cid:105) ) (cid:11) = (cid:10) LM( (cid:104) Φ( G ) (cid:105) ) (cid:11) . This shows that (8) holds.The second part of the lemma is clear from the definition of a reduced Gr¨obner basisand using that Φ( G ) ∪ { y i − rem (cid:0) (cid:80) c ∈C µ i,c x c , Φ( G ) (cid:1) } is also a Gr¨obner basis. (cid:3) From the computational point of view, Theorem 3.4 is very useful. Instead of comput-ing a Gr¨obner basis of (cid:101) I directly, one can first compute a Gr¨obner basis G for the corenetwork, with a smaller number of variables and polynomials, then add the polynomials y i − (cid:80) c ∈C µ i,c x c , and, finally, simplify them using polynomial division by Φ( G ). The sec-ond step involves only linear algebra. A possible issue here is to verify that the rationalfunctions φ c → c (cid:48) are algebraically independent. We provide in Section 4 a list of networkstructures involving intermediates for which the condition is fulfilled. R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 11 X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X X + X κ −− (cid:42)(cid:41) −− κ X κ −−→ X + X X κ −−→ X Figure 1.
Reaction network of Example 3.5.
Example 3.5.
An interesting example to show the advantage of using Theorem 3.4 isExample 4.4 of [4]. We consider the reaction network (cid:101) N with associated digraph given inFigure 1.This reaction network has 29 species and 46 reactions. Therefore the steady state idealis generated by 29 polynomials in 29 variables and 46 parameters. Using Singular [7] andmonomial order grevlex with x > · · · > x (the same monomial order that is used in [4]),it took between 110 and 115 seconds to compute the reduced Gr¨obner basis. This basisconsists of 169 polynomials.Now we consider the monomial order introduced in Theorem 3.4 for the removal of the15 intermediates: X , X , X , X , X , X , X , X , X , X , X , X , X , X , X . The original network (cid:101) N is an extension of the following core network N with 14 speciesand 16 reactions: X + X k −−→ X + X X + X k −−→ X + X k −−→ X + X X + X k −−→ X + X X + X k −−→ X + X k −−→ X + X X + X k −−→ X + X X + X k −−→ X + X k −−→ X + X X + X k −−→ X + X X + X k −−→ X + X k −−→ X + X X + X k −−→ X + X X + X k −−→ X + X k −−→ X + X X k −−→ X . The functions φ c → c (cid:48) are algebraically independent over R by Corollary 4.6 in Section 4.We consider grevlex with x > · · · > x for the monomials corresponding to N . The Information about the processor: Intel(R) Core(TM) i5-3570 CPU @3.4GHz 3.4GHz with 8GB RAM.We report the interval of obtained times after several runs of Singular, computed in milliseconds.
R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 12 monomial order in Theorem 3.4 is then associated with the following matrix (cid:101) Q = . . . . . . −
10 0 . . . − − . . . , and order of variables x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x > x . The reduced Gr¨obner basis of (cid:101) I with this monomial order has 33 polynomials and it takesabout 96 seconds to compute it directly with Singular. Alternatively the strategy outlinedin Theorem 3.4 can be applied. The steady state ideal I of N is generated by 11 polynomialsin 14 variables and 16 parameters. Using Singular, the reduced Gr¨obner basis of I has 18polynomials and its computation takes less than a millisecond. The computation time forthe polynomials H i ( y, x ) is neglectable, since they are found by solving 15 independentlinear equations. Therefore the reduced Gr¨obner basis of the ideal of the original systemhas 18+15=33 polynomials and can be computed in less than a millisecond.We conclude that in general, regarding computational time, the monomial order in-troduced in Theorem 3.4 is a good choice for networks with intermediates, and further,by applying the strategy of Theorem 3.4 we reduce the computation time considerably,compared with direct computation of the reduced Gr¨obner basis. Remark 3.6.
Theorem 3.4 holds regardless the choice of method to compute a Gr¨obnerbasis. Since the computation of the polynomials H i is simple linear algebra, even forthe fastest available methods for the computation of Gr¨obner bases, decomposing thecomputation as in Theorem 3.4 should be faster than direct computation of the basis ofthe steady state ideal of (cid:101) N . Remark 3.7.
For polynomials with integer coefficients, it is usually faster to compute aGr¨obner basis using the so-called p-modular approach, see e.g. [23, 29]. These methodsfirst choose a so-called lucky prime and compute a Gr¨obner basis of the ideal in Z p [ x ]. Thenthe coefficients of this Gr¨obner basis are lifted to a Gr¨obner basis in Q [ x ]. For the sake ofcomparison, we also computed how long it takes to find a Gr¨obner basis using p-modularapproaches on the extended network in Example 3.5 with grevlex and x > · · · > x .Using the largest prime number in Singular, p = 32003, it takes 127 seconds to computethe Gr¨obner basis over Z . Since coefficients in the starting basis are 1 or −
1, one maythink that p = 2 is a lucky prime. It took 97 seconds to compute the Gr¨obner basis over Z . These times are larger than the times reported in Example 3.5 (and these Gr¨obnerbases still need to be lifted to Q ( κ )[ y, x ]).An important consequence of Theorem 3.4 concerns parameter-free model discrimina-tion. In this setting one seeks elements of the steady state ideal (cid:101) I involving only theconcentrations of species that are experimentally measurable. These elements are called invariants . Each invariant implies that there is a set of monomials that lie on a hyperplane,and the hypothesis of coplanary is then tested using experimental data [13, 14, 18–20]. Thisapproach is attractive because it does not require knowing the values of the reaction rateconstants. R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 13
Experimentally measurable species do not typically involve intermediates. In this case,Theorem 3.4 tells us that invariants on the non-intermediate species can be computeddirectly from the core network, using elimination ideals.
Corollary 3.8.
Let N be a reaction network and (cid:101) N an extension of it via the addition of m intermediates Y , . . . , Y m . Let X i , . . . , X i p be non-intermediates. Then (cid:101) I ∩ R ( κ )[ x i , . . . , x i p ] = Φ( I ∩ R ( k )[ x i , . . . , x i p ]) . Proof.
For simplicity, assume { i , . . . , i p } = { n − p + 1 , . . . , n } and let x = ( x n − p +1 , . . . , x n ).Consider the monomial order lex( y , . . . , y m , x , . . . , x n ) on R ( κ )[ y, x ], and lex( x , . . . , x n )on R ( k )[ x ]. Let G be a Gr¨obner basis of I . By Theorem 3.4, (cid:101) G is a Gr¨obner basis of (cid:101) I .By the properties of lex and Lemma 2.12(ii) we have (cid:101) I ∩ R ( κ )[ x ] = (cid:104) (cid:101) G ∩ R ( κ )[ x ] (cid:105) = (cid:104) Φ( G ∩ R ( k )[ x ]) (cid:105) = Φ( I ∩ R ( k )[ x ]) . This concludes the proof. (cid:3)
Note that the monomial order on R ( κ )[ y, x ] given in Theorem 3.4 is of elimination typewith respect to the partition { y , . . . , y m } ∪ { x , . . . , x n } . Example 3.9.
Consider the network in Example 2.6 and its core network in Example2.9. In order to find invariants of the extended network involving the concentration of thenon-intermediate species
E, X , X , X , we consider the ideal I ∩ R ( k )[ e, x , x , x ], whichis generated by the polynomial e ( k k x x − k k x ) . We have φ ( k , k , k , k ) = (cid:16) κ κ κ + κ , κ κ κ + κ , κ κ κ + κ , κ κ κ + κ (cid:17) . The functions φ c → c (cid:48) are algebraically independent over R by Corollary 4.6. By Corollary3.8 the ideal (cid:101) I ∩ R ( κ )[ e, x , x , x ] is generated by the polynomial e (cid:16) κ κ κ + κ κ κ κ + κ x x − κ κ κ + κ κ κ κ + κ x (cid:17) . Detecting binomial steady state ideals. A binomial is a polynomial having atmost two terms. An ideal is said to be binomial if it admits a set of generators consistingof binomials only. By [8, Corollary 1.2], an ideal is binomial if and only if any reducedGr¨obner basis (with respect to any monomial order) consists of binomials.It is of biological relevance in the study of reaction networks to determine whether thereexists a choice of reaction rate constants k for which there are multiple positive steady statesin some coset x + S defined by the vector subspace S that contains the image of F k (seeSection 1). This property is termed multistationarity . If the steady state ideal is binomial,then there exist efficient ways to determine whether the network admits multistationarity[22, 24, 25]. This leads to the problem of determining whether an ideal is binomial, and incase it is, of finding a binomial basis of it. As noted, both questions can be addressed byfinding a Gr¨obner basis of the steady state ideal of the network. Thus, for networks withintermediates, our results can be applied also to detect binomial steady state ideals.Recall that we are assuming that the rational functions φ c → c (cid:48) ( κ ) are algebraicallyindependent over R . Theorem 3.10.
Let N be a reaction network and (cid:101) N an extension of it via the addition of m intermediates Y , . . . , Y m .The steady state (cid:101) I is binomial if and only if • I is binomial, and, • for any reduced Gr¨obner basis G of I and for every i = 1 , . . . , m , the remainder ofthe division of (cid:80) c ∈C µ i,c x c by Φ( G ) has at most one term. R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 14
Proof.
Fix any monomial order on R ( k )[ x , . . . , x n ] associated with an n × n matrix Q andconsider the monomial order with matrix (cid:101) Q from Theorem 3.4. Let G be the reducedGr¨obner basis of I and (cid:101) G (cid:48) = Φ( G ) ∪ (cid:110) y i − rem (cid:0) (cid:88) c ∈C µ i,c x c , Φ( G ) (cid:1)(cid:111) the reduced Gr¨obner basis of (cid:101) I (cf. Theorem 3.4). Using that an ideal is binomial if andonly if any reduced Gr¨obner basis consists of binomials, the theorem is a consequence ofthe following two facts: • By definition, (cid:101) G (cid:48) consists of binomials if and only if Φ( G ) is a set of binomials andthe remainder of the division of (cid:80) c ∈C µ i,c x c by Φ( G ) has at most one term. • By the algebraic independence of φ c → c (cid:48) , Φ( G ) consists of binomials if and only if G does. (cid:3) Since the polynomial (cid:80) c ∈C µ i,c x c has exactly one term for 1-input intermediates, wereadily obtain the following corollary. Corollary 3.11.
Let N be a reaction network and (cid:101) N an extension of it via the additionof m intermediates Y , . . . , Y m . Then (cid:101) I is binomial if and only if I is binomial. Since 1-input intermediates are the most abundant form of intermediates found in real-istic networks, this corollary implies that in order to check whether a steady state ideal isbinomial, we can often remove intermediates and check whether the steady state ideal ofthe core network is binomial.
Example 3.12.
Consider the network in Example 2.5 and its core network given in Exam-ple 2.8. The functions φ c → c (cid:48) are algebraically independent over R by Example 4.1. Sincethe steady state ideal of N is (cid:104)− ( k − k ) x x − k x , ( k − k ) x x + 2 k x (cid:105) , the core network has a binomial steady state ideal. The reduced Gr¨obner basis for thisideal with monomial order lex( x , x , x ) is G = (cid:110) x − ( k − k )2 k x x (cid:111) . We apply Theorem 3.10 to conclude that the steady state ideal of the extended networkis also binomial. The intermediates Y , Y are 1-input intermediates and hence the re-mainder condition of the theorem is automatically fulfilled. For the intermediate Y ,rem (cid:0) µ , X X x x + µ , X X x x , Φ( G ) (cid:1) has a single term with monomial x x . Thereforewe conclude that the extended network also has a binomial steady state ideal.The following example shows that extended networks with multi-input intermediatesmight not have binomial steady state ideals, even though their core networks have. Example 3.13.
Consider the network given in Example 2.6 and its core network givenin Example 2.9. The steady state ideal of the core network is binomial with basis B = { k x e − k x f, k x e − k x f } . The intermediates Y and Y are 2-input intermediates.The remainder of the division of µ , X E x f + µ , X F x f by Φ( G ) for G the reducedGr¨obner basis of I with the monomial order lex( x , x , x , f, e ) is κ κ x f + κ κ κ κ + κ κ x f, which has two terms. Therefore by Theorem 3.10 the steady state ideal of the network inExample 2.6 is not binomial. R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 15
Remark 3.14.
In [4], a method for determining whether a homogeneous ideal is binomial isintroduced. The method avoids the computation of Gr¨obner bases and is regarded as a fastmethod. If the steady state ideal of the core network is homogeneous, then Theorem 3.10 orCorollary 3.11 in combination with this method provide a fast procedure to detect binomialsteady state ideals.Interestingly, steady state polynomials of core networks are often homogeneous of degreetwo, since it is common that non-intermediate species appear in complexes of the form X i + X j , yielding quadratic terms in the steady state polynomials. This is for example thecase for so-called Post-Translational Modification Networks [28].4. Algebraic independence
In this section we discuss how to check whether the functions φ c → c (cid:48) are algebraicallyindependent over R and provide classes of intermediates for which this property holds.Consider a set of rational functions A = (cid:8) f g , . . . , f m g m (cid:9) ⊆ R ( x , . . . , x n ). By § III.7, TheoremIII, in [16], the set A is algebraically independent over R if and only if the rank of theassociated Jacobian matrix (cid:16) ∂ ( f i /g i ) ∂x j (cid:17) i,j over R ( x ) is m .Another way to check algebraic independence that requires the computation of a Gr¨obnerbasis is as follows. Let ϕ be the function on R n minus the zero locus of the product g · · · g m defined by x = ( x , . . . , x n ) (cid:55)→ (cid:18) f ( x ) g ( x ) , . . . , f m ( x ) g m ( x ) (cid:19) . By § ϕ ) is the variety associated with the ideal J := (cid:10) g T − f , . . . , g m T m − f m , − yg · · · g m (cid:11) ∩ R [ T , . . . , T m ] . Since the sets of polynomials vanishing on a set and on its closure agree (see [6] afterDefinition 2 in § A is algebraically independent over R if and only if J = { } . Example 4.1.
The functions φ c → c (cid:48) of Examples 2.5 and 2.8 are φ X X → X ( κ ) = κ µ , X X + κ µ , X X = κ κ κ + κ + κ + κ κ κ ( κ + κ )( κ + κ + κ ) ,φ X X → X ( κ ) = κ µ , X X = κ κ κ ( κ + κ )( κ + κ + κ ) ,φ X → X ( κ ) = κ + κ µ , X = κ + κ κ κ + κ . We find that J = { } . Hence the algebraic independence condition holds for the networkin Example 2.8. Alternatively, one easily checks that the associated Jacobian matrix hasrank 3.The computations above can be simplified by taking into account what parameters occurin each of the rational functions. Definition 4.2.
Let (cid:101) N be an extension of N via the addition of the intermediates { Y , . . . , Y m } .Consider the digraph associated with (cid:101) N . Let Y , . . . , Y t (cid:48) denote the vertex sets of the con-nected components of the subgraph induced by the subset of vertices { Y , . . . , Y m } .Let R (cid:48) ⊆ R be the subset of reactions of the core network that are not in (cid:101) R . Thesereactions arise necessarily from paths through intermediates. We say that two reactions r : c → c (cid:48) , r : c → c (cid:48) ∈ R (cid:48) overlap if there exist paths through intermediates c → Y i → · · · → Y i p → c (cid:48) , c → Y j → · · · → Y j q → c (cid:48) with all intermediates belonging to the same set Y i .Consider the equivalence relation on R (cid:48) generated by the overlap relation: r ∼ r (cid:48) if andonly if there exist r = r, r , . . . , r p = r (cid:48) such that r i , r i +1 overlap for all i = 0 , . . . , p − R (cid:48) , . . . , R (cid:48) t be the equivalence classes of this equivalence relation. R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 16
Example 4.3.
Consider the network in Example 2.8. The set R (cid:48) consists of two reactions X + X → X and X + X → X . The subgraph of the digraph associated with (cid:101) N induced by the set of intermediates is connected. Thus the two reactions of R (cid:48) areequivalent and there is one equivalence class. Lemma 4.4.
The set { φ c → c (cid:48) ( κ ) | c → c (cid:48) ∈ R} is algebraically independent over R ifand only if the set { φ c → c (cid:48) ( κ ) | c → c (cid:48) ∈ R (cid:48) i } is algebraically independent over R for all i = 1 , . . . , t .Proof. Since R (cid:48) i ⊆ R for all i = 1 , . . . , t , the forward implication is clear.To prove the reverse implication, assume that the sets T i = { φ c → c (cid:48) ( κ ) | c → c (cid:48) ∈ R (cid:48) i } are algebraically independent over R for all i = 1 , . . . , t . By construction, the sets ofparameters appearing in the rational functions φ c → c (cid:48) ( κ ) are disjoint for two reactions indifferent equivalence classes. Therefore the union of the sets T , . . . , T t is algebraicallyindependent over R . Furthermore if c → c (cid:48) ∈ R \ R (cid:48) , then the parameter κ c → c (cid:48) appearsonly in φ c → c (cid:48) ( κ ). As a consequence the set t (cid:91) i =1 T i ∪ { φ c → c (cid:48) ( κ ) | c → c (cid:48) ∈ R \ R (cid:48) } = { φ c → c (cid:48) ( κ ) | c → c (cid:48) ∈ R} is algebraically independent over R . (cid:3) Example 4.5.
Consider the network in Example 4.3. The algebraic independence of thefunctions φ c → c (cid:48) ( κ ) for all reactions c → c (cid:48) in R follows in this case from the algebraicindependence of the functions φ c → c (cid:48) ( κ ) for the reactions X + X → X and X + X → X . Corollary 4.6. If R (cid:48) = ∅ or each of the equivalence classes R (cid:48) , . . . , R (cid:48) t consist of onereaction, then the rational functions φ c → c (cid:48) ( κ ) are algebraically independent over R . For the networks in Example 2.6 and Example 3.5, each of the equivalence classes consistof one reaction. Therefore, by Corollary 4.6, the algebraic independence condition holds.We next show that the algebraic independence condition holds for specific classes ofintermediates without the need of doing any extra computation.
Lemma 4.7.
For the following extension networks, with intermediates Y , . . . , Y m , the set { φ c → c (cid:48) ( κ ) | c → c (cid:48) ∈ R} is algebraically independent over R .(i) c ←−→ Y ←−→ Y ←−→ . . . ←−→ Y m ←−→ c (cid:48) , provided { Y , . . . , Y m } is a set ofintermediates and where ←−→ means the reaction can be irreversible or reversible.(ii) c c c Y Y . . . Y m(cid:96) (cid:68) (cid:68) (cid:96) (cid:60) (cid:60) (cid:96) p (cid:34) (cid:34) ... c p with an arbitrary digraph structure among the complexes c , Y , . . . , Y m such that thereexists a directed path from c to Y m . R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 17 (iii) c κ (cid:47)(cid:111) κ Y κ (cid:47)(cid:111) κ (cid:96) , (cid:22) (cid:22) (cid:96) ,t (cid:27) (cid:27) Y κ (cid:47)(cid:111) κ (cid:96) , (cid:22) (cid:22) (cid:96) ,t (cid:27) (cid:27) . . . κ m − (cid:47)(cid:111) κ m Y m(cid:96) m, (cid:22) (cid:22) (cid:96) m,tm (cid:27) (cid:27) c , c , c m, ... ... ... c ,t c ,t c m,t m where some of the reactions with label κ i might not exist, and for each ≤ i ≤ m ,either t i ≥ .Proof. We start by recalling how to find µ i,c using a labeled digraph (see proof of Theorem2 of the electronic supplementary material of [11]). For each non-intermediate complex c ,consider the labeled digraph (cid:98) G c with vertex set { Y , . . . , Y m , (cid:63) } and labeled edges Y i κ Yi → Yj −−−−→ Y j if Y i → Y j ∈ (cid:101) R , (cid:63) κ c → Yi x c −−−−−→ Y i if c → Y i ∈ (cid:101) R and Y i β i −→ (cid:63) with β i = (cid:80) Y i → c (cid:48) κ Y i → c (cid:48) if β i (cid:54) = 0.For every vertex v of (cid:98) G c define θ ( v ) as the set of all spanning trees rooted at v . Givensuch a tree τ , let π ( τ ) be the product of the labels of the edges of τ . Then µ i,c = (cid:80) τ ∈ θ ( Y i ) π ( τ ) (cid:80) τ ∈ θ ( (cid:63) ) π ( τ ) . (11)(i) If one of the reactions is irreversible, then the core network consists of exactly onereaction, either c → c (cid:48) or c (cid:48) → c , and the set { φ c → c (cid:48) ( κ ) | c → c (cid:48) ∈ R} is algebraicallyindependent over R .If all reactions are reversible, we write c κ −− (cid:42)(cid:41) −− κ Y κ −− (cid:42)(cid:41) −− κ Y κ −− (cid:42)(cid:41) −− κ . . . κ m − −−−− (cid:42)(cid:41) −−−− κ m Y m κ m +1 −−−− (cid:42)(cid:41) −−−− κ m +2 c (cid:48) , and we have φ c (cid:48) → c ( κ ) = κ µ ,c (cid:48) , φ c → c (cid:48) ( κ ) = κ m +1 µ m,c . By the expressions for µ ,c (cid:48) and µ m,c in (11), both rational functions have the same denominator and κ m +1 is not part oftheir numerator. Therefore, algebraic independence of κ µ ,c (cid:48) and κ m +1 µ m,c follows fromthe algebraic independence of the numerators of these two rational functions. Since κ m +1 is a factor of φ c → c (cid:48) ( κ ) and is not part of the numerator of φ c (cid:48) → c ( κ ), the two functions φ c → c (cid:48) ( κ ) , φ c (cid:48) → c ( κ ) are algebraically independent over R .(ii) We have φ c → c i ( κ ) = (cid:96) i µ m,c for i = 1 , . . . , p . Thus the set { φ c → c i | ≤ i ≤ p } isalgebraically independent over R if and only if { (cid:96) i | ≤ i ≤ p } is, which clearly holds.(iii) The reactions of the core network are of the form c → c i,j . We consider the graph (cid:98) G c (removing the edges for which there is no reaction): (cid:63) κ x c (cid:48) (cid:48) Y κ (cid:47) κ + (cid:80) t j =1 (cid:96) ,j (cid:114) (cid:114) Y κ (cid:111) κ (cid:47) (cid:80) t j =1 (cid:96) ,j (cid:116) (cid:116) . . . κ (cid:111) κ m − (cid:47) Y mκ m (cid:111) (cid:80) tmj =1 (cid:96) m,j (cid:118) (cid:118) We have φ c → c i,j ( κ ) = (cid:96) i,j µ i,c . The denominators of the rational functions µ i,c as given in(11) agree. Therefore it is enough to check that the polynomials ρ i,j := (cid:96) i,j (cid:80) τ ∈ θ ( Y i ) π ( τ )for all i, j are algebraically independent over R . a spanning tree is rooted at v if v is the only vertex with no outgoing edges R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 18
For each 1 ≤ i ≤ m , there exists a spanning tree rooted at Y i involving an edge of theform Y j → (cid:63) only for j (cid:9) i . Now consider the smallest index i such that there exists acomplex c i,j . The parameter (cid:96) i,j appears in a polynomial ρ i ,i only for i = i . Hencethe polynomials ρ i ,i are algebraically independent if and only if they are for i > i . Weproceed in the same way now considering the smallest index k > i such that there exists acomplex c k,j . This process terminates in at most m steps. (cid:3) Corollary 4.6 and Lemma 4.7(i) show that typical rational functions arising from realisticnetworks, such as those built from the mechanism in Example 2.2, fulfil the algebraicindependence condition.5.
Another class of species: enzymes
In this final section we consider another class of species for which reduction mechanismshave also been defined, namely enzymes, and study how Gr¨obner bases of extended andreduced networks relate.5.1.
Enzymes.
A species E ∈ S is an enzyme if for every reaction the stoichiometriccoefficient of E in the reactant and the product agree [21]. This automatically gives thatthe steady state polynomial of E is identically zero, and implies that the concentration of E is constant in time and only depends on the initial amount e of E . For example, E and F are enzymes in the network of Example 2.9.The core network obtained by removal of E consists of simply removing E from eachside of the reaction (this is an example of an embedded network, see [17]). For example, areaction X + E κ −→ X + E becomes X k −→ X . (12)After fixing the initial amount of enzyme e , the steady states of the extended networksatisfying that the concentration of E is e agree with the steady states of the core networkwith k = e κ .This might lead one to think that enzymes are redundant and that similar properties asthose that hold for intermediates also hold for enzymes. For example, one might think thereis an easy way to obtain a Gr¨obner basis of the steady state ideal of the extended networkfrom one of the core network, or that a binomial steady state ideal remains binomial uponremoval of intermediates. But this is not the case, as the following examples illustrate. Example 5.1.
Let N be the network2 X k (cid:47) (cid:47) k (cid:50) (cid:50) X k (cid:47) (cid:47) X A binomial basis of the steady state ideal is {− k x + ( k − k ) x } . Now consdier thefollowing network by adding one enzyme E :2 X κ −−→ X κ −−→ X X + E κ −−→ X + E. A reduced Gr¨obner basis of its steady state ideal is { x − κ κ x + κ κ x e } , and hence thisideal is not binomial.The previous example suggests the following: Consider a reaction as in (12). One mightobtain a Gr¨obner basis of the steady state ideal of the extended network by considering aGr¨obner basis of the steady state ideal of the core network and substituting the parameter κ by k e . The following example gives a negative answer to this question. R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 19
Example 5.2.
Let N be the following network X k (cid:40) (cid:40) X k (cid:47) (cid:47) X . X k (cid:54) (cid:54) The set of steady state polynomials is {− k x − k x − k x , k x } . With every arbitrary monomial order on R ( k )[ x ], the reduced Gr¨obner basis of the steadystate ideal is { x } .Let now N (cid:48) be the extension of N via the enzyme E : X k (cid:40) (cid:40) X + E k (cid:47) (cid:47) X + E. X k (cid:54) (cid:54) The set of steady state polynomials of N (cid:48) is {− κ x − κ x − κ x e, κ x e } . The steady state ideal is different from (cid:104) x (cid:105) . Thus, there is not a monomial order on R ( κ )[ x, e ] for which the reduced Gr¨obner basis can be obtained from the set { x } bymaking the substitution k = κ e . Example 5.3.
When a binomial basis of the steady state ideal is obtained from linearcombinations of the steady state polynomials (see [4]), then the steady state ideal of thecore network is binomial if and only if that of the extended network is.
Acknowledgements.
This work has been supported by the Danish Research Council forIndependent Research. We thank Martin Helmer, Ang´elica Torres and Carsten Wiuf forcomments on previous versions of this manuscript.
References [1] M. Bardet, J. Faug`ere, and B. Salvy. On the complexity of the F Gr¨obner basisalgorithm.
J. Symb. Comput. , 70:49–70, 2015.[2] D. Bayer and M. Stillman. A theorem on refining division orders by the reverselexicographc order.
Duke Math. J. , 55(2):321–328, 1987.[3] R. Bradford, J. H. Davenport, M. England, H. Errami, V. Gerdt, D. Grigoriev,C. Hoyt, M. Koˇsta, O. Radulescu, T. Sturm, and A. Weber. A case study on theparametric occurrence of multiple steady states. In
Proceedings of the InternationalSymposium on Symbolic and Algebraic Computation, ISSAC , pages 45–52. Associationfor Computing Machinery, 2017.[4] C. Conradi and T. Kahle. Detecting binomiality.
Adv. Appl. Math. , 71(C):52–67, 2015.[5] D. Cox, J. Little, and D. O’Shea.
Using Algebraic Geometry . Springer-Verlag NewYork, 2nd edition, 2005.[6] D. Cox, J. Little, and D. O’Shea.
Ideals, Varieties, and Algorithms . Springer Interna-tional Publishing, 4th edition, 2015.[7] W. Decker, G. Greuel, G. Pfister, and H. Sch¨onemann.
Singular ,2016.[8] D. Eisenbud and B. Sturmfels. Binomial ideals.
Duke Math. J. , 84(1):1–45, 1996.[9] J. Faug`ere. A new efficient algorithm for computing Gr¨obner bases ( F ). J. Pure Appl.Algebr. , 139:61–88, 1999.[10] M. Feinberg. Lectures on chemical reaction networks. Available online at , 1980.
R ¨OBNER BASES OF REACTION NETWORKS WITH INTERMEDIATE SPECIES 20 [11] E. Feliu and C. Wiuf. Simplifying biochemical models with intermediate species.
J.R. Soc. Interface , 10:20130484, 2013.[12] J. Gunawardena. Chemical reaction network theory for in-silico biologists. Availableonline at http://vcp.med.harvard.edu/papers/crnt.pdf , 2003.[13] J. Gunawardena. Distributivity and processivity in multisite phosphorylation can bedistinguished through steady-state invariants.
Biophys. J. , 93(11):3828–3834, 2007.[14] H. A. Harrington, K. L. Ho, T. Thorne, and M. P. H. Stumpf. Parameter-free modeldiscrimination criterion based on steady-state coplanarity.
Proc. Natl. Acad. Sci. U.S. A. , 109(39):15746–15751, 2012.[15] J. Herzog and T. Hibi.
Monomial Ideals . Springer-Verlag London, 1st edition, 2011.[16] W. V. D. Hodge.
Methods of Algebraic Geometry , volume 1. Cambridge UniversityPress, 1st edition, 1953.[17] B. Joshi and Shiu J. Atoms of multistationarity in chemical reaction networks.
J.Math. Chem. , 51(1):153–178, 2013.[18] R. L. Karp, M. P´erez Mill´an, T. Dasgupta, A. Dickenstein, and J. Gunawardena.Complex-linear invariants of biochemical networks.
J. Theor. Biol. , 311(21):130–138,2012.[19] A. L. MacLean, Z. Rosen, H. M. Byrne, and H. A. Harrington. Parameter-free methodsdistinguish wnt pathway models and guide design of experiments.
Proc. Natl. Acad.Sci. U. S. A. , 112(9):2652–2657, 2015.[20] A. K. Manrai and J. Gunawardena. The geometry of multisite phosphorylation.
Bio-phys. J. , 95(12):5533–5543, 2008.[21] M. Marcondes de Freitas, E. Feliu, and C. Wiuf. Intermediates, catalysts, persistence,and boundary steady states.
J. Math. Biol. , 2016.[22] S. M¨uller, E. Feliu, G. Regensburger, C. Conradi, A. Shiu, and A. Dickenstein. Signconditions for injectivity of generalized polynomial maps with applications to chemicalreaction networks and real algebraic geometry.
Found. Comput. Math. , 16(1):69–97,2016.[23] M. Noro and K. Yokoyama. Verification of Gr¨obner basis candidates. In
MathematicalSoftware – ICMS 2014 , pages 419–424. Springer Berlin Heidelberg, 2014.[24] M. P´erez Mill´an and A. Dickenstein. The structure of MESSI biological systems. arXiv:1612.08763 , 2016.[25] M. P´erez Mill´an, A. Dickenstein, A. Shiu, and C. Conradi. Chemical reaction systemswith toric steady states.
Bull. Math. Biol. , 74(5):1027–1065, 2012.[26] L. Robbiano. Term Orderings on the Polynomial Ring.
Lecture Notes in ComputerScience , 204(93):513–517, 1985.[27] E. D. Sontag. Structure and stability of certain chemical networks and applicationsto the kinetic proofreading model of T-cell receptor signal transduction.
Institute ofElectrical and Electronics Engineers. Transactions on Automatic Control , 46(7):1028–1047, 2001.[28] M. Thomson and J. Gunawardena. The rational parameterization theorem for multi-site post-translational modification systems.
J. Theor. Biol. , 261:626–636, 2009.[29] F. Winkler. A p-adic approach to the computation of Gr¨obner bases.