Invariants and Other Structural Properties of Biochemical Models as a Constraint Satisfaction Problem
 Sylvain Soliman^{1}Email author
Affiliated with
DOI: 10.1186/17487188715
© Soliman; licensee BioMed Central Ltd. 2012
Received: 29 April 2011
Accepted: 29 May 2012
Published: 29 May 2012
Abstract
Background
We present a way to compute the minimal semipositive invariants of a Petri net representing a biological reaction system, as resolution of a Constraint Satisfaction Problem. The use of Petri nets to manipulate Systems Biology models and make available a variety of tools is quite old, and recently analyses based on invariant computation for biological models have become more and more frequent, for instance in the context of module decomposition.
Results
In our case, this analysis brings both qualitative and quantitative information on the models, in the form of conservation laws, consistency checking, etc. thanks to finite domain constraint programming. It is noticeable that some of the most recent optimizations of standard invariant computation techniques in Petri nets correspond to wellknown techniques in constraint solving, like symmetrybreaking. Moreover, we show that the simple and natural encoding proposed is not only efficient but also flexible enough to encompass sub/surinvariants, siphons/traps, etc., i.e., other Petri net structural properties that lead to supplementary insight on the dynamics of the biochemical system under study.
Conclusions
A simple implementation based on GNUProlog's finite domain solver, and including symmetry detection and breaking, was incorporated into the BIOCHAM modelling environment and in the independent tool Nicotine. Some illustrative examples and benchmarks are provided.
1 Background
1.1 Introduction
Reaction models like those of reactome.org, KEGG pathway database [1] or biomodels.net represent a growing part of Systems Biology especially for metabolic or signalling pathways, cellcycle and more generally postgenomic regulation systems. They build on established standards like BioPAX or SBML [2] to facilitate the exchange and comparison of models and benefit from a large number of available tools, especially ODE integration based simulators.
The use of Petri nets to represent those models, taking into account the difference between compounds and reactions in the graph, and make available various kinds of analyses is quite old [3], however it remains somehow focused towards mostly qualitative and structural properties. Some have been used for module decomposition, like (I/O) Tinvariants [4, 5], related to dynamical notions of elementary flux modes [6]. However, there is, to our knowledge, very little use of Pinvariant computation, which provides both qualitative information about some notion of module related to the "life cycle" of compounds, and quantitative information related to conservation laws  each Pinvariant defines a conserved moiety of the obtained ODE system, whatever the rate laws  and Jacobian matrix singularity  induced by any Pinvariant since it defines a linear dependency between variables. Conservation law extraction is actually already provided by a few tools, but then using numerical methods, based on the quantitative view of the model, and not integer arithmetic (as in direct Pinvariant analysis).
We present here a very simple way to incorporate invariant computation in an existing biological modelling tool, using constraint programming with symmetry detection and breaking. We compare it to other approaches and evaluate it, for the case of Pinvariants, on some examples of various sizes, like the MAPK cascade models of [7] and [8]. This experimentation is done through an implementation of the described method in the BIOCHAM modelling environment [9, 10], and in the independent tool Nicotine. We benchmark the efficiency against state of the art Petri net tools on various models. Finally we show that the presented approach allows to compute, within the same framework, other interesting structural properties like sub/surinvariants or siphons/traps, bringing even more insight into the dynamics of the biochemical system under study.
1.2 Petri net view of a reaction model
A Petri net is a bipartite oriented (weighted) graph of transitions, usually represented as square boxes, and places, usually represented as circles, that defines a (actually not unique) transition relation on markings of the net, i.e., multisets of tokens associated to places. The relation is defined by firings of transitions, i.e., when there are tokens (as many as the weight of the incoming arc) in all preplaces of a transition, they can be consumed and as many tokens as the weight on the outgoing arc are added to each postplace. The classical Petri net view of a reaction model is simply to associate biochemical species to places and biochemical reactions to transitions.
In this Petri net, starting from a marking with at least one token in A and in E, one can remove one of each to produce one token in AE (firing of t _{1} ) and then either remove it to add again one token to A and one to E (firing of t _{1} ), or to add one B and one E (firing of t _{2} ).
P (resp. T) invariants are defined, as usual, as vectors V representing a multiset of places (resp. of transitions) such that V · I = 0 (resp. I · V = 0) where I is the incidence matrix of the Petri net, i.e. I _{ ij } is the number of arcs from transition i to place j, minus the number of arcs from place j to transition i. Intuitively, a Pinvariant is a multiset representing a weighting of the places and such that any such weighted marking remains invariant by any firing; a Tinvariant represents a multiset of firings that will leave invariant any marking (see also section 2.1). As explained in the introduction, for reaction models these invariants are used for flux analysis, variable simplification through conservation law extraction, module decomposition, etc. Note that we are concerned with the classical invariant problem and thus restrict our study to integer weights. This is an important difference with respect to the aforementioned flux analyses but it arises from the fact that the biochemical models we studied did not come from metabolism but from the modelling of signal transduction pathways, cell cycle, circadian rhythm, etc. In all these cases the stoichiometry was integer and, for instance, the extracted conservation laws will include only integer number of molecules.
1.3 Related work
To compute the invariants of a Petri net, especially if this computation is combined with other Petri net analyses, like sinks and sources, traps, deadlocks, etc. the most natural solution is to use a Petri net dedicated tool like INA, PiNA, or Charlie for instance through the interface of Snoopy [11], which allows the import of SBML models as Petri nets. Standard integer methods like FourierMotzkin elimination will then provide an efficient means to compute P or Tinvariants (see for instance [12] for a review). These methods however generate lots of candidates which are afterwards eliminated and also need to incorporate some means (like equality class definition) to avoid combinatorial explosion at least in some simple cases, as explained in Section 2.2.
Another way to extract the minimal semipositive invariants of a model is to use one of the software tools that provide this computation for biological systems, generally as "conservation law" computation, and based on linear algebra methods like QR factorization [13]. This is the case for instance of the METATOOL [14] and COPASI [15] tools. The idea is to use a linear relaxation of the problem, which suits well very big graphs, but needs again a posteriori filtering of the candidate solutions. Moreover, these methods do not incorporate any means of symmetry elimination (see Section 2.2). A recent technique for elementary mode computation relies on Mixed Integer Programming (MIP) [16] and is thus quite similar in theory to the ideas of thus article, however it is tailormade for elementary modes whereas for invariants pure Integer Programming would be enough, it is focused around the computation of a partial basis of these modes, which is an important problem but not the focus in this article, and  once again  it does not incorporate any symmetry breaking.
Finally, the most recent developments in invariant computation rely on a symbolic encoding through Binary decision Diagrams [17]. The tools based on this technique can prove quite efficient and are not unrelated to the symbolic encoding we present here through constraints. However they do not seem to integrate symmetry detection, also rely on filtering for minimality and thus, though they provide a symbolic solution very fast in some cases, might also benefit from some of the ideas we present. See section 2.5 for a more precise evaluation.
2 Results and Discussion
2.1 Finding invariants as a Constraint Satisfaction Problem
We will illustrate our new method for computing the invariants with the case of Pinvariants (but Tinvariants, being dual, would work in the same fashion). Consider a Petri net with p places and t transitions, these transitions represent reactions L _{ i } → R _{ i }, where L _{ i } encodes the stoichiometry of the reactants as a vector over places, and R _{ i } the same for the products of the reaction. A Pinvariant is a vector s.t. V ^{ T } · I = 0, i.e. ∀1 ≤ i ≤ t V · L _{ i } = V · R _{ i }. Since those vectors all live in , it is quite natural to see this as a Constraint Satisfaction Problem (CSP) [18–20] with t (linear) equality constraints on p finite domain (FD) variables.
where obviously equation (2) is redundant.
The task is actually to find invariants with minimal support, with respect to set inclusion (a linear combination of invariants belonging to also being an invariant), i.e., having as few nonzero components as possible, these components being as small as possible, but of course non trivial, we thus add the constraint that V · 1 > 0.
Example 3 In our running example we thus add A + E + AE + B > 0.
Now, to ensure minimality the labelling is invoked from small to big values. This means that for each variable, if an enumeration remains necessary after constraint propagation, values are tried in an increasing order starting at 0. This is closely related to the enumeration strategy used in the mixed integer programming method of [16] that allows them to look for shortest elementary modes. Such a restriction in the construction of the basis might thus also be possible in our approach. Then, a branch and bound procedure is wrapped around this search for solutions, maintaining a partial base of Pinvariant vectors and adding the constraint that a new vector V is solution if , which means that its support is not bigger than that of any vector of the base.
Unfortunately, even with the last constraint, no search heuristic was found that makes removing subsumed Pinvariants unnecessary. Thus, if a new vector is added to , previously found vectors with a bigger support must be removed. Section 2.6 will demonstrate other structural properties for which this step is not necessary.
The algorithm can be summarized as follows:
Algorithm 1 Minimal invariants computation
1: post the CSP for invariant V: ∀1 ≤ i ≤ t V · L _{ i } = V · R _{ i } and V · 1 > 0
2: repeat
3: find a solution, enumerating from low to high
4: add the solution to the basis
5: remove nonminimal invariants from the basis if there are any
6: post the new constraint
7: until no solution found
8: expand symmetrical solutions of
This algorithm was implemented directly into Nicotine^{1} and then added to BIOCHAM [9], which are both programmed in GNUProlog, and allowed for immediate testing.
Example 4 In our running example we find two minimal semipositive Pinvariants:

