Faten Nabli
email: faten7@gmail.com

Thierry Martinez

Franois Fages

Sylvain Soliman
email: sylvain.soliman@inria.fr

Francois Martinez

Fages

On Enumerating Minimal Siphons in Petri nets using CLP and SAT solvers: Theoretical and Practical Complexity

Keywords: Petri nets, Siphons, Traps, Systems Biology, SAT, Constraint Logic Programming

Petri nets are a simple formalism for modeling concurrent computation. They are also an interesting tool for modeling and analysing biochemical reaction systems, bridging the gap between purely qualitative and quantitative models. Biological networks can indeed be complex, large, and with many unknown kinetic parameters, which makes the development of quantitative models difficult. In this paper, we focus on the Petri net representation of biochemical reactions and on two structural properties of Petri nets, siphons and traps, that bring us information about the persistence of some molecular species, independently of the kinetics. We first study the theoretical time complexity of minimal siphon decision problems in general Petri nets, and present three new complexity results: first, we show that the existence of a siphon of a given cardinality is NP-complete; second, we prove that deciding the Siphon-Trap property is co-NP-complete; third, we prove that deciding the existence of a minimal siphon containing a given set of places, deciding the existence of a siphon of a given cardinality and deciding the Siphon-Trap property can be done in linear time in Petri nets of bounded tree-width. Then, we present a Boolean model of siphons and traps, and two methods for enumerating all minimal siphons and traps of a Petri net, by using a SAT solver and a Constraint Logic Program (CLP) respectively. On a benchmark of 345 Petri nets of hundreds of places and transitions, extracted from biological models from the BioModels repository, as well as on a benchmark composed of 80 Petri nets from the Petriweb database of industrial processes, we show that both the SAT and CLP methods are overall faster by one or two orders of magnitude compared to the state-of-the-art algorithm from the Petri net community, and are in fact able to solve all the enumeration problems of our practical benchmarks.

Introduction

Petri nets were introduced in the 60's as a simple formalism for describing and analyzing information processing systems that are characterized as being concurrent, asynchronous, non-deterministic and possibly distributed [START_REF] Peterson | Petri Net Theory and the Modeling of Systems[END_REF]. The use of Petri nets for representing biochemical reaction systems, by mapping molecular species to places and reactions to transitions, was introduced quite late in [START_REF] Reddy | Petri net representations in metabolic pathways[END_REF], together with some Petri net concepts and tools for the analysis of metabolic networks. In particular, the traditional Petri net concepts of place-invariants (P-invariants) and transition-invariants (T-invariants) have important interpretations in biochemical networks: P-invariants express structural conservation laws between molecular species, which correspond to linear invariants and possible variable eliminations in systems of ordinary differential equations, while T-invariants revealed to be equivalent to the notion of extremal fluxes in metabolic networks [START_REF] Zevedei-Oancea | Topological analysis of metabolic networks based on petri net theory[END_REF], one of the main tools for analyzing and optimizing metabolic networks [START_REF] Larhlimi | A new constraint-based description of the steady-state flux cone of metabolic networks[END_REF][START_REF] De Figueiredo | Can sugars be produced from fatty acids? a test case for pathway analysis tools[END_REF][START_REF] Kamp | Metatool 5.0: fast and flexible elementary modes analysis[END_REF]. Constraint logic programs have been proposed to compute P-invariants and T-invariants in [START_REF] Soliman | Invariants and other structural properties of biochemical models as a constraint satisfaction problem[END_REF] and [START_REF] Nabli | Steady-state solution of biochemical systems, beyond S-Systems via T-invariants[END_REF] respectively. Constraint programming methods have also been applied successfully to many other biology related problems. For instance by Devloo et al. to discover efficiently the steady-state of large gene regulation networks [START_REF] Devloo | Identification of all steady states in large biological systems by logical analysis[END_REF]. Fanchon and al. use constraints to infer ranges of parameter values from observations [START_REF] Fanchon | Modeling the molecular network controlling adhesion between human endothelial cells: Inference and simulation using constraint logic programming[END_REF] and for analysing discrete genetic regulatory networks [START_REF] Corblin | A declarative constraintbased method for analyzing discrete genetic regulatory networks[END_REF]. Bockmayr and Courtois [START_REF] Bockmayr | Using hybrid concurrent constraint programming to model dynamic biological systems[END_REF] use Hybrid Concurrent Constraint to model a variety of biological phenomena, such as reaching thresholds, kinetics, gene interaction or biological pathways. In [START_REF] Larhlimi | A new constraint-based description of the steady-state flux cone of metabolic networks[END_REF], Larhlimi and Bockmayr take advantage of the implicit representation that constraints bring, to describe the elementary flux cone of some metabolic pathways. Backofen et al. pioneered in [START_REF] Backofen | Application of constraint programming techniques for structure prediction of lattice proteins with extended alphabets[END_REF] the use of constraints and symmetry breaking for predicting the structure of proteins, etc.

In this paper, we consider the Petri net concepts of siphons and traps. A siphon is a set of places that, once unmarked, remains unmarked. A trap is a set of places that, once marked, can never loose all its tokens. Thus, siphons and traps have opposed effects on the token distribution in a Petri net. These structural properties provide sufficient conditions for reachability (whether the system can produce a given protein or reach a given state from a given initial state) and liveness (deadlock freedom from a given initial state) properties in ordinary Petri nets. It is proved that in order for a net to have all its transitions live, it is necessary that each siphon remains marked. Otherwise (i.e., once a siphon is empty), transitions having their input places in a siphon cannot be live. One way to keep each siphon marked is to have a marked trap inside it.

In fact, this condition is necessary and sufficient for a free-choice net to be live [START_REF] Peterson | Petri Net Theory and the Modeling of Systems[END_REF],

We first study the theoretical time complexity of siphon extraction problems in general Petri nets. It has been shown in [START_REF] Tanimoto | Finding minimal siphons in general petri nets[END_REF] and [START_REF] Yamauchi | Time complexity analysis of the minimal siphon extraction problem of petri nets[END_REF] that the problems of existence of a minimal siphon of a given cardinality (k-MinimalSiphon), or containing a given place (Q-MinimalSiphon), are NP-complete, and recently in [START_REF] Oanea | New algorithms for deciding the siphon-trap property[END_REF] that the siphon-trap property (stp) is in co-NP. Here we provide new theoretical complexity results. First we show the NP-completeness of the existence problem of a siphon of a given cardinality (k-Siphon) and thus the NP-hardness of MinimalCardinalitySiphon. Second, we prove that deciding the siphon-trap property is in fact co-NP-complete. Third, we prove that deciding the existence of a minimal siphon containing a given set of places, deciding the existence of a siphon of a given cardinality and deciding the Siphon-Trap property are of linear time complexity in Petri nets with bounded tree-width. These latter results follow from Courcelle's theorem.

Then we consider a simple Boolean model for defining siphons and traps, and two methods for enumerating the set of all minimal siphons and traps of a Petri net. The first method iterates the resolution of the Boolean model executed with a SAT solver, while the second proceeds by backtracking with a Constraint Logic Program over Booleans (CLP(B)).

