A Graph Theoretical Approach for Testing Binomiality of Reversible Chemical Reaction Networks
AA Graph Theoretical Approach for Testing Binomiality ofReversible Chemical Reaction Networks
Hamid Rahkooy Cristian Vargas Montero
CNRS, Inria, and the University of Lorraine CNRS, Inria, and the University of LorraineNancy, France Nancy, [email protected] [email protected]
Abstract
We study binomiality of the steady state idealsof chemical reaction networks. Considering rateconstants as indeterminates, the concept of un-conditional binomiality has been introduced andan algorithm based on linear algebra has beenproposed in a recent work for reversible chem-ical reaction networks, which has a polynomialtime complexity upper bound on the numberof species and reactions. In this article, usinga modified version of species–reaction graphs,we present an algorithm based on graph the-ory which performs by adding and deleting edgesand changing the labels of the edges in order totest unconditional binomiality. We have imple-mented our graph theoretical algorithm as wellas the linear algebra one in Maple and made ex-periments on biochemical models. Our experi-ments show that the performance of the graphtheoretical approach is similar to or better thanthe linear algebra approach, while it is drasticallyfaster than Gr¨obner basis and quantifier elimina-tion methods.
Chemical reactions are transformations betweenchemical species where a creation or eliminationof species may happen with respect to changes intime, pressure, temperature, etc. A set of chem-ical reactions is called a chemical reaction net- work (CRN) and if all the reactions in a CRN arereversible, it is called a reversible chemical reac-tion network (RCRN). We assume mass–actionkinetics in this article. The following is an ex-ample of an RCRN.A + B k k C k k A + 2 D,with species A, B, C and D and rate constants k , k , k and k .Ordinary differential equations (ODE) can beused to study the changes in the concentrationof species of a CRN. The ODEs associated to theabove CRN are˙ x = p , p = − k x x + k x + k x − k x x ˙ x = p , p = − k x x + k x ˙ x = p , p = k x x − k x − k x + k x x ˙ x = p , p = 2 k x − k x x (1)The ideal generated by p , p .p , p is calledthe steady state ideal of the CRN. The real posi-tive zeros of the above ideal are called the steadystates . Finding steady states is a fundamentalproblem in CRN theory. For a thorough intro-duction to CRN theory, we refer to Feinberg’sBook [12].A CRN is called binomial if its steady stateideal is binomial. Following the work in [26], in1 a r X i v : . [ c s . S C ] O c t his article we investigate binomiality over thering Q [ k ij , x , . . . , x n ], which we call uncondi-tional binomiality . Some authors have consid-ered binomiality of CRNs over Q ( k ij )[ x , . . . , x n ]when k ij are specialised in R , e.g., in [24], whichwe call conditional binomiality .Binomial ideals and toric varieties are his-toric topics in thermodynamic and go back toBoltzmann and Einstein. Binomiality and toric-ity have been widely studied in mathematics[13, 29, 10]. Binomiality is a hard problem andthe typical approach by computing a Gr¨obnerbasis is EXPSPACE-complete [21].In CRN theory, various articles have beenwritten on binomiality and toricity, e.g., by Gor-ban et al. [14] and Grigoriev and Weber [17]and also other authors [9, 27]. Feinberg [11] andCraciun, et al. [7] have studied toric dynamicalsystems.Dickenstein et al. have presented sufficient lin-ear algebra conditions with inequalities for con-ditional binomiality in [24], which lead to MESSICRNs [23]. For homogeneous ideals, it has beenshown that Dickenstein et al.’s condition is nec-essary as well [6]. The latter has been imple-mented in Maple and Macaulay II in [19, 18].A geometric view towards toricity has been con-sidered by Grigoriev et al. in [16], introducingshifted toricity and presenting algorithms usingquantifier elimination [8, 15, 30] and Gr¨obnerbases [3, 4] A first order logic test for toricityhas been given in [25]. In [26] unconditional bi-nomiality has been introduced and for reversiblereactions a polynomial time algorithm, based ona linear algebra approach, has been presented.Graph theory has been used in the study ofCRNs, e.g., for detecting concordance and weakmonotonicity [12, 28].The main contribution of this article is a graphtheoretical approach for testing unconditional bi-nomiality of an RCRN. We use a modificationof species–reaction graphs and present an algo-rithm equivalent to the linear algebra approachpresented in [26]. We have implemented thegraph theoretical algorithm as well the linearalgebra algorithm in Maple [20] and compared them with each other and with the algorithmsbased on Gr¨obner basis and quantifier elimina-tion presented in [16] via experiments on biolog-ical models from the BioModels repository [5].Our experiments showed that the graph theoret-ical and linear algebra approaches are not onlydrastically faster, but they can also handle manycases that were not possible to compute usingGr¨obner basis and quantifier elimination.The plan of this article is as follows. In Section2 we briefly review the linear algebra approachin [26]. Section 3 presents the graph theoreticalalgorithm and the proof of its correctness, whichare the main contributions of this article. Imple-mentations of and experiments using both thegraph theoretical and linear algebra algorithmsare presented in Section 4. In this section we briefly review the linear al-gebra method in [26] for testing unconditionalbinomiality in RCRNs.
Definition 1 (Definition 1, [26]) . Let C be a re-versible chemical reaction network with the com-plexes C , . . . ,C s , let k ij , ≤ i (cid:54) = j ≤ s, bethe constant rate of the reaction from C i to C j ,and let x , · · · , x n be the concentrations of thespecies. b ij := − k ij m i + k ji m j is called the bi-nomial associated to the reaction from C i to C j .If there is no reaction between C i and C j , set b ij := 0. Example 1.
The binomials associated to thereversible reactions presented in the introductionsection are b = − k x x + k x and b = − k x + k x x .Following Definition 1, one can write the cor-responding ODEs for an RCRN as˙ x k = p k = s (cid:88) j =1 c ( k ) ij b ij , k = 1 , . . . , n, (2)2here c ( k ) ij is the difference between the stoichio-metric coefficients of the k − th species in the re-action C i C j . It has been shown that if rateconstants are indeterminates then the above sumof binomials representation is unique [26] . Definition 2 (Definition 6, [26]) . Let C be anRCRN as in Definition 1 and assume that p i , thegenerators of its steady state ideal, are writtenas sum of binomials as in Equation 2. We definethe binomial coefficient matrix of C to be thematrix whose rows are labeled by p , . . . , p n andwhose columns are labeled by nonzero b ij andthe entry in row p k and column b ij is c ( k ) ij .The binomial coefficient matrix can be used totest the unconditional binomiality of an RCRN. Theorem 3 (Theorem 8, [26]) . A RCRN is un-conditionally binomial (i.e., assuming the rateconstants to be indeterminates) if and only if thereduced row echelon form(RREF) of its binomialcoefficient matrix has at most one nonzero entryat each row.
The theorem leads to a linear algebra approachfor testing unconditional binomiality [26, Algo-rithm 1], which is polynomial time in terms ofthe size of matrix, which is the maximum of thenumber of species and reactions in the RCRN.In Section 4, we present an implementation ofthe algorithm in Maple and compare it with theother approaches.
In this section we present a graph theoreticalalgorithm equivalent to the linear algebra ap-proach in [26]. A CRN can be represented asa graph in several ways. A first idea is thewell–known complex–reaction graph of a CRN,in which the complexes are considered as thevertices and the reactions as the directed edges.Another idea is species–reaction graphs , whichis used to study concordance and weakly mono-tonicity of kinetics in a CRN [12, Theorem Figure 1: Mod. Spec–Reac. Graph (Example 2)10.5.5, Theorem 11.5.1, Theorem 11.6.1]. We usea modified version of the species–reaction graph,defined only for RCRNs, for detecting uncondi-tional binomiality. For definition and notationson graph theory we rely on [31].
Definition 4 (Modified Species–ReactionGraph of an RCRN) . Let S be the set of speciesand R the set of reactions of a given RCRN.The modified species–reaction graph G of theRCRN is defined as follows. • For each species s ∈ S , consider a vertex of G (species vertices) • For each reaction r ∈ R , consider a vertexof G (reactions vertices) • For each reaction vertex, there exists anundirected edge to the species vertices of thespecies that appear in the reaction. • Each edge of the graph is labeled by aninteger number which is the difference be-tween the stoichiometric coefficients of thespecies(present at one end of the edge) inthe reactant and product complexes.
Example 2.
The species–reaction graph of thefollowing RCRN is showed in Figure 1.2 A + B C A 2 B.In order to check unconditional binomiality ofan RCRN using the modified species–reaction3raph, we present an algorithm that modifies thegraph by adding and deleting the edges, and up-dating the adge labelsso that unconditional bi-nomiality can be easily read off from the finalgraph. The algorithm simulates the proceduredescribed in Theorem 3 from Section 2 and Al-gorithm 1 from [26] by reducing a binomial co-efficient matrix to its RREF in order to test un-conditional binomiality of the network.The idea of the algorithm is as follows. First,we iterate through the reaction vertices, select-ing and marking an (random) unmarked speciesvertex connected to the current reaction vertex(initially all species vertices are unmarked). Ifthere are no unmarked species vertices, then thecurrent reaction vertex is skipped. Then, wedelete all the edges incident to the current reac-tion vertex that are different from the edge go-ing to the marked species vertex. Furthermore,if an edge exists from the marked species vertexto a reaction vertex, then we add the edge fromthe reaction vertex to the current species vertex.Nevertheless, if the edges already exists, then weupdate the label accordingly and if the new la-bel is zero then we eliminate the edge. The finalgraph is reached when all reaction vertices havebeen visited. At the end of the algorithm, uncon-ditional binomiality is checked by testing if eachcomponent of the final graph contains either onlyone species vertex or one species vertex and onereaction vertex. The algorithm is fully describedbelow.The functions in Algorithm 1 are self-explanatory. We just mention that
ElimEdge eliminates the edge connecting a species vertexand reaction vertex, and
UpdCf updates the co-efficient of the edge that goes from a species ver-tex to a reaction vertex (Detailed description ofthe functions can be found in the Git repository[22]).
Theorem 5.
Algorithm 1 is correct, terminatesand its asymptotic worst case time complexitycan be bounded by O (max( r, n ) ω ) , where ω is theconstant appearing in the complexity of matrixmultiplication, r is the number of reactions and n is the number of species of an RCRN. Function
BinomialityTestViaGraph( S , R , G ) Input: S : set of species of the RCRN R : set of reactions of the RCRN Output:
UnconditionallyBinomial orNotUnconditionallyBinomial G := CreateGraph( S , R ) SV := GetSpeciesVertices( G ) RV := GetReactionVertices( G ) foreach s ∈ SV do SetMark( s,false ) foreach r ∈ RV do speciesF ound := false cs := null rSpecies := GetConnectedSpecies( G , r ) foreach sr ∈ rSpecies do if IsNotMarked( sr ) then speciesF ound := true cs := sr break if speciesF ound then SetMark( cs,true ) foreach s (cid:48) ∈ rSpecies do if s (cid:48) (cid:54) = cs then multX := − GetCoeff( G , s (cid:48) , r )GetCoeff( G ,cs,r ) ElimEdge( G , s (cid:48) , r ) sReactions := GetConnectedReactions( G , cs ) foreach r (cid:48) ∈ sReactions do if r (cid:48) (cid:54) = r then if AreConnected( G , s (cid:48) , r (cid:48) ) then cf := GetCoeff( G , cs, r (cid:48) ) ∗ multiX + GetCoeff( G , s (cid:48) , r (cid:48) ) if cf (cid:54) = 0 then UpdCf( G , s (cid:48) , r (cid:48) , cf ) else ElimEdge( G , s (cid:48) , r (cid:48) ) else cf := GetCoeff( G , cs, r (cid:48) ) ∗ multiX AddEdge( G , s (cid:48) , r (cid:48) , cf ) if IsUnconditionallyBinomial( G ) then R := UnconditionallyBinomial else R := NotUnconditionallyBinomial return R Algorithm 1:
Testing unconditional binomi-ality via graphs4 roof. (Proof of the correctness)
Assumethat the algorithm output is unconditionally bi-nomial, then we must prove that the steady stateideal of the RCRN is binomial. To do so, wewill show that the steps of the graph theoreticalAlgorithm 1 are equivalent to the steps of thelinear algebra approach in [26, Algorithm 1] us-ing Reduced Row Echelon Form (RREF) in thebinomial coefficient matrix. Based on Theorem3, Algorithm 1 from [26] initially selects a pivotin the matrix which is equivalent to marking anunmarked species vertex that is connected to areaction vertex in steps 9 to 14. Reducing thenonzero entries to zero in the column of the se-lected pivot in the matrix is equivalent to elimi-nating the edges from the reaction vertex to allother species vertices that are not the one se-lected in step 20. While performing the reduc-tion, some entries of other rows may be affected.The equivalent to this in Algorithm 1 is the up-date of the coefficients in some edges in step 27or the elimination of edges in step 29 or the ad-dition of edges in step 32. Then, obtaining theRREF of the matrix is equivalent to a combina-tion of the following items: • species vertices without edges, which areequivalent to the zero rows and/or • reaction(species) vertices connected to atleast one species(reaction) vertex, which areequivalent to the final columns(rows) of thematrix.Finally, testing unconditional binomiality in thematrix (via checking that it has at most onenonzero entry at each row) is equivalent to check-ing that the components of the final graph con-tain either: • only one species vertex or, • one species vertex and one reaction vertexconnected to one another. (Proof of termination) For the proof of ter-mination, note that the graph has finitely manyedges and vertices. Hence, each of the loops ter-minates at some point. The loops at lines 6 and 22 eventually terminate as they iterate throughreaction vertices which are finite and there is nocreation of such vertices anywhere. Likewise, theloops at lines 4,10 and 17 terminate as they iter-ate through species vertices which are also finite,and no new vertices are created. (Proof of the complexity bound)
Thiscomes from the fact that each operation in thegraph theoretical algorithm corresponds to anoperation in the linear algebra approach.As a sidenote, for any transformed graph onecan generate a binomial coefficient matrix bytaking the species vertices as rows, reactions ver-tices as columns and the labels of edges as entriesand vice versa, from any binomial coefficient ma-trix one can generate a graph by applying thereverse procedure.
Example 3.
Graphs generated at each step ofAlgorithm 1 for the following RCRN has beenshown in Figure 2.3 B 2 C + A 2 D + 2 B 3 B.The final graph has a component that containstwo species vertices and three reaction verticesand therefore the RCRN is not unconditionallybinomial. On the other hand, Algorithm 1 in[26] gives the following matrix which shows thatthe RCRN is not unconditionally binomial sincethe first and second rows contain more than onenonzero entries.
We have implemented our graph theoretical ap-proach in Algorithm 1 in Section 3 as well asthe linear algebra approach in [26, Algorithm 1]in Maple. Both algorithms are available in theGit repository [22]. The performance of the im-plementations is tested on the biomodels fromthe BioModels repository [5] and the results arecompared to the experiments done on the samebiomodels using Gr¨obner basis and quantifierelimination in [16]. Note that in order to be ableto test unconditional binomiality of biomodels,5igure 2: Algorithm 1 on Example 3.we have assumed reversibility with free reverserate constants, while in [16] binomiality with pre-assigned values of rate constants are tested. Theresult of the experiments can be found in the Gitrepository [22].One can notice that overall, Algorithm 1 per-formed equal or better than [26, Algorithm 1].For a vast amount of models the difference ofthe performance of those two algorithms is al-most zero, which was predictable as the proof ofTheorem 5 suggests that the steps in the graphtheoretical algorithm and the linear algebra ap-proach are equivalent. However, for some mod-els (e.g., 205, 293 and 574) the graph theoreti-cal approach is faster. It is worth noting thatthese models have large binomial coefficient ma-trices, which may suggest that for RCRN witha large number of species and/or reactions thegraph theoretical approach can be faster. Thiswill be investigated further in future as the num-ber of models with large matrices is not enoughfor a through comparison at this stage.Comparing Algorithm 1 with the Gr¨obner ba- sis and quantifier elimination computations in[16] shows that for the great majority of thecases, the graph theoretical approach is muchfaster. More importantly, there are twentybiomodels that Gr¨obner basis and/or quantifierelimination methods in [16] run out of time (fora six hour limit of computations), while thosecases were handled in less than three seconds viaboth graph theory and linear algebra approaches.For six biomodels Gr¨obner basis computationsterminate in less than six hours, however, ourgraph theory algorithm as well as the linear alge-bra approach are at least 1000 times faster thanthose. For ten biomodels the graph theory (andthe linear algebra) approach is at least 500 timesfaster than the computations via Gr¨obner basisand quantifier elimination.An interesting observation is the quite highnumber (sixty nine) of the biomodels that arenot considered in [16] for the unclear numericvalues of their rate constants, whereas our graphtheoretical (and linear algebra) approaches on al-most all of those cases terminated in less than asecond.A comparison of the performance of our graphtheory algorithm with the algorithms proposedin [24, 6] are similar to the performance of thelinear algebra algorithm in [26] vs algorithms in[24, 6]. This is because, as it is mentioned ear-lier in this section, the graph theoretical algo-rithm has a similar (or better) performance tothe linear algebra algorithm in [26]. In particu-lar, for two reversible biomodels in the database(Biomodels 491 and 492), the graph theoreticalmethod is more than twice faster than the lin-ear algebra in [26], which means that it is muchfaster than the algorithms in [24, 6]. More de-tails of some comparisons between the linear al-gebra methods in the literature for testing un-conditional binomiality and conditional binomi-ality can be found in [26, Section 3].
The present work introduces an efficient graphtheoretical algorithm for testing unconditional6inomiality in an RCRN, which is different fromthe conditional binomiality considered by sev-eral other authors. The algorithm is essentiallyequivalent to the linear algebra approach pre-sented earlier for testing unconditional binomi-ality. Implementations of the graph theoreticalalgorithm as well the linear algebra approach aredone in Maple and experiments are carried onover several biomodels.While the graph theoretical algorithm seem tohave a slight advantage over the linear algebraapproach, the experiments reveal drastic advan-tage of both of those methods over the existingalgorithms based on Gr¨obner basis and quanti-fier elimination. Additionally, many cases thatcould not be handle by the Gr¨obner basis andquantifier elimination methods in a reasonabletime, were tested in less than few seconds viathe graph theoretical approach.
Acknowledgments
This work has been supported by the bilateralproject ANR-17-CE40-0036/DFG-391322026SYMBIONT [1, 2].
References [1] F. Boulier, F. Fages, O. Radulescu,S. Samal, A. Schuppert, W. Seiler,T. Sturm, S. Walcher, and A. Weber. TheSYMBIONT project: Symbolic methods forbiological networks.
ACM Communicationsin Computer Algebra , 52(3):67–70, Decem-ber 2018.[2] F. Boulier, F. Fages, O. Radulescu,S. Samal, A. Schuppert, W. Seiler,T. Sturm, S. Walcher, and A. Weber. TheSYMBIONT project: Symbolic methodsfor biological networks.
F1000Research ,7(1341), 2018.[3] Bruno Buchberger.
Ein Algorithmus zumAuffinden der Basiselemente des Rest-klassenringes nach einem nulldimension-alen Polynomideal . Doctoral dissertation, Mathematical Institute, University of Inns-bruck, Innsbruck, Austria, 1965.[4] Bruno Buchberger. Ein algorithmisches kri-terium f¨ur die l¨osbarkeit eines algebraischengleichungssystems.
Aequationes Mathemat-icae , 4:374–383, 01 1970.[5] Vijayalakshmi Chelliah, Nick Juty, IshanAjmera, Raza Ali, Marine Dumousseau, Mi-hai Glont, Michael Hucka, Ga¨el Jalowicki,Sarah Keating, Vincent Knight-Schrijver,Audald Lloret-Villas, Kedar Nath Natara-jan, Jean-Baptiste Pettit, Nicolas Ro-driguez, Michael Schubert, Sarala M.Wimalaratne, Yangyang Zhao, HenningHermjakob, Nicolas Le Nov`ere, and CamilleLaibe. BioModels: Ten-year anniversary.
Nucl. Acids Res. , 2015.[6] Carsten Conradi and Thomas Kahle. De-tecting 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. Realquantifier elimination is doubly exponential.
J. Symb. Comput. , 5(1–2):29–35, February–April 1988.[9] Alicia Dickenstein, Mercedes P´erez Mill´an,Anne Shiu, and Xiaoxian Tang. Mul-tistationarity in structured reaction net-works.
Bulletin of Mathematical Biology ,81(5):1527–1581, Feb 2019.[10] David Eisenbud and Bernd Sturmfels. Bi-nomial ideals.
Duke Mathematical Journal ,84(1):1–45, Jul 1996.[11] Martin Feinberg. Complex balancing in gen-eral kinetic systems.
Arch. Ration. Mech.An. , 49(3):187–194, January 1972.[12] Martin Feinberg.
Foundations of Chemi-cal Reaction Network Theory , volume 2027f
Applied Mathematical Sciences . Springer,2019.[13] William Fulton.
Introduction to Toric Vari-eties , volume 131 of
Annals of MathematicsStudies . Princeton University Press, 1993.[14] Alexander N. Gorban and Vassili N.Kolokoltsov. Generalized mass action lawand thermodynamics of nonlinear Markovprocesses.
Math. Model. Nat. Phenom. ,10(5):16–46, 2015.[15] Dima Grigoriev. Complexity of decidingTarski algebra.
J. Symb. Comput. , 5(1–2):65–108, February–April 1988.[16] Dima Grigoriev, Alexandru Iosif, HamidRahkooy, Thomas Sturm, and Andreas We-ber. Efficiently and effectively recognizingtoricity of steady state varieties.
CoRR ,abs/1910.04100, 2019.[17] Dima Grigoriev and Andreas Weber. Com-plexity of solving systems with few indepen-dent monomials and applications to mass-action kinetics. In
Proceedings of the CASC2012 , volume 7442 of
LNCS , pages 143–154.Springer, 2012.[18] Alexandru Iosif and Hamid Rahkooy. Anal-ysis of the Conradi-Kahle algorithm for de-tecting binomiality on biological models.
CoRR , abs/1912.06896, 2019.[19] Alexandru Iosif and Hamid Rahkooy.Maplebinomials, a Maple package for test-ing binomiality, August 2019. Available athttps://doi.org/10.5281/zenodo.3564428.[20] Maplesoft, a division of Waterloo MapleInc., Waterloo, Ontario, 2019.[21] Ernst W. Mayr and Albert R. Meyer. Thecomplexity of the word problems for com-mutative semigroups and polynomial ideals.
Adv. Math. , 46(3):305–329, 1982.[22] Cristian Vargas Montero and HamidRahkooy. TestBinomialityRCRN, Maple modules for testing unconditional bi-nomiality of reversible chemical reac-tion networks, June 2020. Available athttps://doi.org/10.5281/zenodo.3901559.[23] Mercedes P´erez Mill´an and Alicia Dicken-stein. The structure of MESSI biologi-cal systems.
SIAM J. Appl. Dyn. Syst. ,17(2):1650–1682, 2018.[24] Mercedes P´erez Mill´an, Alicia Dickenstein,Anne Shiu, and Carsten Conradi. Chemicalreaction systems with toric steady states.
Bull. Math. Biol. , 74(5):1027–1065, October2012.[25] Hamid Rahkooy and Thomas Sturm.First-order tests for toricity.
CoRR ,abs/2002.03586, 2020.[26] Hamid Rahkooy and Thomas Sturm. Alinear algebra approach for detecting bi-nomiality of steady state ideals of re-versible chemical reaction networks.
CoRR ,abs/2002.12693, 2020.[27] AmirHosein Sadeghimanesh and ElisendaFeliu. The multistationarity structure ofnetworks with intermediates and a binomialcore network.
Bulletin of Mathematical Bi-ology , 81(7):2428–2462, May 2019.[28] Guy Shinar and Martin Feinberg. Concor-dant chemical reaction networks and thespecies-reaction graph.
Mathematical bio-sciences , 241, 08 2012.[29] Bernd Sturmfels.
Gr¨obner bases and con-vex polytopes , volume 8 of
University Lec-ture Series . AMS, Providence, RI, 1996.[30] Volker Weispfenning. The complexity of lin-ear problems in fields.
J. Symb. Comput. ,5(1–2):3–27, February–April 1988.[31] Douglas B. West.