E = AE = 1 and A = B = 0

A = B = AE = 1 and E = 0
2.2 Equality classes
The problem of finding minimal semipositive invariants is clearly EXPSPACEhard since there can be an exponential number of such invariants. For instance the model given in Example 5 (described in [12] among others, and called "classic XY" in [17], where × is the number of places between each pair of transitions and Y the number of transitions) has 2^{ n } minimal semipositive Pinvariants (each one with either A _{ i } or B _{ i } equal to 1 and the other equal to 0).
A first remark is that in this example, there is a variable symmetry between all the pairs (A _{ i }, B _{ i }) of variables corresponding to places. This symmetry is easy to detect (purely syntactical) and can be eliminated through the usual ordering of variables, by adding the constraints A _{ i } ≤ B _{ i }.
This classical CSP optimization is enough to avoid most of the trivial exponential blowups and corresponds to the initial phase of parallel places detection and merging of the equality classes optimization [21] for the standard FourierMotzkin algorithm. Note however that in that method, classes of equivalent variables are detected and eliminated before and during the invariant computation, which would correspond to local symmetry detection and was not implemented in our prototype.
Moreover, in [21], equality class elimination is done through replacement of the symmetric places by a representative place. The full method reportedly improves by a factor two the computation speed. Even if in the context of the original article this is done only for ordinary Petri nets (Petri nets where the weights are only 0 or 1), we can see that it can be even more efficient to use this replacement technique in our case:
Example 6
...
A + B ⇒ 4*C
...
Instead of simply adding A ≤ B to our constraints, which will lead to 3 solutions when C = 1 before symmetry expansion: (A, B) ∈ {(0, 4), (1, 3), (2, 2)}, replacing A and B by D will reduce to a single solution D = 4 before expansion of the subproblem A + B = D.
This partial detection of independent subproblems, which can be seen as a complex form of symmetry identification, can once again be done syntactically at the initial phase, and can be stated as follows: replace ∑_{ i } k _{ i } * A _{ i } by a single variable A if all the A _{ i } occur only in the context of this sum i.e., in our Petri net all pretransitions of A _{ i } are connected to A _{ i } with k _{ i } edges and to all other A _{ j } with k _{ j } edges and same for posttransitions. For a better constraint propagation, another intermediate variable can be introduced such that A = gcd(k _{ i }) · A'. In our experiments the simple case of parallel places (i.e., all k _{ i } equal to 1 in the sum) was however the one encountered most often.
2.3 Example, the MAPK Cascade
The MAPK signal transduction cascade is a well studied system that appears in lots of organisms and is very important for regulating cell division [22]. It is composed of layers, each one activating the next, and in detailed models shows two intertwined pathways conveying EGF and NGF signals to the nucleus.
A simple MAPK cascade model, that of [23] without scaffold, is used here as an example to show the results of Pinvariant computation.
Pinvariants of the MAPK cascade model of [23]
RAFK, RAFRAFK 
RAFPH, RAFPHRAF~{p1} 
RAF, MEKRAF~{p1}, RAFRAFK, RAFPHRAF~{p1}, MEK~{p1}RAF~{p1}, RAF~{p1} 
MEKPH, MEKPHMEK~{p1}, MEKPHMEK~{p1, p2} 
MEK, MAPKMEK~{p1, p2}, MEKRAF~{p1}, MEKPHMEK~{p1}, MEKPHMEK~{p1, p2}, MAPK~{p1}MEK~{p1, p2}, MEK~{p1}RAF~{p1}, MEK~{p1}, MEK~{p1, p2} 
MAPKPH, MAPKPHMAPK~{p1}, MAPKPHMAPK~{p1, p2} 
MAPK, MAPKMEK~{p1, p2}, MAPKPHMAPK~{p1}, MAPK~{p1, p2} MAPK~{p1}MEK~{p1, p2}, MAPK~{p1}, MAPKPHMAPK ~{p1, p2}, 
In the next section other examples are used as benchmarks of this method, they are all much bigger than this one, which had only about 30 compounds, however note that one of those is still a model of the MAPK signalling cascade.
Note that these 7 Pinvariants define 7 algebraic conservation rules (i.e., mass conservation) and thus decrease the size of the corresponding ODE model from 22 variables and equations to only 15.
2.4 Evaluation on other biochemical examples
Minimal semipositive Pinvariant computation on bigger models of biochemical reaction networks
All the curated models in the September 2010 release of biomodels.net were also tested and none of them required more than 1s to compute all its minimal Pinvariants.
We could not compare our results with those provided in [13] since the models they use, coming from metabolic pathways flux analyses, do not have an integer stoichiometry matrix, however the examples of Table 2 show the feasibility of Pinvariant computation by constraint programming for quite big networks. Note that for networks of this size, the upper bound of the domain of variables had to be set manually. It was actually set to the value 8, which is about the double of the maximum value in all the biological models we have encountered up to now. The only overapproximation of the upper bound found was the product of the l.c.m. of stoichiometric coefficients of each reaction, which explodes really fast and leads to unnecessarily long computation. The manual bound results in a loss of completeness, but it is not enforced either by QRfactorization methods, and does not seem to miss anything on real life examples.
Though they are not specifically suited for this task (i.e., finding integer invariants), we tried some of the most well known Elementary Flux Modes computing packages on these examples. METATOOL [14] and efmtool [25] were chosen, since both can be run as Matlab packages. The results are not included in Table 2 but are summarized with the nonbiochemical examples of next section.
2.5 Nonbiochemical benchmarks
Even if our main purpose is to use the insight on the dynamics gained from the structural properties computed by our CSP, an evaluation of the proposed method on nonbiochemical models remains of interest.
The literature on invariant computation is quite large, however there does not seem to exist any standardized benchmark. Each author selects some examples with different properties (see for instance [12] from which only a few examples are used in [17], even though it is cited as reference) and few reuse the previously published sets of examples.
Moreover, even when the software used in these articles is available, usually only binary implementations are available, and only for some specific architectures and through a specific request process. In some cases none is provided at all.
Therefore, using a machine comparable in specifications, we chose to reuse the data published in the most recent work, that of Ciardo et al. [17]. Since we had to reencode ourselves the selected examples, only a subset of their benchmarks is covered, namely the classical dining philosophers problem [26], the standard exponential invariant case [12] and the circular trains [27]. These seem to cover the whole range of different schemes appearing in [17].
Note that there are usually many symmetries in these parametric examples and thus that a more powerful (or manual) symmetry detection would be called for in these specific cases. Nevertheless, since (intracellular) biochemical systems usually do not generate such structure, we did not push further the integration of more advanced symmetry detection/breaking in our tools.
All the models used for the biochemical and nonbiochemical benchmarks can be found at: http://contraintes.inria.fr/~soliman/nicotine_data/
METATOOL's "CONSERVATION RELATIONS" were used when possible, but that only allows to find  as expected  91 out of the 10 billion invariants for the classic example, in 0.33s. Models were thus transposed such that METATOOL and efmtool's EFM search correspond to Pinvariant computation. Transposed models appear with a 'b' ending in the data repository. efmtool was given the SBML files as input whereas some .dat files were generated for METATOOL. For all the examples of this section as well as Kohn's map, METATOOL gave the error message "Cannot sort modes with more than 52 rows" that was interpreted as some kind of "out of memory" error. For efmtool, in the same cases (all examples of this section plus Kohn's map) the computation was stopped after 10 minutes or more, with messages like "iteration 43/116: 224850 modes, dt = 2040206ms." that were interpreted as overtime. Note however that as already stated, these packages do not focus on integer stoichiometric matrices and thus have a much broader scope that might explain their poor performance on our benchmarks.
Minimal semipositive Pinvariant computation on general (nonbiochemical) benchmarks of the literature
model  BDD V2  BDD V4  GreatSPN  Nicotine  Metatool CR  Metatool EFM  efmtool 