We compare this Boolean constraint solving approach to the state-of-theart algorithms from the Petri net community described in [START_REF] Cordone | Enumeration algorithms for minimal siphons in petri nets based on place constraints[END_REF] for computing minimal sets of siphons and traps, which have already been shown to outperform Mixed Integer Linear Programs previously proposed in [START_REF] Murata | Petri nets: properties, analysis and applications[END_REF][START_REF] Cordone | Characterization of minimal and basis siphons with predicate logic and binary programming[END_REF]. On a benchmark composed of 345 curated biological models of hundreds of species and reactions each, from the BioModels1 repository [START_REF] Le Novère | BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems[END_REF], and of 80 Petri nets from the Petriweb2 [START_REF] Goud | Petriweb: A repository for petri nets[END_REF] database of industrial processes, we show that the SAT solver MiniSAT and CLP(B) solver GNU-Prolog are both faster by one or two orders of magnitude than the dedicated algorithms, and can in fact enumerate all solutions for all the intances of those benchmarks in a seconds. Finally, we question ourselves why these programs perform so well in practice, and provide some elements of explanation related to our theoretical complexity results.

Preliminaries on Petri Nets

Petri Nets

Definition 1 A Petri net graph N is a weighted bipartite directed graph N = (P, T, W ), where P is a non-empty finite set of vertices called places, T is a non-empty finite set of vertices called transitions, P ∩ T = ∅, and W : (P × T ) ∪ (T × P ) → N is a weight function attached to the arcs.

Intuitively, the weight of a P ×T arc represents the minimum number of tokens (molecules) necessary to enable a transition and the weight on a T × P arc represents the quantity produced. By abuse of notation, the weight is equal to one if it is omitted. The weight of a P × T (resp. T × P ) arc equals zero if there exist no arcs from place to transition (resp. from transition to place). Definition 2 A marking for a Petri net graph is a mapping m : P → N which assigns a number of tokens to each place. We say that a place p is marked if m(p) > 0, otherwise it is said to be unmarked. Definition 3 A Petri net is a 4-tuple (P, T, W, m 0 ) where (P, T, W ) is a Petri net graph and m 0 is an initial marking.

Let N = (P, T, W ) be a Petri net graph.

Definition 4

The set of predecessors (resp. successors) of a transition t ∈ T is the set of places

• t = {p ∈ P | W (p, t) > 0} (resp. t • = {p ∈ P | W (t, p) > 0}).
Similarly, the set of predecessors (resp. successors) of a place p ∈ P is the set of transitions

• p = {t ∈ T | W (t, p) > 0} (resp. p • = {t ∈ T | W (p, t) > 0}).
The set of predecessors (resp. successors) • S (resp. S • ) of a set of places S is the union of sets of predecessors (resp. successors) of each place p ∈ S:

• S = p∈S • p (resp. S • = p∈S p • ). The set of predecessors (resp. successors) • Q (resp. Q • ) of a set of transi- tions Q is the union of sets of predecessors (resp. successors) of each transition t ∈ Q: • Q = t∈Q • t (resp. Q • = t∈Q t • ).
Definition 5 N is ordinary if for all p ∈ P and for all t ∈ T , W (p, t) ≤ 1 and

W (t, p) ≤ 1. Definition 6 A transition t is enabled at marking m when ∀p ∈ • t : m(p) ≥ W (p, t).
For every two markings m, m : P → N and every transition t ∈ T , there is a transition step m t → m if for all p ∈ P , m(p) ≥ W (p, t) and m (p) = m(p) -W (p, t) + W (t, p). m t → m means that the transition t is enabled in m and its firing leads to m . An enabled transition may or may not fire if there are other transitions enabled.

Example 1 The classical Petri net view of a reaction model is to associate biochemical species to places and biochemical reactions to transitions. The wellknown system of Michaelis-Menten enzymatic reactions can be represented by the Petri net depicted in Figure 1. It consists of three reactions that take place in two discrete steps: the first involves the reversible formation of a complex (AE) between the enzyme (E) and its substrate (A), and the second step involves an irreversible transformation of the product (B) with release of the enzyme.

A -A transition t ∈ T is dead at marking m if it is not enabled in any marking m reachable from m. -A marking m is dead if there is no transition enabled in m.

+ E AE → B + E A B E AE t 1 t -1 t 2
-A Petri net is deadlock free (weakly live) if there is no reachable dead markings.

Siphons and Traps

Let N = (P, T, W ) be a Petri net graph.

Definition 10 A trap is a non-empty set of places P ⊆ P whose successors are also predecessors:

P • ⊆ • P .
A siphon is a non-empty set of places P ⊆ P whose predecessors are also successors:

• P ⊆ P • .
A siphon (resp. a trap) is proper if its predecessor set is strictly included in its successor set: • P P • (resp. P • • P ). A siphon (resp. a trap) is minimal if it does not contain any other siphon (resp. trap).

It is worth remarking that a siphon in N is a trap in the dual Petri net graph, obtained by reversing the direction of all arcs in N . Note also that since predecessors and successors of an union are the union of predecessors (resp. successors), the union of two siphons (resp. traps) is a siphon (resp. a trap).

The Petri net of Example 1 has two minimal siphons, {A, AE} and {E, AE}. {E, AE} is both a minimal siphon and a minimal trap since

• {E, AE} = {E, AE} • = {t 1 , t -1 , t 2 }
Although siphons and traps are stable under union, it is worth noting that minimal siphons do not form a generating set of all siphons. A siphon is called a basis siphon if it cannot be represented as a union of other siphons [START_REF] Murata | Petri nets: properties, analysis and applications[END_REF]. Obviously, a minimal siphon is also a basis siphon, however, not all basis siphons are minimal. For instance, in Example 1, the generating set of siphons is formed by {A, AE}, {E, AE}, {A, AE, B} and {E, AE, B}, but only {A, AE} and {E, AE} are minimal, the two others cannot be obtained by union of minimal siphons.

The following propositions show that traps and siphons provide a structural characterization of some particular dynamical properties on markings.

Proposition 1 [START_REF] Peterson | Petri Net Theory and the Modeling of Systems[END_REF] For every subset P ⊆ P of places, P is a trap if and only if for any marking m ∈ N P with m p ≥ 1 for some place p ∈ P , and any marking m ∈ N P such that m σ → m for some sequence σ of transitions, there exists a place p ∈ P such that m p ≥ 1.

Proposition 2 [START_REF] Peterson | Petri Net Theory and the Modeling of Systems[END_REF] For every subset P ⊆ P of places, P is a siphon if and only if for any marking m ∈ N P with m p = 0 for all p ∈ P , and any marking m ∈ N P such that m σ → m for some sequence σ of transitions, we have m p = 0 for all p ∈ P .

Application to Deadlock Detection

One reason to consider minimal siphons is that they provide a sufficient condition for the non-existence of deadlocks. A deadlock occurs in a marked Petri net if no transition is enabled. It has been shown indeed that in a deadlocked (and marked) Petri net, all unmarked places form a siphon [START_REF] Chu | Deadlock analysis of petri nets using siphons and mathematical programming[END_REF]. The siphonbased approach for deadlocks detection checks if the net contains a proper siphon that can become unmarked by some firing sequence. In parrticular the following stp property provides a sufficient condition for ordinary Petri nets to be deadlock free.

