A Linear Algebra Approach for Detecting Binomiality of Steady State Ideals of Reversible Chemical Reaction Networks
AA Linear Algebra Approach for DetectingBinomiality of Steady State Ideals of ReversibleChemical Reaction Networks
Hamid Rahkooy , Ovidiu Radulescu , and Thomas Sturm , CNRS, Inria, and the University of Lorraine, Nancy, France [email protected] LPHI CNRS UMR 5235, and University of Montpellier, Montpellier, France [email protected] MPI-INF and Saarland University, Saarbr¨ucken, Germany [email protected]
Abstract.
Motivated by problems from Chemical Reaction NetworkTheory, we investigate whether steady state ideals of reversible reactionnetworks are generated by binomials. We take an algebraic approachconsidering, besides concentrations of species, also rate constants as in-determinates. This leads us to the concept of unconditional binomiality,meaning binomiality for all values of the rate constants. This concept isdifferent from conditional binomiality that applies when rate constantvalues or relations among rate constants are given. We start by repre-senting the generators of a steady state ideal as sums of binomials, whichyields a corresponding coefficient matrix. On these grounds we proposean efficient algorithm for detecting unconditional binomiality. That al-gorithm uses exclusively elementary column and row operations on thecoefficient matrix. We prove asymptotic worst case upper bounds on thetime complexity of our algorithm. Furthermore, we experimentally com-pare its performance with other existing methods.
Keywords:
Binomial Ideals, Linear Algebra, Reversible Chemical Re-action Networks A chemical reaction is a transformation between two sets of chemical objectscalled chemical complexes . The objects that form a chemical complex are chemi-cal species . In other words, complexes are formal sums of chemical species repre-senting the left hand and the right hand sides of chemical reactions. A chemicalreaction network is a set of chemical reactions. For example CO + H k k CO + H O ,2 CO k k CO + C a r X i v : . [ c s . S C ] J un H. Rahkooy, O. Radulescu & T. Sturm is a chemical reaction network with two reversible reactions.A kinetics of a chemical reaction network is an assignment of a rate func-tion, depending on the concentrations of chemical species at the left hand side,to each reaction in the network. A kinetics for a chemical reaction network iscalled mass-action if for each reaction in the chemical reaction network, therate function is a monomial in the concentrations of the chemical species withexponents given by the numbers of molecules of the species consumed in thereaction, multiplied by a constant called rate constant. Reactions are classifiedas zero-order, first-order, etc. according to the order of the monomial giving therate. For reversible reactions, the net reaction rate is a binomial, the differencebetween the forward and backward rates. In the example above k , k , k , k are the rate constants. In this article we generally assume mass-action kinetics.We furthermore assume that reactions are reversible, unless explicitly specifiedotherwise.The change in the concentration of each species over time in a reaction canbe described via a system of autonomous ordinary differential equations. Forinstance, consider the chemical reaction network above and let x , x , x , x , x be the indeterminates representing the concentrations of the species CO , H , CO , H O and C , respectively. The corresponding differential equations are˙ x = p , p = − k x x + k x x − k x − k x x , (1)˙ x = p , p = − k x x + k x x , (2)˙ x = p , p = k x x − k x x + − k x + 2 k x x , (3)˙ x = p , p = k x x − k x x , (4)˙ x = p , p = − k x x + k x x + k x − k x x . (5)Each zero of the polynomials p , p , p , p , p gives a concentration of speciesin which the system is in equilibrium. The zeros of p , p , p , p , p are called thesteady states of the chemical reaction network. Accordingly, the ideal h p , p .p , p , p i ⊆ Q [ k , k , k , k , x , x , x , x , x ] is called the steady state ideal of the chemi-cal reaction network. We consider the coefficient field Q because of computabilityissue. Otherwise, theoretically, our results hold for any coefficient field. The so-lutions of these polynomials can be in R or in C .For a thorough introduction to chemical reaction network theory, we referto Feinberg’s Book [17] and his lecture notes [16]. We follow the notation ofFeinberg’s book in this article.An ideal is called binomial if it is generated by a set of binomials. In thisarticle we investigate whether the steady state ideal of a given chemical reac-tion network is binomial. We are interested in efficient algorithms for testingbinomiality. Consider the steady state ideal I = h p , p , p , p , p i ⊆ Q [ k , k , k , k , x , x , x , x , x ] , (6)given by Equations (1)–(5). Reducing p , p and p with respect to p and p ,we have I = h− k x x + k x x , − k x + k x x i , (7) itle Suppressed Due to Excessive Length 3 which shows that the ideal I is binomial. In this article, we work over the ring Q [ k ij , x , . . . , x n ] and investigate binomiality over this ring.Note that in the literature there exist also slightly different notions of bino-miality. Eisenbud and Sturmfels in [12] call an ideal binomial if it is generatedby polynomials with at most two terms. Following this definition, some au-thors, e.g., Dickenstein et al. in [32] have considered the steady state ideal as anideal in the ring Q ( k ij )[ x , . . . , x n ] and studied the binomiality of these ideals in R [ x , . . . , x n ] after specialising k ij with positive real values. In order to distin-guish between the two notions, we call unconditionally binomial a steady stateideal that is binomial in Q [ k, x ] (the notion used in this paper) and condition-ally binomial a steady state ideal that is binomial in Q ( k )[ x ], i.e. for specifiedparameters k (the notion used in [32]).The notions of binomial ideals and toric varieties have roots in thermody-namics, dating back to Boltzmann. Binomiality corresponds to detailed balance,which for reaction networks means that at thermodynamic equilibrium the for-ward and backward rates should be equal for all reactions. Detailed balance isa very important concept in thermodynamics, for instance it has been used byEinstein in his Nobel prize winning theory of the photoelectric effect [11], byWegscheider in his thermodynamic theory of chemical reaction networks [36]and by Onsager for deriving his famous reciprocity relations [33]. Because de-tailed balance implies time reversal symmetry, systems with detailed balance cannot produce directed movement and can only dissipate heat. This is importantin applications, for instance in molecular biology, where molecular motors cannot function with detailed balance. Although most interesting molecular devicesfunction without detailed balance and binomiality, some of their subsystems cansatisfy these conditions. The interest of studying binomiality relies in the sim-plicity of the analysis of such subsystems. For instance, important propertiessuch as multistationarity and stability are easier to establish for binomial sys-tems. Toricity, also known as complex, or cyclic, or semi-detailed balance is alsoknown since Boltzmann that has used it as a sufficient condition for deriving hisfamous H-theorem [1]. Binomiality implies toricity, but the converse is not true:in order to have binomiality, a toric system must obey constraints on the ratesconstants, such as the well known Weigscheider-Kolmogorov condition asking forthe equality of the products of forward and backward rates constants in cyclesof reversible reactions. In this paper we focus on the situation when detailedbalance is satisfied without conditions on the rate constants.Detecting binomiality of an ideal, particularly of a steady state ideal, is adifficult problem, both from a theoretical and a practical point of view. Theproblem is typically solved by computing a Gr¨obner basis, which is EXPSPACE-complete [29]. Recent linear algebra approaches for solving the problem in adifferent setting than our problem construct large matrices which also points atthe difficulty of the problem [30,6].There is quite comprehensive literature on chemical reaction network the-ory. An excellent reference to this topic is [17,16]. As mathematical concepts,binomiality and toricity have been widely studied and their properties have H. Rahkooy, O. Radulescu & T. Sturm been investigated by various authors, e.g., Fulton [18], Sturmfels [35], Eisenbudet al. [12]. Binomiality and toricity show up quite often in chemical reactionnetworks. Binomiality in the case of detailed balancing of reversible chemical re-actions has been studied by Gorban et al. [20,21] and Grigoriev and Weber [25].Feinberg [15] and Horn and Jackson [26] have studied toric dynamical systems.Gatermann et al. studied deformed toricity in [19]. Craciun, et al. have consid-ered the toricity problem over the real numbers in [7] and have presented severalinteresting results in this regard, among them, they have shown that complexbalanced systems are the same as toric dynamical systems, although toric steadystates are different from that. It has been shown in [9,10] that the binomialstructure will imply much simpler criteria for multistationarity. These resultsgive strong motivation for one to study algorithms for detecting binomial net-works. Especially, in [9], the authors defined linearly binomial network and theyproposed sufficient conditions for a network to be linearly binomial. The proofis constructive even though it has not been presented as an algorithm. Theirmethod is also quite straightforward and can handle more general networks inmany applications.Dickenstein et al. have presented sufficient linear algebra conditions with in-equalities for binomiality of the steady state ideals in [30]. Their idea has beendeveloped in [31], where the concept of MESSI reactions has been introduced.Conradi and Kahle have proved in [6] that for homogenous ideals (i.e. for chem-ical reaction networks without zero-order reactions), the sufficient condition ofDickenstein et al. is necessary as well and also introduced an algorithm for test-ing binomiality of homogenous ideals. As many biochemical networks are nothomogeneous, the algorithm requires heuristics in such cases. The algorithmhas been implemented in Maple and Macaulay II in [28,27] and experimentshave been carried out on several biological models. Grigoriev et al. in [23] haveconsidered the toricity of steady state ideals from a geometric point of view. In-troducing shifted toricity, they presented algorithms, complexity bounds as wellas experimental results for testing toricity using two important tools from sym-bolic computation, quantifier elimination [8,22,37] and Gr¨obner bases [4,5,13,14].Recently, first order logic test for toricity have been introduced [34].The main idea of this article is to consider the generators of the steadystate ideal as sums of the binomials associated to the reactions rather than themonomials associated to the complexes. This is feasible for a reversible chemicalreaction network. Following the above observation and assigning a binomial toeach reaction, one can write the generators of the steady state ideal as sums ofthose binomials with integer coefficients.As our main result, we have proved that a reversible chemical reaction net-work is unconditionally binomial if and only if it is “linearly” binomial (i.e., thereexist linear combinations of the generators such that these combinations are bi-nomials). More precisely, having represented of the generators of the steady stateideal as sum of binomials, one can test the binomiality exclusively using elemen-tary row and column operations on the coefficient matrix of these binomials. itle Suppressed Due to Excessive Length 5
This can be done by computing the reduced row echelon form of the coefficientmatrix, which yields an efficient method for testing binomiality.Our main contributions in this article are the following.1. We introduce a new representation of the generators of the steady state idealof a reversible chemical reaction as a sum of certain binomials rather thanmonomials.2. Using that representation, we assign a matrix with entries in Z to a reversiblechemical reaction network, such that the binomiality of the steady state idealcan be tested by computing the reduced row echelon form of this matrix.3. We prove a worst-case upper bound on the time complexity of our binomi-ality test. We experimentally compare our test with the existing binomialitytests in the literature, which demonstrates the applicability of our method.Our representation of the steady state ideal as a sum of certain binomials, aswell as the matrices associated to them are further original ideas presented in thispaper. While typically complex-species matrices are used for testing binomiality,we use reaction-species matrices for this purpose.The plan of the article is as follows. Section 1 gives an introduction to thenecessary concepts of chemical reaction network theory, reviews the literatureand presents the idea of this work. Section 2 includes the main definitions andresults. In this section we show our representation of the generators of the steadystate ideal of a reversible chemical reaction network and present our algorithmfor testing binomiality. In Section 3, we discuss the complexity of our method.We furthermore compare our algorithm with other existing algorithms in theliterature via experiments. In Section 4 we summarise our results and drawsome conclusions. In this section, we present our main result based on which we present an algo-rithm for testing unconditional binomiality of reversible chemical reaction net-works. In Subsection 2.1 we introduce a representation for the generators of thesteady state ideal of a chemical reaction network as sum of binomials. We showthat this representation is unique for reversible reaction networks, consideringrate constants as indeterminates. In Subsection 2.2, we define a matrix asso-ciated to a chemical reaction network which is essentially the species–reactionmatrix, rather than the stoichiometric matrix which is the species–complex ma-trix. Having considered constant rates as indeterminates, the uniqueness of ourmatrix for reversible reactions comes from the uniqueness of representing of thegenerators of the steady state ideal as sum of binomials.
Consider the following reversible reaction between two complexes C and C . H. Rahkooy, O. Radulescu & T. Sturm C k k C .Let m i , i = 1 ,
2, be the product of the concentrations of the species in C i with the stoichiometric coefficients as the powers. We call m i the monomialassociated to C i . Also let x be the concentration of a species that is in C with the stoichiometric coefficient α and is not in C . The differential equationdescribing the kinetics of this species is˙ x = − α k m + α k m . (8)For a species in C with stoichiometric coefficient α which is not in C with theconcentration x , the differential equation will be˙ x = α k m − α k m . (9)For a species with concentration x that appears in both C and C , the dif-ferential equation will be ˙ x = c ( k m − k m ), where c ∈ Z is the differ-ence between the corresponding stoichiometric coefficients in C and C . Set b := − k m + k m and b := k m − k m . The steady state ideal of theabove chemical reaction network is h α b , α b i , which is equal to h b i , since b = − b .For a reversible reaction network with more than one reaction, one can as-sociate a binomial of the form b ij := k ij m i − k ji m j to each reaction. Then thepolynomials generating the steady state ideal can be written as sums of b ij withinteger coefficients. We make this more precise in the following definition. Definition 1.
Let C be a reversible chemical reaction network with the com-plexes C , . . . , C s , let k ij , ≤ i = j ≤ s , be the rate constant of the reactionfrom C i to C j , and let x , . . . , x n be the concentrations of the species in the chem-ical reaction network. We call a monomial m i the monomial associated to C i if m i is the product of the concentrations of those species that appear in C i withthe stoichiometric coefficients of the species as the powers. If there is a reactionbetween C i and C j , then b ij := − k ij m i + k ji m j is called the binomial associatedto the reaction from C i to C j , otherwise b ij := 0 .Example 1. Recall the following chemical reaction network form Section 1: CO + H k k CO + H O ,2 CO k k CO + C .Following the notation in Section 1, let x , x , x , x , x be the concentrationsof CO , H , CO , H O and C , respectively. The monomials associated to thecomplexes CO + H , CO + H O , 2 CO and CO + C are x x , x x , x and x x , respectively. The binomials associated to the two reactions in this networkare b = − k x x + k x x and b = − k x + k x x . As there is noreaction between the first and third complexes we have b = b = 0. Similarly, b = b = 0, b = b = 0 and b = b = 0. Also, by definition, b = − b , itle Suppressed Due to Excessive Length 7 b = − b , etc.. Using the binomials associated to the reactions, one can writethe polynomials generating the steady state ideal as p = b − b (10) p = b (11) p = − b + 2 b (12) p = − b (13) p = − b . (14)Hence, the steady state ideal can be written as h p , p , p , p , p i = h b , b i . (15)As Example 1 and the definition of the binomials b ij in Definition 1 suggestsone can write the generators of the steady state ideal of every reversible chemicalreaction networks as sums of b ij with integer coefficients, i.e., assuming that R is the set of reactions in the chemical reaction network˙ x k = p k = X C i → C i ∈R c ( k ) ij b ij , (16)for k = 1 . . . n and c ( k ) ij ∈ Z .For clarification, we may remind the reader that in this article we assumeworking over Q [ k ij , x , . . . , x n ]. This is the case, in particular, for Definition 1and the discussion afterwards. In [32], the authors specialise k ij with positivereal values, in which case, the steady state ideal may or may not be binomial over R [ x , . . . , x n ]. Similarly, specialising k ij in Equation 16 can result in writing p k assum of different binomials. In other words, if k ij specialised, the representationof p k as sum of binomials in 16 is not necessarily unique. This is illustrated inthe following example. Example 2. [32, Example 2.3] Let C = 2 A , C = 2 B and C = A + B . Considerthe reversible chemical reaction network given by the following reactions:2 A k k B A k k A + BA + B k k B .Assuming x and x to be the concentrations of A and B , respectively, by Defi-nition 1, b = − k x + k x (17) b = − k x + k x x (18) b = k x − k x x . (19) H. Rahkooy, O. Radulescu & T. Sturm
It can be checked that the generators of the steady state ideal can be written as p = 2 b + b + b (20) p = − b − b − b . (21)If k = k then k x x = k x x , hence k x x will occur in b and b with opposite signs which will be cancelled out in b + b , resulting in writing p as sum of b and − k x + k x . This is another way of writing p as sumof binomials. Because binomiality relies here on the condition k = k , this isan example of conditional binomiality.If we consider the rate constants k ij as indeterminates, i.e., if we consider thesteady state ideal as an ideal over the ring Q [ k ij , x , . . . , x n ], then the represen-tation in Equation (16) as sum of binomials b ij will be unique. More precisely,we have the following. Lemma 1.
Given a reversible chemical reaction network with the notation ofDefinition 1, if k ij are indeterminates then the generators of the steady stateideal can be uniquely written as sum of the binomials presented in Equation 16.Proof. Assuming that k ij , 1 ≤ i, j ≤ s are indeterminates, they will be alge-braically independent over Q [ x , . . . , x n ]. Therefore for monomials m t and m t in Q [ x , . . . , x n ] associated to two distinct complexes and for all 1 ≤ i, j, i , j ≤ s , k ij m t and k i j m t will be distinct monomials in Q [ k ij , x , . . . , x n ]. Hence bino-mials b ij associated to the reversible reactions are not only pairwise distinct,but also their monomials are pairwise distinct in Q [ k ij , x , . . . , x n ]. This im-plies that the generators of the steady state ideal have unique representations in Q [ k ij , x , . . . , x n ] as sum of b ij with integer coefficients.Having a unique representation as in Equation 16 enables us to representour binomial coefficient matrix, defined later, which is the base of our efficientalgorithm for testing unconditional binomiality of reversible chemical reactionnetworks.Considering rate constants k ij as indeterminates, if a steady state ideal isunconditionally binomial, i.e., binomial in the ring Q [ k ij , x , . . . , x n ], then itselimination ideal is binomial in the ring Q [ x , . . . , x n ]. Indeed, the elimination ofa binomial ideal is a binomial ideal. This can be seen from elimination property ofGr¨obner bases. Authors of [12] have studied binomial ideals and their propertiesintensively. In particular Corollary 1.3 in the latter article state the binomiality ofthe elimination ideal of a binomial ideal. We remind the reader that the definitionof binomiality in this article is different from [12]. In the latter, binomial idealshave binomial and monomial generators, however in the current article, we onlyconsider binomial generators. Restricting the definition of binomial ideal to theideals with only binomial generators, most of the result in [12] still holds, inparticular the one about the elimination of binomial ideals. Therefore, if thesteady state ideal of a chemical reaction network is binomial in Q [ k ij , x , . . . , x n ],then its elimination I ∩ Q [ x , . . . , x n ] is also a binomial ideal. itle Suppressed Due to Excessive Length 9 Geometrically, the above discussion can be explained via projection of thecorresponding varieties. Given a chemical reaction network, assume that reactionrates k ij are indeterminates and let the number of k ij be t . Let V denote thesteady state variety, i.e., the variety of the steady state ideal. V is a Zariskiclosed subset of K t + n , where K is an appropriate field (e.g., C ). If V is a cosetof a subgroup of the multiplicative group ( K ∗ ) t + n , then the projection of V intothe space generated by x . . . , x n , i.e., V ∩ ( K ∗ ) n is also a coset. In particular,the projection of a group is a group. Since the variety of a binomial ideal is acoset [23,24], the projection of the variety of a binomial ideal is the variety of abinomial ideal. As special cases, the projection of a toric variety, a shifted toricvariety and a binomial variety (defined in [23,24]) is a toric, a shifted toric anda binomial variety, respectively. For a detailed study of toricity of steady statevarieties, we refer to [23]. Remark 1. – We may mention that in [7], the authors have studied toric dy-namical systems , where they have considered working over Q [ k ij , x , . . . , x n ]and presented several interesting results. In particular, Theorem 7 in thatarticle states that a chemical reaction network is toric if and only if the rateconstants lie in the variety of a certain ideal in Q [ k ij ], called the moduli ideal . – Toric dynamical systems are known as complex balancing mass action sys-tems [7].
Let C be a reversible chemical reaction network as in Definition 1and assume that the generators of its steady state ideal are written as the linearcombination of the binomials associated to its reactions as in Equation 16, i.e., p k = s X C i → C i ∈R c ( k ) ij b ij for k = 1 , . . . , n. We define the binomial coefficient matrix of C to be the matrix whose rows arelabeled by p , . . . , p n and whose columns are labeled by non-zero b ij and the entryin row p k and column b ij is c ( k ) ij ∈ Z . By the definition, the binomial coefficient matrix of a reversible chemicalreaction network is the coefficient matrix of the binomials that occur in therepresentation of the generators of the steady state ideal as sum of binomials.As we consider k ij indeterminates, the representation of the generators of thesteady state ideal of a given complex is unique, which implies that the binomialcoefficient matrix of a given complex is unique too. Example 3.
Consider the chemical reaction network in Example 1, with genera-tors of the steady state ideal as follows. p = b − b (22) p = b (23) p = − b + 2 b (24) p = − b (25) p = − b . (26)The binomial coefficient matrix of this chemical reaction network is M = b b p − p p − p − p − . (27)Another simple example is the reaction4 A k k A + B ,with the binomial associated to it as b := − k x + k x x , where x is theconcentration of A and x is the concentration of B . The steady state ideal isgenerated by { b, − b } , and the binomial coefficient matrix for this network is (cid:0) − (cid:1) .One can test binomiality of the steady state ideal of a reversible reactionnetwork using its binomial coefficient matrix. Theorem 1.
The steady state ideal of a reversible chemical reaction networkis unconditionally binomial, i.e., binomial in Q [ k ij , x , . . . , x n ] , if and only ifthe reduced row echelon form of its binomial coefficient matrix has at most onenon-zero entry at each row.Proof. Let G = { p , . . . , p n } ⊆ Q [ k ij , x , . . . , x n ] be a generating set for thesteady state ideal of a given reversible chemical reaction network C , and let { b ij | ≤ i = j ≤ s } be the ordered set of non-zero binomials associated to thereactions. Fix a term order on the monomials in Q [ k ij , x , . . . , x n ].First we prove that if the reduced row echelon form of the binomial coefficientmatrix has at most one non-zero entry at each row, then the steady state ideal isbinomial. The proof of this side of the proposition comes from the definition ofreduced row echelon form. In fact, the reduced row echelon form of the binomialcoefficient matrix of C can be computed by row reduction in that matrix, whichis equivalent to the reduction of the generators of the steady state ideal withrespect to each other. Therefore, computing the reduced row echelon form ofthe binomial coefficient matrix and multiplying it with the vector of binomials b ij , one can obtain another basis for the steady state ideal. Having this, if the itle Suppressed Due to Excessive Length 11 reduced row echelon form has at most one non-zero entry at each row, then thenew basis for the steady state ideal will only include b ij . Therefore the steadystate ideal will be binomial.Now we prove the “only if” part of the proposition, that is, if the steadystate ideal of C is binomial, then the reduced row echelon form of the binomialcoefficient matrix has at most one non-zero entry at each row. We claim that foreach pair of polynomials p t , p m ∈ G , p t is reducible with respect to p m if andonly if there exists a binomial b ij that occurs in both p t and p m and includestheir leading terms. The “only if” part of the claim is obvious. To prove the “if”part of the claim, let p m be reducible with respect to p t . Then the leading termof p m divides the leading term of p t . Since the leading terms are multiples of k ij and these are disjoint indeterminates, this is only possible if both of the leadingterms are equal. If the leading terms are equal, then b ij in which the leadingterms occur, must itself occur in both p t and p m . Therefore p t and p m share abinomial associated to a reaction, which is in contradiction with our assumption.From the above claim and the definition of the reduced row echelon formone can see that p , . . . , p n are pairwise irreducible if and only if the binomialcoefficient matrix of C is in reduced row echelon form.Now we prove that p , . . . , p n are pairwise irreducible if and only if they forma Gr¨obner basis in which polynomials are pairwise irreducible. Note that thisdoes not necessarily imply that G is a a reduced Gr¨obner basis, as p i are notnecessarily monic. Assume that p , . . . , p n are pairwise irreducible. We prove thatthe greatest common divisor of each pair of the leading terms of the p , . . . , p n is1. By contradiction, assume that there exists a monomial not equal to 1 whichdivides the leading terms of both p t , p m , for 1 ≤ t, m ≤ n . Then there exists avariable x l such that x l divides the leading terms of p t and p m . Since each leadingterm is the monomial associated to a complex, the species with concentration x occurs in two complexes with associated monomials as the leading terms of p t and p m . Then both p t and p m have as their summand the binomials that areassociated to the reactions including those complexes. As for each complex thereexists at least one binomial associated, both p m and p t have as a summand onecommon binomial b ij . However, we had already proved that this implies that p t and p m are not pairwise irreducible, which is a contradiction to the assumptionthat the greatest common divisor of the leading terms of p t and p m is not 1.Now by Buchberger’s first criterion if the greatest common divisor of the leadingterms of each pair of polynomials in G is 1 then G is a Gr¨obner basis. The otherside of this claim is obvious.From what we have proved until now, we can conclude that the binomialcoefficient matrix of C is in reduced row echelon form if and only if G is aGr¨obner basis with pairwise irreducible elements. On the other hand, by a resultof Eisenbud and Sturmfels [12], the steady state ideal of C is binomial if andonly if every Gr¨obner basis of it includes binomials. Therefore we conclude thatthe steady state ideal is binomial if and only if the reduced row echelon form ofthe binomial coefficient matrix has at most one non-zero entry in each row. Function
BinomialityTest( C ) Input: C = { ( m , . . . , m s ) ∈ [ X ] n , k ij } Output:
Binomial or NotBinomial b ij := − k ij m i + k ji m j , ≤ i = j ≤ s B := ( b ij , ≤ i = j ≤ m ) p k := P c kij b ij , ≤ k ≤ n M := Matrix( c kij ) ˜ M = ReducedRowEchelonForm( M ) G := ˜ MB if IsBinomial( G ) then R := Binomial else R := NotBinomial return R Algorithm 1:
Testing Unconditional Binomiality of Reversible Chemical Re-action Networks
Example 4.
Following Example 3, one case easily see that the reduced row ech-elon form of the binomial coefficient matrix (27) is M = b b p p p p p , (28)which means that the steady state ideal is unconditionally binomial and is gen-erated by { b , b } .Theorem 1 yields Algorithm 1 for testing unconditional binomiality. Theinput of the algorithm is a reversible chemical reaction network, given by thevector of monomials associated to its complexes, ( m , . . . , m s ), and the rates k ij .It uses a function IsBinomial which takes a set of polynomials and checks if allof them are binomial.
Generalisation to Non-Reversible Networks
The unconditional binomiality testvia the binomial coefficient matrix for a reversible chemical reaction networkcan be used as a subroutine for testing unconditional binomiality of an arbitrarychemical reaction network. In order to do so, partition a given chemical reactionnetwork C into a reversible reaction network C and a non–reversible reactionnetwork C . Apply Algorithm 1 to C , construct its binomial coefficient matrix,say M . Construct the stoichiometric coefficient matrix of C , say M , and con-sider the block matrix M := ( ˜ M | M ). Compute the row reduced echelon formof M , say ˜ M . If all the rows of ˜ M have at most one non-zero entry, then thesteady state ideal is binomial. itle Suppressed Due to Excessive Length 13 Otherwise, one can consider computing ˜ M as a preprocessing step and runanother method, e.g., Gr¨obner bases, quantifier elimination as in [23], or themethod in Dickenstein, et al [30]. Proposition 1.
Let r be the number of reactions and n be the number of speciesof a reversible chemical reaction network C . The asymptotic worst case timecomplexity of testing unconditional binomiality of the steady state ideal of C viaAlgorithm 1 can be bounded by O (max( r, n ) ω ) where ω ≈ . , which is alsothe complexity of matrix multiplication.Proof. The operations in steps 1–4 and 7–11 are at most linear in terms of r and n . Since M is a matrix of size n × r , where r = | b ij | , and B is a vectorof size r , computing reduced row echelon form in step 5 and also the matrixmultiplication in step 6 will cost at most O (max( r, n ) ω ). Therefore the totalnumber of operations in the algorithm can be bounded by O (max( r, n ) ω ).In [23, Section 4] it has been shown that there exists an exponential asymp-totic worst case upper bound on the time complexity of testing toricity. An imme-diate consequence of that result is that the time complexity of testing binomialitycan be bounded by the same exponential function. Following the arguments in[23, Section 4], one can show that there exists an algorithm for testing binomial-ity over Q [ k ij , x , . . . , x n ] and Q [ x , . . . , x n ] simultaneously, with an exponentialupper bound for the worst case time complexity.As mentioned earlier in Section 2, the reduced Gr¨obner basis of a binomialideal, with respect to every term order, includes only binomials. This directly canbe seen from running Buchberger’s algorithm and that S–polynomials and theirreductions by binomials are binomial. Eisenbud et. al’ article [12], with a slightlydifferent definition of binomial ideals, investigates many properties of binomialideals using the latter fact. Following this fact, a typical method for testingbinomiality is via computing a reduced Gr¨obner basis of a steady state ideal I ⊆ Q [ k ij , x , . . . , x n ] The drawback of computing Gr¨obner bases is that thisis EXPSPACE-complete [29]. So our algorithm is asymptotically considerablymore efficient than Gr¨obner basis computation. Example 5 (Models from the BioModels Repository ). – There are twenty non–reversible biomodels in which Gr¨obner basis compu-tations done in [23] for testing conditional binomiality do not terminate ina six–hour time limit, however our algorithm terminates in less than threeseconds. Also there are six cases in which Gr¨obner basis computations ter-minate in less than six hours, but are at least 1000 times slower than ouralgorithm. Finally there are ten models in which Gr¨obner basis is at least500 times slower than our computations. – There are sixty nine biomodels that are not considered for computation in[23] because of the of the unclear numeric value of their rate constants. Ourcomputations on almost all of those cases terminated in less than a second. – (Reversible models from the BioModels Repository) Biomodels 491 and 492are both reversible. Biomodel 491 has 52 species and 86 reactions. The bino-mial coefficient matrix of this biomodel has size 52 ×
86 and has ± .
344 secondsthat it is unconditionally binomial, while a Gr¨obner basis computation takesmore than 12 seconds to check its conditional binomiality. BioModel 492 hasalso 52 species, and includes 88 reactions. The binomial coefficient matrixhas entries ± ×
88. This biomodel is also unconditionallybinomial. It takes 0 .
25 seconds for Maple to check its unconditional bino-miality via Algorithm 1 in Maple, while a Gr¨obner basis computation takesnear 18 seconds, as one can see in the computations in [23, Table 3], whichshow the group structure of the steady state varieties of the models.Dickenstein et al. in [30] have proposed a method for testing toricity ofa chemical reaction network. The definitions and purpose of that work areslightly different from our article, hence comparisons between those two meth-ods should be treated with caution. While we focus on unconditional binomialityof the steady state ideals of reversible reaction networks, , i.e., binomiality in Q [ k ij , x , . . . , x n ], with the aim of efficiency of the computations, the authors ofthe above article are interested in conditional binomiality with algebraic depen-dencies between k ij such that the elimination ideal is binomial. Having men-tioned that, our method leads to the computation of reduced row echelon formof a matrix of size n × r with integer entries which is polynomial time, whileTheorem 3.3. in [30] requires constructing a matrix of size n × s with entriesfrom Z [ k ij ] and finding a particular partition of its kernel.Considering Example 2.3 in [32], our algorithm constructs the matrix M andits reduced row echelon form ˜ M : M = (cid:18) − − − (cid:19) , ˜ M = (cid:18) (cid:19) , (29)and we see that the steady state ideal is not unconditionally binomial over Q [ k ij , x , . . . , x n ]. The method in [32] constructs (cid:18) − k − k k + k k − k k + k − k − k − k + k (cid:19) , (30)and finds an appropriate partition, which shows that the steady state ideal isbinomial in Q [ x , . . . , x n ] if and only if k = k . As a larger example, considerthe chemical reaction network given in Example 3.13 in [32] and assume that itis a reversible chemical reaction network. Our method constructs a matrix withentries ± × ± ×
10 matrix with entries as linear polynomials in Z [ k ij ] and computes aparticular partition of the kernel of the matrix. itle Suppressed Due to Excessive Length 15 For homogeneous ideals, Conradi and Kahle have shown in [6] that the suf-ficient condition for conditional binomiality in [32] is necessary, too. Their Al-gorithm 3.3 tests conditional binomiality of a homogeneous ideal, which can begeneralised by homogenising. The algorithm computes a basis for the ideal degreeby degree and performs reductions with respect to the computed basis elementsat each degree step. Since our algorithm is intended for steady state ideals ofreversible chemical reaction networks, which are not necessarily homogeneous,our following comparison with the Conradi–Kahle algorithm bears a risk of be-ing biased by homogenisation. We discuss the execution of both algorithms onExample 3.15 in [32]. This chemical reaction network does not satisfy the suf-ficient condition presented in [32, Theorem 3.3]. Testing this condition leads tothe construction of a 9 ×
13 matrix with entries in Z [ k ij ], followed by furthercomputations, including finding a particular partition of its kernel. Theorem3.19 in [32] is a generalisation of Theorem 3.3 there, which can test conditionalbinomiality of this example by adding further rows and columns to the matrix.Conradi and Kahle also treat this example with their algorithm. This requiresthe construction of a coefficient matrix of size 9 ×
13 with entries in Z [ k ij ] andcertain row reductions. If we add reactions so that the reaction network becomesreversible, our algorithm will construct a matrix of size 9 × ± Q [ k ij , x , . . . , x ]. Binomiality of steady state ideals is an interesting problem in chemical reac-tion network theory. It has a rich history and literature and is still an activeresearch area. For instance, recently MESSI systems have been introduced [31]following the authors’ work on binomiality of a system. Finding binomiality andtoricity is computationally hard from both a theoretical and a practical point ofview. It typically involves computations of Gr¨obner bases, which is EXPSPACE-complete.In a recent work [23] we investigated toricity of steady state varieties andgave efficient algorithms. In particular, we experimentally investigated toricityof biological models systematically via quantifier elimination. Besides that, wepresented exponential theoretical bounds on the toricity problem. The currentarticle, restricting to reversible reaction networks, aims at an efficient linearalgebra approach to the problem of unconditional binomiality, which can beconsidered as a special case of the toricity problem.In that course, considering rate constants as indeterminates, we assign aunique binomial to each reaction and construct the coefficient matrix with re-spect to these binomials. Our algorithm proposed here computes a reduced rowechelon form of this matrix in order to detect unconditional binomiality. Thealgorithm is quite efficient, as it constructs comparatively small matrices whoseentries are integers. It is a polynomial time algorithm in terms of the number ofspecies and reactions. While other existing methods for testing conditional bino- miality have different settings and purposes than our algorithm, for the commoncases, our algorithm has advantages in terms of efficiency.
Acknowledgments
This work has been supported by the bilateral project ANR-17-CE40-0036 andDFG-391322026 SYMBIONT [2,3].
References
1. Ludwig Boltzmann.
Lectures on Gas Theory . University of California Press, Berke-ley and Los Angeles, CA, 1964.2. Fran¸cois Boulier, Fran¸cois Fages, Ovidiu Radulescu, Satya Samal, Andreas Schup-pert, Werner Seiler, Thomas Sturm, Sebastian Walcher, and Andreas Weber. TheSYMBIONT project: Symbolic methods for biological networks.
ACM Communi-cations in Computer Algebra , 52(3):67–70, 2018.3. Fran¸cois Boulier, Fran¸cois Fages, Ovidiu Radulescu, Satya Samal, Andreas Schup-pert, Werner Seiler, Thomas Sturm, Sebastian Walcher, and Andreas Weber. TheSYMBIONT project: Symbolic methods for biological networks.
F1000Research ,7(1341), 2018.4. Bruno Buchberger.
Ein Algorithmus zum Auffinden der Basiselemente des Rest-klassenringes nach einem nulldimensionalen Polynomideal . Doctoral dissertation,Mathematical Institute, University of Innsbruck, Austria, 1965.5. Bruno Buchberger. Ein Algorithmisches Kriterium f¨ur die L¨osbarkeit eines alge-braischen Gleichungssystems.
Aequationes Mathematicae , 3:374–383, 1970.6. Carsten Conradi and Thomas Kahle. Detecting binomiality.
Advances in AppliedMathematics , 71:52–67, 2015.7. Gheorghe Craciun, Alicia Dickenstein, Anne Shiu, and Bernd Sturmfels. Toricdynamical systems.
J. Symb. Comput. , 44(11):1551–1565, November 2009.8. James H. Davenport and Joos Heintz. Real quantifier elimination is doubly expo-nential.
J. Symb. Comput. , 5(1–2):29–35, 1988.9. Alicia Dickenstein, Mercedes P´erez Mill´an, Anne Shiu and Xiaoxian Tang. Mul-tistatonarity in Structured Reaction Networks.
Bulletin of Mathematical Biology .81: 1527–1581, 2019.10. AmirHossein Sadeghimanesh and Elisenda Feliu. The multistationarity structure ofnetworks with intermediates and a binomial core network.
Bulletin of MathematicalBiology . 81: 2428–2462, 2019.11. Albert Einstein. Strahlungs-emission und -absorption nach der Quantentheorie.
Verh. Dtsch. Phys. Ges. , 18:318–323, 1916.12. David Eisenbud and Bernd Sturmfels. Binomial ideals.
Duke Math. J. , 84(1):1–45,1996.13. Jean-Charles Faug`ere. A new efficient algorithm for computing Gr¨obner bases(F4).
J. Pure Appl. Algebra , 139(1):61–88, 1999.14. Jean-Charles Faug`ere. A new efficient algorithm for computing gr¨obner baseswithout reduction to zero (F5). In
Proc. ISSAC 2002 , pages 75–83. ACM, NewYork, NY, 2002.15. Martin Feinberg. Complex balancing in general kinetic systems.
Arch. Ration.Mech. An. , 49(3):187–194, 1972.itle Suppressed Due to Excessive Length 1716. Martin Feinberg. Lectures on chemical reaction networks, 1979.17. Martin Feinberg.
Foundations of Chemical Reaction Network Theory , volume 202of
Applied Mathematical Sciences . Springer, 2019.18. William Fulton.
Introduction to Toric Varieties , volume 131 of
Annals of Mathe-matics Studies . Princeton University Press, 1993.19. Karin Gatermann. Counting stable solutions of sparse polynomial systems inchemistry. In
Symbolic Computation: Solving Equations in Algebra, Geometry,and Engineering , volume 286 of
Contemporary Mathematics , pages 53–69. AMS,Providence, RI, 2001.20. A. N. Gorban and V. N. Kolokoltsov. Generalized mass action law and thermody-namics of nonlinear Markov processes.
Math. Model. Nat. Phenom. , 10(5):16–46,2015.21. A. N. Gorban and G. S. Yablonsky. Three waves of chemical dynamics.
Math.Model. Nat. Phenom. , 10(5):1–5, 2015.22. D. Yu. Grigoriev. Complexity of deciding Tarski algebra.
J. Symb. Comput. ,5(1–2):65–108, 1988.23. Dima Grigoriev, Alexandru Iosif, Hamid Rahkooy, Thomas Sturm, and AndreasWeber. Efficiently and effectively recognizing toricity of steady state varieties.
CoRR , abs/1910.04100, 2019.24. Dima Grigoriev and Pierre D. Milman. Nash resolution for binomial varieties asEuclidean division. A priori termination bound, polynomial complexity in essentialdimension 2.
Adv. Math. , 231(6):3389–3428, 2012.25. Dima Grigoriev and Andreas Weber. Complexity of solving systems with fewindependent monomials and applications to mass-action kinetics. In
Proc. CASC2012 , volume 7442 of
LNCS , pages 143–154. Springer, 2012.26. F. Horn and R. Jackson. General mass action kinetics.
Arch. Ration. Mech. An. ,47(2):81–116, 1972.27. Alexandru Iosif and Hamid Rahkooy. Analysis of the Conradi–Kahle algorithm fordetecting binomiality on biological models.
CoRR , abs/1912.06896, 2019.28. Alexandru Iosif and Hamid Rahkooy. MapleBinomials, a Maple package for testingbinomiality of ideals. http://doi.org/10.5281/zenodo.3564428 , 2019.29. Ernst W. Mayr and Albert R. Meyer. The complexity of the word problems forcommutative semigroups and polynomial ideals.
Adv. Math. , 46(3):305–329, 1982.30. Mercedes P´erez Mill´an, Alicia Dickenstein, Anne Shiu, and Carsten Conradi.Chemical reaction systems with toric steady states.
Bulletin of Mathematical Bi-ology , 74(5):1027–1065, 2012.31. Mercedes P´erez Mill´an and Alicia Dickenstein. The structure of MESSI biologicalsystems.
SIAM J. Appl. Dyn. Syst. , 17(2):1650–1682, 2018.32. Mercedes P´erez Mill´an, Alicia Dickenstein, Anne Shiu, and Carsten Conradi.Chemical reaction systems with toric steady states.
Bull. Math. Biol. , 74(5):1027–1065, 2012.33. Lars Onsager. Reciprocal relations in irreversible processes. I.
Phys. Rev. ,37(4):405, 1931.34. Hamid Rahkooy and Thomas Sturm. First-order tests for toricity.
CoRR ,abs/2002.03586, 2020.35. Bernd Sturmfels.
Gr¨obner Bases and Convex Polytopes , volume 8 of
UniversityLecture Series . AMS, Providence, RI, 1996.36. Rudolf Wegscheider. ¨Uber simultane Gleichgewichte und die Beziehungen zwis-chen Thermodynamik und Reactionskinetik homogener Systeme.
Monatsh. Chem.Verw. Tl. , 22(8):849–906, 1901.8 H. Rahkooy, O. Radulescu & T. Sturm37. Volker Weispfenning. The complexity of linear problems in fields.