trains 1010  4.81  om  0.03  3.26  na (20)  om  ot 
classic 1010  0.01  0.01  ot  0.15  na (91)  om  ot 
philo 30  1.04  0.01  0.01  2.68  3.04  om  ot 
The CSP approach can therefore be seen as a kind of intermediate between purely implicit (i.e., solutions encoded, for instance as a decision diagram, and needing to be decoded to be displayed) and purely explicit methods. It also remains very flexible as next section will prove and could incorporate many more optimizations (variable ordering heuristics, more symmetry elimination, etc.) at a quite low cost.
All the 80 Petri nets of http://www.petriweb.org/ were also tested. Only one took more than 1s: model 1516, which took about 3s to compute 1133 minimal Pinvariants. Since we do not have data for the other approaches on these models they were not added to the table of results but they confirm the feasibility and generality of our approach.
We think that the structure of this kind of net is however very different (average degree, arc weights, etc.) from that of usual biochemical reaction models and intend to explore this distinction further in the future.
2.6 Generalizing the approach to other structural properties
An interesting feature of the presented method is that it is actually flexible enough to encompass other structural properties than place or transition invariants. This is, to our knowledge, not the case of other alternative methods.
If for the Petri net of Example 1 one obtained the constraints shown in Example 2 to compute Pinvariants, one can notice that they can easily be adapted to compute sur or subinvariants, i.e., weighted sums that can only grow (resp. decrease) during the evolution of the system (see [28], for instance, for a formal definition). Indeed the following CSP describes exactly all the subinvariants of the system and is obtained in the same manner but with ≤ instead of =.
Surinvariants would be obtained with ≥ instead of ≤.
Now, getting a basis of minimal sub/surinvariants can be done with the same branch and bound technique used for invariants, allowing to obtain information on pools of species of the biochemical system that, for instance, never increase during any ODE simulation.
One can go slightly farther and once again reuse the same machinery, including symmetry breaking, to compute siphons and traps of the Petri net (see [29] for definition and example of use in biology). This time a boolean CSP is obtained with the following constraints for the example of traps:
To compute siphons one simply need to reverse ⇒ into ⇐.
Note that in the boolean domain, the support minimality can be imposed by enumerating in increasing (lexicographic) order, there is no need for any a posteriori check of minimality (step 5 of Algorithm 1). The algorithm thus becomes:
Algorithm 2 Minimal traps computation
1: post the CSP for trap V
2: repeat
3: find a solution, enumerating from low to high
4: add the solution to the basis
5: post the new constraint
6: until no solution found
7: expand symmetrical solutions of
This computation of traps and siphons can actually bring information about the dynamics of the model, including temporal logic formulae that it satisfies^{2}, together with other structural properties [4, 30] they provide an interesting toolkit to analyze structurally the dynamics of a Systems Biology model.
3 Conclusion
Pinvariants of a biological reaction model are not so difficult to compute in most cases. They carry information about conservation laws that are useful for efficient and precise dynamical simulation of the system, and provide some notion of module, which is related to the life cycle of molecules. Tinvariants are already used more commonly, and get more and more focus recently.
We introduced a new method to efficiently compute P and Tinvariants of a reaction network, based on FD constraint programming. It includes symmetry detection and breaking and scales up well to the biggest reaction networks found. Completeness is lost on the biggest examples but we still look for a better upper bound on domains to restore it.
The idea of applying constraint based methods to classical problems of the Petri net community is not new, but seems currently mostly applied to the modelchecking. We argue that structural problems (invariants, sinks, attractors, etc.) can also benefit from the knowhow developed for finite domain CP solving, like symmetry breaking, search heuristics, flexibility, etc. and thus intend to generalize our approach to other problems of this category.
Endnotes
^{1} http://contraintes.inria.fr/~soliman/nicotine.html
^{2}This is the topic of a paper currently submitted to the CMSB 2011 conference. Depending on the outcome, a reference or a short explanation will be added.
Declarations
Acknowledgements
We thank the French ANR project BioTempo (ANR10BLANC0218) for its support.
Authors’ Affiliations
References
 Kanehisa M, Goto S: KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Research 2000, 28:27–30.PubMedView Article
 Hucka M, et al.: The Systems Biology Markup Language (SBML): A Medium for Representation and Exchange of Biochemical Network Models. Bioinformatics 2003,19(4):524–531.PubMedView Article
 Reddy VN, Mavrovouniotis ML, Liebman MN: Petri net representations in metabolic pathways. In Proceedings of the 1st International Conference on Intelligent Systems for Molecular Biology (ISMB). Edited by: Hunter L, Searls DB, Shavlik JW. AAAI Press; 1993:328–336.
 Gilbert D, Heiner M, Lehrack S: A Unifying Framework for Modelling and Analysing Biochemical Pathways Using Petri Nets. In CMSB'07: Proceedings of the fifth international conference on Computational Methods in Systems Biology, of Lecture Notes in Computer Science. Volume 4695. Edited by: Calder M, Gilmore S. Edinburgh, Scotland: SpringerVerlag; 2007.
 GrafahrendBelau E, Schreiber F, Heiner M, Sackmann A, Junker BH, Grunwald S, Speer A, Winder K, Koch I: Modularization of biochemical networks based on a classification of Petri net by Tinvariants. BMC Bioinformatics 2008.,9(90):
 Schuster S, Fell DA, Dandekar T: A general definition of metabolic pathways useful for systematic organization and analysis of complex metabolic networks. Nature Biotechnology 2002, 18:326–332.View Article
 Chickarmane V, Kholodenkob BN, Sauro HM: Oscillatory dynamics arising from competitive inhibition and multisite phosphorylation. Journal of Theoretical Biology 2007, 244:68–76.PubMedView Article
 Schoeberl B, EichlerJonsson C, Gilles E, Muller G: Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nature Biotechnology 2002,20(4):370–375.PubMedView Article
 Calzone L, Fages F, Soliman S: BIOCHAM: An Environment for Modeling Biological Systems and Formalizing Experimental Knowledge. Bioinformatics 2006,22(14):1805–1807.PubMedView Article
 Fages F, Soliman S, ChabrierRivier N: Modelling and Querying Interaction Networks in the Biochemical Abstract Machine BIOCHAM. Journal of Biological Physics and Chemistry 2004,4(2):64–73.View Article
 Rohr C, Marwan W, Heiner M: Snoopy  a unifying Petri net framework to investigate biomolecular networks. Bioinformatics 2010,26(7):974–975.PubMedView Article
 Colom JM, Silva M: Convex geometry and semiflows in P/T nets. A comparative study of algorithms for computation of minimal Psemiflows. In Advances in Petri Nets 1990, of Lecture Notes in Computer Science. Volume 483. Edited by: Rozenberg G. London, UK: SpringerVerlag; 1991:79–112.View Article
 Vallabhajosyulaa RR, Chickarmane V, Sauro HM: Conservation Analysis of Large Biochemical Networks. Bioinformatics 2005, 22:346–353.View Article
 von Kamp A, Schuster S: Metatool 5.0: fast and flexible elementary modes analysis. Bioinformatics 2006,22(15):1930–1931.PubMedView Article
 Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, Singhal M, Xu L, Mendes P, Kummer U: COPASI  a COmplex PAthway SImulator. Bioinformatics 2006,22(24):3067–3074.PubMedView Article
 de Figueiredo LF, Podhorski A, Rubio A, Kaleta C, Beasley JE, Schuster S, Planes FJ: Computing the shortest elementary flux modes in genomescale metabolic networks. Bioinformatics 2009,25(23):3158–3165.PubMedView Article
 Ciardo G, Mecham G, PaviotAdet E, Wan M: PSemiflow Computation with Decision Diagrams. In PETRI NETS '09: Proceedings of the 30th International Conference on Applications and Theory of Petri Nets. Berlin, Heidelberg: SpringerVerlag; 2009:143–162.
 Mackworth AK: Consistency in Networks of Relations. Artifcial Intelligence 1977, 8:99–118.View Article
 Meseguer P: Constraint satisfaction problems: an overview. AI Communications 1989, 2:3–17.
 Kumar V: Algorithms for ConstraintSatisfaction Problems: A Survey. AI Magazine 1992, 13:32–44.
 Law CF, Gwee BH, Chang J: Fast and memoryefficient invariant computation of ordinary Petri nets. IEE Proceedings: Computers and Digital Techniques 2007,1(5):612–624.View Article
 Roovers K, Assoian RK: Integrating the MAP kinase signal into the G1 phase cell cycle machinery. BioEssays 2000,22(9):818–826.PubMedView Article
 Levchenko A, Bruck J, Sternberg PW: Scaffold proteins may biphasically affect the levels of mitogenactivated protein kinase signaling and reduce its threshold properties. PNAS 2000,97(11):5818–5823.PubMedView Article
 Hardy S, Robillard PN: Visualization of the simulation data of biochemical network models: a painted Petri net approach. In Proceedings of the SCSC'07 summer computer simulation conference. Edited by: San Diego, CA. USA: Society for Computer Simulation International; 2007:802–808.
 Terzer M, Stelling J: Largescale computation of elementary flux modes with bit pattern trees. Bioinformatics 2008,24(19):2229–2235.PubMedView Article
 Graf S, Steffen B, Lüttgen G: Compositional Minimisation of Finite State Systems Using Interface Specifications. Journal of Formal Aspects of Computing 1996,8(5):607–616.View Article
 Genrich H: Predicate/transition nets. In Advances in Petri nets, of Lecture Notes in Computer Science. Volume 254. Edited by: Brauer W, Reisig W, Rozenberg G. SpringerVerlag; 1987:207–247.
 Sifakis J: Structural Properties of Petri Nets. In Proceedings of MFCS'78, 7th Symposium on the Mathematical Foundations of Computer Science, of Lecture Notes in Computer Science. Volume 64. Edited by: Winkowski J. SpringerVerlag; 1978:474–483.
 ZevedeiOancea I, Schuster S: Topological analysis of metabolic networks based on Petri net theory. In Silico Biology 2003.,3(29):
 Nabli F, Soliman S: Steadystate solution of biochemical systems, beyond SSystems via Tinvariants. In CMSB'10: Proceedings of the 8th International Conference on Computational Methods in Systems Biology. Edited by: Quaglia P. CoSBi, ACM; 2010:14–22.
 Calzone L, Gelay A, Zinovyev A, Radvanyi F, Barillot E: A Comprehensive Modular Map of Molecular Interactions in RB/E2F Pathway. Molecular Systems Biology 2008.,4(173):
 Kohn KW: Molecular Interaction Map of the Mammalian Cell Cycle Control and DNA Repair Systems. Molecular Biology of the Cell 1999,10(8):2703–2734.PubMed