Definition 11 Given a Petri net (P, T, W, m 0 ), the siphon-trap property (STP) holds when every siphon contains a marked trap.

Theorem 1 ( [START_REF] Commoner | Deadlocks in petri nets[END_REF]) An ordinary Petri net in which the STP holds is deadlock free. Proof. We just have to show that STP is preserved by transition firing and that a dead marking does not satisfy STP. The preservation follows from the fact that a marked trap remains marked after the firing of a transition. If a marking m is dead, then the set S = {p ∈ P | m(p) = 0} is such that S • = T since every transition should have an empty predecessor. Then • S ⊆ S • , and moreover S is non-empty since T is non-empty. Therefore S is a siphon and does not contain a marked trap.

The relevance of siphons and traps for other liveness properties in systems biology are reported in [START_REF] Heiner | Petri nets for systems and synthetic biology[END_REF].

Application to Systems Biology

One example of the relevance of traps and siphons in biology was given in [START_REF] Zevedei-Oancea | Topological analysis of metabolic networks based on petri net theory[END_REF] for the analysis of the production and accumulation of starch in potato tubers during growth, while starch is consumed after the tubers are deposited after the harvest. This can be seen by a purely structural analysis of the Petri net graph of the metabolic network, since starch and several of its precursors form traps in the reaction net during growth, while starch and possible intermediates of degradation form siphons after the harvest. A simplified version of this Petri net is depicted in Figure 2, where G 1 stands for glucose-1-phosphate, G u is UDP-glucose, S is the starch, I stands for intermediary species and P 1 and P 2 represent external metabolites [45]. In this network, either the branch producing starch (t 3 and t 4 ) or the branch consuming it (t 5 and t 6 ) is operative, as it is shown in Figure 3 and Figure 4 respectively. This is realized by a switching mechanism in the gene regulatory network with synthesis of the corresponding enzymes. Two Petri nets are thus derived from this model: one Petri net where t 5 and t 6 are removed (in this Petri net, t 3 and t 4 are said to be operative) and one Petri net where t 3 and t 4 are removed.

It can be easily observed that the set {G u , S} is a trap when t 3 and t 4 are operative: once a token arrives in S, no transition can be fired and the token remains there independently of the evolution of the system. Dually, {S, I} is a siphon when t 5 and t 6 are operative: once the last token is consumed from S and I, no transition can generate a new token in these places, so they remain empty.

Another interesting example of use of the concept of siphons and traps, also from [START_REF] Zevedei-Oancea | Topological analysis of metabolic networks based on petri net theory[END_REF], deals with the analysis of the role of the triosephosphate isomerase (TPI) in Trypanosoma brucei metabolism. Earlier, Helfert et al. [START_REF] Helfert | Roles of triosephosphate isomerase and aerobic metabolism in trypanosoma brucei[END_REF] supposed that glycolysis could proceed without TPI, but unexpected evidence that all system fluxes (Pyruvate, Glycerol) decrease was found which lead the authors to build a kinetic model for explaining that phenomenon. However, a purely structural explanation for the necessary presence of TPI in glycolysis and glycerol production was provided in [START_REF] Zevedei-Oancea | Topological analysis of metabolic networks based on petri net theory[END_REF] simply by showing the existence of siphons and traps in the model.

G 1 P 1 P 2 Gu S I t 1 t 2 t 3 t 4 t 5 t 6
Fig. 2: Petri net graph modeling the growth metabolism of the potato plant [START_REF] Zevedei-Oancea | Topological analysis of metabolic networks based on petri net theory[END_REF].

G 1 P 1 P 2 Gu S t 1 t 2 t 3 t 4
Fig. 3: Petri net graph modeling the growth metabolism of the potato plant [START_REF] Zevedei-Oancea | Topological analysis of metabolic networks based on petri net theory[END_REF] with the branch producing starch operative.

3 Theoretical Complexity of Siphon and Trap Properties

Preliminaries on SAT, FO and MSO

We shall prove some NP-completeness results using the following characterization of Boolean satisfiability. For a finite set of variables V = {x 1 , . . . x m }, let ¬V = {¬x 1 , . . . , ¬x m } denotes the set of negative literals constructed upon V . For a Boolean formula in conjunctive normal form

φ = c 1 ∧ • • • ∧ c n over V , we have for all 1 ≤ i ≤ n, c i = i,1 ∨ • • • ∨ i,ni
, and for all 1

≤ j ≤ n i , i,j ∈ V ∪ ¬V . Let us write C φ = {i ∈ N | 1 ≤ i ≤ n} and L φ = {(i, j) ∈ N 2 | 1 ≤ i ≤ n, 1 ≤ j ≤ n i }.
Lemma 1 A Boolean formula φ in conjunctive normal form is satisfiable if and only if there exists a subset X ⊆ L φ such that for all 1 ≤ i ≤ n, there exists 1 ≤ j ≤ n i such that (i, j) ∈ X,

G 1 P 1 P 2 S I t 1 t 2 t 5 t 6
Fig. 4: Petri net graph modeling the growth metabolism of the potato plant [START_REF] Zevedei-Oancea | Topological analysis of metabolic networks based on petri net theory[END_REF] with the branch consuming starch operative.

for all (i, j), (i , j ) ∈ X, i ,j = ¬ i,j .

Proof.

If µ : V → {0, 1} satisfies φ, we pose for all x ∈ V , µ(¬x) = 1 -µ(x)
and then it suffices to observe that X = {(i, j) ∈ L | µ( i,j ) = 1} satisfies the conditions of the lemma. Conversely, given a subset X ⊆ L φ satisfying these conditions, we pose µ : V → {0, 1} such that for all x ∈ V , µ(x) = 1 if there exists (i, j) ∈ X such that i,j = x and 0 otherwise. Then, we observe that µ satisfies φ.

We say that φ is satisfied by X if X is a subset of L φ satisfying the conditions of Lemma 1.

We shall also use the language of first-order logic (FO) to express properties over finite set and graph structures assumed to be fixed. For each finite set S, we assume a unary predicate S(x), also written x ∈ S by abuse of notation, which tests membership. Non-emptyness, S = ∅, is expressed by the FO formula ∃x(x ∈ S), and set inclusion, S ⊆ S , by ∀x(x ∈ S → x ∈ S ). Similarly for each Petri net graph, we assume two unary (set) predicates, place(x) and transition(x), which distinguish between places and transitions, and a binary predicate, edge(x, y), which tests incidence. This leads to the following characterization of siphons and traps: Definition 12 The siphon and trap properties can be defined in FO by the following schemas of formulae:

Siphon(S) : S = ∅ ∧ ∀p ∈ S place(p) ∧ ∀t(∃p ∈ S edge(t, p) → ∃p ∈ S edge(p, t)) Trap(S) : S = ∅ ∧ ∀p ∈ S place(p) ∧ ∀t(∃p ∈ S edge(p, t) → ∃p ∈ S edge(t, p))
We shall also provide a series of linear time complexity results by showing that some problems can be expressed in monadic second-order logic (MSO) over finite graph structures, in order to use Courcelle's theorem. MSO extends FO by the addition of second-order quantifiers over predicates, also noted ∃ and ∀, with the restriction to apply to unary predicates only, i.e. sets. An example of MSO formula is given in the proof of Theorem 4 in the next section.

Definition 13 ([42]

) A tree decomposition of a non-oriented graph G = (V, E) is a pair (X, T ) where T = (I, A) is a tree and X = (X i ) i∈I is a family of subsets of V such that

-i∈I X i = V , -for all {v, v } ∈ E, there exists i ∈ I such that {v, v } ⊆ X i , -for all i, j, k ∈ I, if j lies on the path from i to k in T , then X i ∩ X k ⊆ X j .
The tree-width of a tree decomposition is max i∈I |X i |-1. The tree-width tw(G) of G is the minimum tree-width taken over all possible tree decompositions of G.

The tree-width of an ordinary Petri net graph N = (P, T, W ) follows that definition given for graphs, by reading the Petri net as the non-oriented graph

G = (P ∪ T, E) where E = {{p, t} ⊆ P ∪ T | W (p, t) > 0 or W (t, p) > 0}.
Example 2 Let us consider the Petri net graph depicted in Figure 5. The nonoriented graph associate to this Petri net is depicted in Figure 6. Two possible tree decompositions of this non-oriented graph are shown in Figure 7. In both decompositions, each graph edge connects two vertices which belong to the same tree node. Graph vertices are adjacent only when the corresponding subtrees intersect. In the first tree decomposition (Figure 7 on the left), each tree node contains at most four vertices, hence the width of this decomposition is three. In the second tree (Figure 7 on the right), each tree node contains at most three vertices, hence the width of this decomposition is two, which is the optimal tree-width over all possible tree decompositions. Hence, the tree-width of the Petri net graph depicted in Figure 5 equals two.

A B C D E F G t 1 t 2 t 3 t 4 t 5
Fig. 5: Example of a Petri net with its associated non-oriented graph depicted in Figure 6 and two tree decompositions given in Figure 7.

Courcelle's theorem states that every graph property definable in MSO can be decided in linear time on graphs of bounded tree-width. Theorem 2 (Courcelle [START_REF] Courcelle | The monadic second-order logic of graphs i. recognizable sets of finite graphs[END_REF]) For a given formula P in monadic secondorder logic on the structure of graphs and for a given positive integer k, there exists a linear time algorithm that given a finite graph G of tree-width at most k decides whether G satisfies P .

Given a tree-automaton constructed from the formula, and a tree constructed from the graph decomposition that minimizes the tree-width, deciding if the automaton accepts the tree can be done in linear time. Courcelle's theorem shows the existence of such automata for MSO properties. However, the automata can be of hyper-exponential sizes. This makes the linear algorithm usually unusable in practice. There is some recent work in [START_REF] Courcelle | Automata for the verification of monadic second-order graph properties[END_REF][START_REF] Kneis | Courcelle's theorem -a game-theoretic approach[END_REF] on flyautomata and game-theoretic methods, aiming at making this result practical for certain properties of graphs with bounded clique-width. The applicability of these new approaches to our problem is however beyond the scope of this paper.

The Q-MinimalSiphon Problem

Definition 14 The problem FindMinimal is the following problem: "given a Petri net graph N = (P, T, W ) and a subset of places Q ⊆ P , find a minimal siphon S in N such that Q ⊆ S if there is any, or fail." Theorem 3 ([46]) The problem FindMinimal is polynomial.

Definition 15 Q-MinimalSiphon is the following decision problem: "given a Petri net graph N = (P, T, W ) and a subset of places Q ⊆ P , does there exist a minimal siphon S in N such that Q ⊆ S? ". Theorem 4 Q-MinimalSiphon is decidable in linear time for any class of Petri net graphs with bounded tree-width.

Proof. By Definition 12, siphons and traps are expressible in FO, and thus in MSO. We just have to verify that Q-MinimalSiphon can be expressed in MSO as follows:

Q-MinimalSiphon : ∃S(Siphon(S) ∧ Q ⊆ S ∧ ∀S (Siphon(S ) ∧ S ⊆ S → S ⊆ S ))
The linear time complexity then follows from Courcelle's theorem.

In the general case, Q-MinimalSiphon has been shown NP-complete in [START_REF] Yamauchi | Time complexity analysis of the minimal siphon extraction problem of petri nets[END_REF]. We find it useful here to provide a simpler proof of this result, by showing the NP-hardness of the following equivalent problem.

Definition 16 For a given Petri net graph N = (P, T, W ) and Q ⊆ P , a Q-hitting siphon is a siphon S ⊆ P of N such that for every siphon S ⊆ S of N , Q ⊆ S .

Definition 17 Q-HittingSiphon is the following decision problem: "given a Petri net graph N = (P, T, W ) and a subset of places Q ⊆ P , does there exist a Q-hitting siphon in N ? ".

Proposition 3 Given a Petri net graph N = (P, T, W ) and a subset of places Q ⊆ P , there exists a minimal siphon containing Q in N if and only if there exists a Q-hitting siphon in N . Proof. If S is a minimal siphon containing Q, then S is a Q-hitting siphon. Conversely, if S is a Q-hitting siphon, then there exists a minimal siphon S ⊆ S and since S is Q-hitting, then Q ⊆ S .

Theorem 5 Q-MinimalSiphon is NP-complete. Proof. Q-MinimalSiphon is in NP since FindMinimal is polynomial. We just have to show that Q-HittingSiphon is NP-hard.
Let us assume a sat instance described by a set V of variables and a Boolean formula φ over V in conjunctive normal form. Let N = (P, T, W ) be the ordinary Petri net graph, depicted in Figure 8, where

-P = L φ ∪ {q} and T = V ∪ ¬V ∪ C φ where q / ∈ L φ ∪ T , -W satisfies • (i, j) = {¬ i,j } (i, j) • = {i, i,j } for all (i, j) ∈ L φ • q = C φ q • = V ∪ ¬V q x ¬x y ¬y z ¬z x (1, 1)
¬y

(1, 2) z (1, n 1 ) 1 (2, 1) ¬x (2, 2) (2, n 2 ) 2 . . . . . . . . . . . . (n, 1) (n, 2) ¬x (n, nn)
n Fig. 8: Petri net used in the reduction of sat to Q-HittingSiphon in Theorem 5. In this Petri net, a siphon S contains the place q on top of the figure if and only if S contains at least one literal occurrence per clause. Moreover, S is {q}-hitting if and only if S does not contain a pair of contradictory literals.

Lemma 2 φ is satisfiable if and only if there exists a {q}-hitting siphon in N , i.e., a subset S ⊆ P such that S = ∅ and • S ⊆ S • and for all S ⊆ S, if S = ∅ and • S ⊆ S • , then q ∈ S . Proof. Suppose that φ is satisfied by X ⊆ L φ . We verify that X ∪ {q} is a {q}-hitting siphon. Indeed, • X ⊆ V ∪ ¬V = q • and, by hypotheses on X,

• q = C φ ⊆ X • , therefore • (X ∪ {q}) ⊆ (X ∪ {q}) • .
For any S ⊆ S such that S = ∅ and • S ⊆ S • , suppose that q / ∈ S , then since S = ∅, there exists (i, j) ∈ S and • (i, j) = {¬ i,j } ⊆ S • , therefore there exists (i , j ) ∈ S such that ¬ i,j ∈ (i , j ) • , therefore i ,j = ¬ i,j , which contradicts that S ⊆ X; thus, q ∈ S . That proves that X ∪ {q} is a {q}-hitting siphon.

Conversely, suppose that S is a {q}-hitting siphon. We verify that φ is satisfied by X = S ∩ L φ . Indeed, for all 1 ≤ i ≤ n, i ∈ • q, then i ∈ X • , therefore there exists 1 ≤ j ≤ n i such that (i, j) ∈ X. Suppose that there exist (i, j), (i , j ) ∈ X such that i ,j = ¬ i,j , then S = {(i, j), (i , j )} would be such that S = ∅ and • S ⊆ S • , despite q / ∈ S . Therefore, for all (i, j), (i , j ) ∈ X, i ,j = ¬ i,j . Thus, φ is satisfiable.

The theorem is then an immediate consequence of the lemma.

The k-Siphon Problem

Definition 18 The problem k-Siphon is the following decision problem: "given a Petri net graph N = (P, T, W ) and a positive integer k, does there exist a siphon S in N of cardinality k? ".

The linear time complexity immediately follows from Courcelle's theorem.

Theorem 6 k-Siphon is decidable in linear time with respect to the size of the Petri net for any class of Petri net graphs with a bounded tree-width. Proof. Given a Petri net graph N = (P, T, W ) and a positive integer k, there exists a siphon S in N of cardinality k if and only if ∃S(Siphon(S)∧card k (S)) is satisfied, where the formula card k (S)

card k (S) : ∃x 1 . . . x k   1≤i<j≤k x i = x j ∧ ∀x   x ∈ S → 1≤i≤k x = x i     checks that the cardinality of S is k.
We prove the NP-completeness of k-Siphon by polynomial reduction from the set covering problem, one of Karp's original NP-complete problems [START_REF] Karp | Reducibility among combinatorial problems[END_REF]. Let us recall that the problem k-SetCovering is the following decision problem: "given a finite set U (the universe), a subset S of P(U) and an integer k, does there exist a subset S ⊆ S of cardinality k such that U = S? ". Theorem 7 k-Siphon is NP-complete. Proof. k-Siphon is in NP since checking that a given set of places is a siphon of cardinality k is polynomial. NP-hardness comes by polynomial reduction from k-SetCovering: given a finite set U (the universe), a subset S of P(U) and an integer k, let N = (P, T, W ) be the ordinary Petri net graph such that P = S, T = U and for all t ∈ T , t • = P and • t = {S ∈ S | t ∈ S}. Then for every subset S ⊆ P , S is a siphon if and only if U = S. This shows that the optimization problem MinimalCardinalitySiphon is NP-hard.

The stp Problem

Definition 19 The Siphon-Trap Property stp is the following decision problem: "given a marked Petri net N = (P, T, W, m 0 ), does every siphon in N include a trap that contains a marked place? ". Theorem 8 stp is decidable in linear time for any class of Petri nets with bounded tree-width.

Proof. Here again, it suffices to remark that stp is expressible in MSO, with a unary predicate marked(p) for distinguishing marked places, and the following MSO formula:

stp : ∀S(Siphon(S) → ∃T (Trap(T ) ∧ T ⊆ S ∧ ∃p(marked(p) ∧ p ∈ T )))
The linear time complexity immediately follows from Courcelle's theorem.

In the general case, stp has been shown to be in co-NP in [START_REF] Oanea | New algorithms for deciding the siphon-trap property[END_REF], by reducing ¬stp to sat. Indeed, ¬stp expresses the existence of a siphon S such that every trap included in S does not intersect M . The encoding in sat focuses on the maximal trap included in S (the union of all the traps included in S), which is computed by removing iteratively places in S that cannot belong to a trap.

In fact, we can prove Theorem 9 stp is co-NP-complete.

Proof. Since stp is in co-NP [START_REF] Oanea | New algorithms for deciding the siphon-trap property[END_REF], it suffices to show that ¬stp is NP-hard.

Let us assume a sat instance described with a set V of variables and a Boolean formula φ over V in conjunctive normal form. Let N = (P, T, W, m 0 ) be the ordinary Petri net, depicted in Figure 9, where

-P = L φ ∪ {0} × V and T = V ∪ ¬V ∪ C φ ∪ {t} where t / ∈ P ∪ V ∪ ¬V ∪ C φ , -W satisfies • (i, j) = C φ ∪ { i,j } (i, j) • = {i, t} for all (i, j) ∈ L φ • (0, x) = {t} (0, x) • = {x, ¬x} for all x ∈ V -m 0 = 1 {0}×V .
Note that the set {0} × V is introduced in places to ensure that P ∩ T = ∅.

Lemma 3 φ is satisfiable if and only if (N, {0} × V ) satisfies ¬stp, i.e., there exists a subset S ⊆ P such that S = ∅ and • S ⊆ S • and for all T ⊆ S, if

T = ∅ and T • ⊆ • T , then T ∩ {0} × V = ∅.
Proof. Suppose that φ is satisfied by a set X ⊆ L φ . We verify that X ∪ {0} × V is a siphon and that it does not contain any trap intersecting {0} × V . Indeed, by hypotheses on X, C φ ⊆ X • , and

• (X ∪ {0} × V ) ⊆ C φ ∪ V ∪ ¬V ∪ {t} ⊆ (X ∪ {0} × V ) • . For any T ⊆ X ∪ {0} × V such that T = ∅ and T • ⊆ • T ,
suppose that there exists x ∈ T ∩ {0} × V , then {x, ¬x} ⊆ T • , therefore {x, ¬x} ⊆ • T and there exist (i, j), (i , j ) ∈ T ∩ X such that i ,j = ¬ i,j , (0, x) 2) S contains all the places (0, x) such that the variable x occurs in these literal occurrences; moreover, a marked place (0, x) belongs to a trap included in S if and only if S contains the pair of contradictory literal occurrences x and ¬x.

(0, y) (0,
that would contradict the hypotheses on X. Therefore (N, {0} × V ) satisfies ¬stp.

Conversely, suppose that there exists a subset S ⊆ P such that S = ∅ and • S ⊆ S • and for all T ⊆ S, if T = ∅ and T • ⊆ • T , then T ∩ {0} × V = ∅. We verify that φ is satisfied by X = S∩L φ . We have X = ∅ because otherwise there would exist (0, x) ∈ S ∩ {0} × V , then t ∈ • S, therefore t ∈ S • and there would exist (i, j) ∈ S ∩ L φ , which contradicts X = ∅. Thus, C φ ⊆ • S ⊆ S • and for all 1 ≤ i ≤ n, i ∈ S • so there exists 1 ≤ j ≤ n i such that (i, j) ∈ X. Suppose that there exist (i, j), (i , j ) ∈ X such that i ,j = ¬ i,j . Then, there exists (0, x) ∈ S ∩ {0} × V such that (0, x) • = { i,j , i ,j }. Then T = {(0, x), (i, j), (i , j )} will be such that T

• = { i,j , i ,j , i, i , t} ⊆ C φ ∪ { i,j , i ,j , t} = • T , but T ∩ {0} × V = ∅. Therefore, for all (i, j), (i , j ) ∈ X, i ,j = ¬ i,j .
The theorem is now an immediate consequence of the lemma.

Boolean Model for Minimal Siphons

The definition of siphons in FO given in Definition 12 directly leads to the following Boolean Constraint Satisfaction Problem (CSP):

Definition 20 Given a Petri net graph N = (P, T, W ) , the CSP Siphon(N ) is the triple (V, D, C) where -V = P , i.e. a variable is introduced for each place, -D(p) = {0, 1} for all p ∈ V , i.e. the variables are Boolean, -p∈S p = 1, i.e. a siphon is not empty,

-C = {(p = 1 → p ∈ • t p = 1) | p ∈ P, t ∈ • p}.
Proposition 4 For every Petri net graph N = (P, T, W ), the CSP Siphon(N ) is satisfied by a valuation ν if and only if {p ∈ P | ν(p) = 1} is a siphon. Proof. It suffices to check that for every non-empty set of places S, we have ∀p, ∀t ∈

• p, p ∈ S → p ∈ • t p ∈ S if and only if • S ⊆ S • .
The encoding in a SAT solver is short and direct. For each transition t in the set of predecessors of a place p, a clause C is added to the satisfiability problem. C is formed by a negated boolean variable associated to p or-ed with boolean variables in the set of predecessors of t. To avoid the trivial case of the empty siphon, one clause is added.

The set of all minimal siphons (w.r.t. set inclusion) can be enumerated in the set inclusion order, by restarting search each time a siphon S is found, with the additional constraint p∈S p = 0, for disallowing any superset of that siphon.

For enumerating siphons in set inclusion order, we compared two techniques: one by labeling an auxiliary cardinality variable in increasing order ( [START_REF] Nabli | Finding minimal siphons as a CSP[END_REF]), one by labeling directly the Boolean variables with increasing value selection (to test first the absence, then the presence of a place in the candidate solution). The second technique has revealed to be much more efficient. The correctness of this technique comes from the fact that if a solution S is found after a solution S, then the two paths in the search tree leading to these solutions have a least common ancestor node: this node corresponds to the labeling of a place p. By construction, p belongs to S and not to S, thus S ⊆ S.

We shall also consider a variant of the above CSP where the constraints p = 1 → p ∈ • t p = 1 are decomposed by introducing an intermediary Boolean variable for each transition. For every Petri net graph N = (P, T, W ), the CSP Siphon'(N ) = (V, D, C) is defined as follows.

-V = P ∪T , i.e. one variable is introduced for each place and each transition,

-D(x) = {0, 1} for all x ∈ V , -p∈S p = 1, i.e. a siphon is not empty, -C = {(p = 1 → t = 1) | p ∈ P, t ∈ • p} ∪ {(t = 0 → p = 0) | t ∈ T, p ∈ • t}.
It is immediate that if Siphon'(N ) is satisfied by a valuation ν, then Siphon(N ) is satisfied by ν| P , and conversely that if Siphon(N ) is satisfied by ν, then ν can be extended to a valuation satisfying Siphon'(N ).

This variant Siphon'(N ) enjoys the following property.

Proposition 5 For every Petri net graph N = (P, T, W ), the CSP Siphon'(N ) has the same tree-width as N .

Proof. It suffices to notice that the primal graph of Siphon'(N ) is isomorphic to the bipartite graph induced by N .

Implementations with SAT and CLP Solvers

This section describes two implementations of the above model with appropriate search strategy, one using an iterated SAT procedure and the other based on Constraint Logic Programming with Boolean constraints. The source code of our implementation is available on the web3 .

Iterative SAT Algorithm

The Boolean model can be directly interpreted using a SAT solver to check the existence of a siphon or trap. However, the enumeration of all minimal siphons requires a particular strategy for iterating the use of the SAT solver with new constraints ensuring the enumeration of the set of all minimal siphons. We use sat4j 4 , an efficient library of SAT solvers in Java for both satisfaction and optimization. It includes an implementation of the MiniSAT algorithm and implements the value selection strategy mentioned above: for each variable, the value 0 is tried before the value 1, for the enumeration of all solutions.

Example 3 Let us consider the Petri net depicting the enzymatic reaction of Example 1. In the first iteration, the problem amounts to solve the following encoding of Horn-dual clauses:

(¬A ∨ AE) ∧ (¬AE ∨ E ∨ A) ∧ (¬E ∨ AE) ∧ (¬E ∨ AE) ∧ (¬B ∨ AE)
The problem is satisfied with the values: E = B = 0 and A = AE = 1, which means that {A, AE} is a minimal siphon. In the second iteration, the clause ¬A ∨ ¬AE is added to ensure minimality, and the problem is satisfied with A = B = 0 and E = AE = 1 meaning that {E, AE} is also a minimal siphon. A new clause is added stating that either E or AE does not belong to the siphon and no more variable assignment can satisfy the problem.

Therefore, this model contains two minimal siphons: {A, AE} and {E, AE}. The enzyme E is a catalyst protein for the transformation of the substrate E in a product B. Such a catalyst increases the rate of the reaction but is conserved in the reaction.

Backtrack Replay CLP(B) Algorithm

The search for siphons can also be implemented with a Constraint Logic Program with Boolean constraints (CLP(B)). We use GNU-Prolog5 [START_REF] Diaz | Design and implementation of the GNU Prolog system[END_REF] for its efficient low-level implementation of Boolean constraint propagators, indexicals and search heuristics.

The enumeration strategy for enumerating all minimal siphons is indeed a variation of the branch-and-bound algorithm, where the search is restarted to find a non-superset siphon each time a new siphon is found. We have tried two versions of the branch-and-bound algorithm: with restart from scratch and by backtracking. In the branch-and-bound with restart method, it is essential to choose a variable selection strategy which ensures diversity. Indeed, an enumeration method with a fixed variable order accumulates failures by always trying to enumerate the same sets first and these failures are only lately pruned by the non-superset constraints. As a consequence, the developed search tree gets more and more dense after each iteration since the previous forbidden sets are repeatedly tried again. It is worth nboting that this phenomenon does not exist in SAT solvers thanks to no-good recording. In CLP, random variable order selection strategy provides a good diversity and performs much better than any uniform heuristics.

However, branch-and-bound by backtracking gives better performance when care is taken for posting the non-superset constraint only once, since reposting it at each backtrack step proved to be inefficient. We have implemented a backtrack replay strategy, i.e. a customized branch and bound procedure where the search is performed as follows:

1. each time a siphon is found, the path leading to this solution is memorized, 2. then the search is fully backtracked to the root in order to add to the model the new non-superset constraint, 3. and then the memorized path is rolled back and replayed to continue the search at the point it was stopped.

Rolling back the memorized path is conceptually equivalent to adding a lexicographic constraint that imposes to be greater than the last siphon found. However the performance of this constraint is not as good as an explicit rollback, because of lack of propagation.

The backtracking procedure can use any variable selection strategy, regardless of the diversity, since it ensures by construction that the already explored Fig. 10: Search tree developed with the backtrack replay strategy for enumerating the 64 minimal siphons of the model 239 of BioModels which contains 51 species and 72 reactions. The branches ending with a cross are solutions. It is remarkable that very few backtracks are necessary to enumerate all solutions using the backtrack replay strategy. sets will not be explored again. In our experiments, the most-constrained variable selection strategy is the one that provided the best performance and is therefore the one that was selected.

Figure 10 (generated with CLPGUI6 [START_REF] Fages | CLPGUI: a generic graphical user interface for constraint logic programming[END_REF]) depicts the search tree that is developed for enumerating the 64 minimal siphons of a biological model of 51 species and 72 reactions. Each sub-tree immediately connected to the root corresponds to the replay of the path with a minimality constraint added. The small number of backtrack points shows the remarkable efficiency of the backtrack replay strategy combined with a simple Boolean constraint propagator.

Evaluation

In the literature, many algorithms have been proposed to compute minimal siphons and traps of Petri nets. Since a siphon in a Petri net N is a trap of the dual net N , it is enough to focus on siphons, the traps are obtained by duality. Some algorithms are based on linear programming [START_REF] Murata | Petri nets: properties, analysis and applications[END_REF][START_REF] Cordone | Characterization of minimal and basis siphons with predicate logic and binary programming[END_REF], Horn clause satisfaction [START_REF] Kinuyama | Generating siphons and traps by petri net representation of logic equations[END_REF][START_REF] Minoux | Deadlocks and traps in petri nets as horn-satisfiability solutions and some related polynomially solvable problems[END_REF] or algebraic approaches [START_REF] Lautenbach | Linear algebraic calculation of deadlocks and traps[END_REF]. More recent state-of-the-art methods are presented in [START_REF] Cordone | Some results on the computation of minimal siphons in petri nets[END_REF][START_REF] Cordone | Enumeration algorithms for minimal siphons in petri nets based on place constraints[END_REF] and show the better performance of the dedicated algorithm of [START_REF] Cordone | Enumeration algorithms for minimal siphons in petri nets based on place constraints[END_REF]. This algorithm uses a recursive problem partitioning procedure to reduce the original search problem to multiple simpler search sub-problems. Each sub-problem has specific additional place constraints with respect to the original problem. This algorithm can be applied to enumerate minimal siphons, place-minimal siphons, or even siphons that are minimal with respect to a given subset of places. In this section we compare our results to that dedicated algorithm.

BioModels Benchmark

The BioModels [START_REF] Le Novère | BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems[END_REF] database is a repository of peer-reviewed, published, computational models of biological processes. These models are written in the Systems Biology Markup Language (SBML) exchange format and are mainly composed of reaction rules, with or without kinetics, and events. This resource allows biologists to store, search and retrieve models referenced in publications.

We consider the curated part of the repository BioModels version February 2013. Among these 404 curated models, 59 models do not hold any reaction (but events only), we thus consider the 345 reaction models from which a Petri net graph can be extracted from the structure of the reactions.

In addition, we consider the following complex biochemical models:

a model of E2F/Rb signaling from [START_REF] Calzone | A comprehensive modular map of molecular interactions in RB/E2F pathway[END_REF] which contains 408 molecular species and 534 reactions, and will be shown to contain 74 minimal siphons, -Kohn's map of the mammalian cell cycle control [START_REF] Kohn | Molecular interaction map of the mammalian cell cycle control and DNA repair systems[END_REF][START_REF] Chabrier-Rivier | Modeling and querying biochemical interaction networks[END_REF], which contains 509 species and 775 reactions, and will be shown to contain 80 minimal siphons.

Petriweb Benchmark

The database Petriweb7 [START_REF] Goud | Petriweb: A repository for petri nets[END_REF] is an attempt to make available a significant body of Petri nets examples in a public repository. They are written in the Petri Net Markup Language (PNML), an emerging standard format supported by many tools PetriWeb supports a restricted form of PNML, including flat, uncoloured nets, plus limited support for hierarchy The repository can be browsed with a web browser, and individual nets can be retrieved and uploaded in PNML.

The repository contains 80 Petri nets given with some interesting properties. The properties are defined by the repository administrator. The properties are intended to be checked automatically by external analysis tools.

Computational Results and Comparisons

In this section we compare the SAT method, using the MiniSAT solver implementation included in the SAT4J library, the CLP method, using GNU-Prolog, both described in Section 5, and the state-of-the-art dedicated algorithm of [START_REF] Cordone | Enumeration algorithms for minimal siphons in petri nets based on place constraints[END_REF].

Table 1 presents the CPU times in milliseconds for enumerating the sets of all minimal siphons of the Petri nets in our benchmark in Petriweb and BioModels (except Model 175 as explained below). The CPU times have been obtained on a PC with an intel Core processor 2.20 GHz and 8 GB of memory. For each benchmark, we provide the total number of models, the minimal, maximal and average numbers of siphons, and the total computation time in milliseconds for enumerating all of them. Quite surprisingly, on all these practical instances, both MiniSAT and GNU-Prolog solve the minimal siphon enumeration problem, in less than one millisecond in average, with a slightly better average performance for the CLP(B) program over the SAT solver. Furthermore, MiniSAT and GNU-Prolog outperform the dedicated algorithm by one or two orders of magnitude.

However, one particular model, number 175 in BioModels, has very high computation time and was excluded from Table 1. Table 2 presents the performance figures obtained on this model and on three other hardest instances for which we also provide the number of places and transitions. Even if the model is quite large, e.g. for Kohn's map of the cell cycle control with 509 species and 775 reactions, the computation time for enumerating its 81 minimal siphons is astonishingly short: one millisecond only. On these hard instances, the SAT solver is faster than the CLP(B) program by one or two orders of magnitude, and is the only algorithm to solve the problem for model 175, in 137 seconds. Model 175 represents a quantitative model that relates the EGF and HRG stimulations of the ErbB receptors to the activation of ERK and AKt in MCF-7 breast cancer cells [START_REF] Birtwistle | Ligand-dependent responses of the ErbB signaling network: experimental and modeling analysis[END_REF]. This is the first model to take into account both the ERK and AKt pathways, four ErbB receptors, and their simultaneous activation by two ligands. Previous models of ErbB (e.g. the model developed in [START_REF] Schoeberl | Computational modeling of the dynamics of the map kinase cascade activated by surface and internalized egf receptors[END_REF]) were limited to a single ErbB receptor because of combinatorial complexity. As a result, the ErbB signaling network is highly connected and indeed the underlying Petri net contains the highest number of arcs, and of organizations as remarked in [START_REF] Dittrich | Chemical organisation theory[END_REF][START_REF] Kaleta | Using chemical organization theory for model checking[END_REF], of the BioModels repository.

Clause Density Analysis

The enumeration of the set of all minimal siphons is a problem of enumeration of all the solutions of an NP-complete problem, so the question is: why are the CSP-based algorithms for enumerating siphons so efficient on the existing benchmarks, even on large graphs from systems biology? One possible explanation could be obtained by considering the well-known phase transition phenomenon in 3-SAT. The probability that a random 3-SAT problem is satisfiable has been shown to undergo a sharp phase transition as the ratio of the number of clauses over the number of variables crosses the critical value of about 4.26 [START_REF] Mitchell | Hard and easy distributions of sat problems[END_REF][START_REF] Crawford | Experimental results on the crossover point in satisfiability problems[END_REF], going from satisfiability to unsatisfiability with probability one, when the number of variables grows to infinity. It is in this region of the density that the SAT instances are difficult to decide, while before and after that density the instances are usually easy.

This density of SAT instances associated to the enumeration of all minimal siphons, grows during enumeration since clauses are added for each solution found. We can thus check whether the initial density of the 3-SAT instances associated to BioModels instances are greater than the critical value of 4.26. If the initial density is above the critical value, it will remain above, and the instances are thus easy because there will be a small number of solutions. On the other hand, if the initial density is below the threshold value of 4.26, the computation time may be long because the threshold value were the 3-SAT hardest instances are may be traversed, and the clauses may be satisfiable with an exponential number of solutions.

The density of a SAT instance is: The initial density distribution of all BioModels instances are illustrated in the histogram of Figure 11.

This histogram shows that the initial density is in fact below the critical value for the majority of models. The initial 3-SAT density of our hardest model number 175 equals 2.39. Since the density grows during enumeration by adding the clauses for minimality, the possibility to traverse the critical region of density exists. Density considerations thus do not suffice to explain why we are so efficient in enumerating all the solutions of an NP-complete problem in large classes of practical Petri nets.

Tree-Width Analysis

There is also a rich literature about the polynomial-time complexity of CSPs when the constraint hypergraph is bounded relatively to a variety of graph measures, including cutwidth and tree-width [START_REF] Gottlob | A comparison of structural CSP decomposition methods[END_REF]. Since we have shown that both Q-MinimalSiphon and k-Siphon are expressible as CSPs, the existence of a polynomial algorithm for deciding these properties for Petri net graphs of bounded tree-width follows by Prop. 5. However, since these problems are also expressible in Monadic Second Order (MSO) logic, they are in fact of linear time complexity in this case, as shown by Theorems 4 and 6.

It is thus interesting to measure the tree-widths of the Petri nets of our benchmark. QuickBB8 is a program for computing the tree-width of a graph. When given enough time, this algorithm yields the exact tree-width of the We have applied QuickBB on the 432 curated models of our BioModels benchmark. For 31 models, the exact tree-width could not be computed in a time-out of one hour, but the tree-width was bounded by 23. For the remaining 342 models, the exact tree-width was computed and was always less than 10 as shown in Figure 12. The tree-width of our three hardest instances are given in Table 3, and have been determined to be less than 10 and 15. These treewidth values are relatively small values for graphs of hundreds of species and reactions. The BioModels benchmark thus seems to indicate that biochemical reaction networks have a bounded tree-width less than ten. This is in agreement with the idea that even very large biochemical processes are not fully intercon-nected as in a grid, but composed through interface molecular species. These considerations suggest that biochemical networks are of bounded tree-width, in which case Theorems 4, 6 and 8 show that the minimal siphon decision problems are indeed tractable, and in fact of linear time theoretical complexity.

Conclusion

Siphons and traps in Petri nets are meaningful pools of places that display a specific behaviour in the Petri net dynamics, and that guarantee some persistence properties in the simulation of a system of biochemical reactions, independently of the kinetics.

We have described a Boolean model for the problem of enumerating all minimal siphons in a Petri net and have compared two Boolean methods to the state-of-the-art algorithm from the Petri net community [START_REF] Cordone | Enumeration algorithms for minimal siphons in petri nets based on place constraints[END_REF]. On the benchmark of 345 biological models from the curated part of the BioModels repository, the Boolean method for enumerating all minimal siphons using MiniSAT is very efficient. It also scales very well in the size of the net. The CLP(B) program also solves all but one instances of the benchmark, with a better performance than MiniSAT in average, but does not scale-up as well on the largest size Petri nets, such as for instance on Kohn's map with 509 species and 775 reactions. The MiniSAT solver and the CLP(B) program both outperform the dedicated algorithms of the Petri net wommunity by one or two orders of magnitude and solve instances out of reach of these other algorithms.

The efficiency of the MiniSAT and CLP(B) methods for enumerating all solutions of an NP-complete problem for all, including large, instances of our practical benchmarks was quite surprising and lead us to study the theoretical complexity of these problems. Besides the proofs of NP-completeness of the existence of a siphon of a given cardinality, and of co-NP-completeness of the siphon-trap property, we have shown that these decision problems are tractable in Petri nets of bounded tree-width. Then we have shown that the BioModels benchmark of large biochemical networks have indeed a relatively small treewidths.

These various results militate for the analysis of biochemical networks with Petri net concepts and Constraint Programming tools.

Fig. 1 : 1 . 7 m 1 t1→

 1171 Fig. 1: Petri net associated to the example of Michaelis-Menten of Example 1, displayed here with a marking enabling transition t 1 .

Fig. 6 :BFig. 7 :

 67 Fig. 6: Non-oriented graph associated to the Petri net of Figure5

Fig. 11 :

 11 Fig. 11: Distribution of the initial density of the initial 3-SAT problems considered to enumerate all minimal siphons of the Petri net from the BioModels benchmark.

Fig. 12 :

 12 Fig. 12: Variation of the computed tree-widths as a function of size (places and transitions) of the Petri nets associated to BioModels. The computed exact tree-widths are depicted by red circles. When not known, the computed upper bounds are depicted by blue bars.

Table 1 :

 1 Computation time in milliseconds on the BioModels and Petriweb benchmarks.

	Benchmark	#	# siphons	siphons size	total time (ms)
		model min-max (avg.) min-max (avg.)	dedicated	MiniSAT	GNU
					algorithm		Prolog
	BioModels	345	0-64 (4.21)	1-413 (3.10)	19734	611	195
	Petriweb	80	0-11 (2.85)	0-7 (2.03)	2325	156	6

Table 2 :

 2 Computation time in milliseconds on the hardest instances of biochemical networks.

Table 3 :

 3 Tree-width of the hardest instances of BioModels database.

http://www.biomodels.net/

http://www.petriweb.org/

http://lifeware.inria.fr/wiki/software/

http://www.sat4j.org/

http://www.gprolog.org/

http://contraintes.inria.fr/∼fages/CLPGUI

http://www.petriweb.org/

http://graphmod.ics.uci.edu/group/quickbb/

Acknowledgements This work was supported by the French OSEO project Biointelligence.