Niyogi

Weinberger Dominique Smale

Attali

Christopher Fillmore

Ishika Ghosh

André Lieutier

Elizabeth Stephenson

Mathijs Wintraecken

Optimal homotopy reconstruction results à la

Keywords: Subject Classification Theory of computation → Computational geometry Homotopy, Inference, Sets of positive reach

In this article we show that the proof of the homotopy reconstruction result by Niyogi, Smale, and Weinberger can be streamlined considerably using Federer's work on the reach and several geometric observations. While Niyogi, Smale, and Weinberger restricted themselves to C 2 manifolds (M) with positive reach, our proof extends to sets (S) of positive reach.

The sample we consider does not have to lie directly on the set of positive reach. Instead, we assume that the two one-sided Hausdorff distances (δ and ε) -between the sample P to the set S, are bounded. We provide explicit bounds in terms of ε and δ, that guarantee that there exists a parameter r such that the union of balls p∈P B(p, r) deformation retracts to S. We provide even better bounds for the manifold case. In both cases, our bounds improve considerably on the state-of-the-art in almost all settings. In fact the bounds are optimal.

Introduction

Can we infer the topology of a subset of the Euclidean space if we are only given a partial geometric information about it? Under which conditions is such inference possible? These questions were first motivated by the shape reconstruction of 3-dimensional objects.

There, the partial geometric information was represented by a finite, in general noisy, set of points obtained from photogrammetric or lidar measurements [START_REF] Amenta | A simple algorithm for homeomorphic surface reconstruction[END_REF][START_REF] Berger | A survey of surface reconstruction from point clouds[END_REF][START_REF] Boissonnat | Geometric structures for three-dimensional shape representation[END_REF][START_REF] Boissonnat | Shape reconstruction from planar cross sections[END_REF][START_REF] Chazal | Smooth manifold reconstruction from noisy and nonuniform approximation with guarantees[END_REF].

More recently, the same questions have arisen in the context of learning and topological data analysis (TDA). In these fields, one seeks to recover a (relatively) low-dimensional support of a probability measure in a high-dimensional space, given a (finite) data set drawn from this probability measure [START_REF] Boissonnat | Geometric and Topological Inference[END_REF][START_REF] Chazal | Geometric inference for measures based on distance functions[END_REF][START_REF] Edelsbrunner | Computational topology: an introduction[END_REF]. Assuming the support is a manifold, one calls this process manifold learning [START_REF] Pless | A survey of manifold learning for images[END_REF].

In [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF], Niyogi, Smale, and Weinberger showed that, given a C 2 manifold of positive reach and a sufficiently dense point sample on (or near) the manifold, the union of balls of certain radii centered on the point sample captures the homotopy type of the manifold. By the nerve theorem, the homotopy type of the union of balls is shared by the Čech complex of the point sample [START_REF] Björner | Topological methods[END_REF][START_REF] Edelsbrunner | Triangulating topological spaces[END_REF]. From the Čech complex we can then infer the topological information such as the homology groups of the underlying manifold. Niyogi, Smale, and Weinberger's homotopy reconstruction result has led to generalizations including [START_REF] Attali | Vietoris-rips complexes also provide topologically correct reconstructions of sampled shapes[END_REF][START_REF] Chazal | A sampling theory for compact sets in Euclidean space[END_REF][START_REF] Kim | Homotopy Reconstruction via the Cech Complex and the Vietoris-Rips Complex[END_REF][START_REF] Wang | Topological inference of manifolds with boundary[END_REF].

Related work and contribution

Niyogi, Smale, and Weinberger [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF] seem to have been unaware of Federer's theory of sets of positive reach [START_REF] Federer | Curvature measures[END_REF]. In this article we discuss the extension of Niyogi, Smale, and Weinberger's result to sets of positive reach and use Federer's theory to streamline their proof (for both the manifold case and general sets of positive reach). We recall the relevant results from Federer's work [START_REF] Federer | Curvature measures[END_REF] in the appendix. Attali et al., Chazal et al., and Kim et al. studied even more general settings [START_REF] Attali | Vietoris-rips complexes also provide topologically correct reconstructions of sampled shapes[END_REF][START_REF] Chazal | A sampling theory for compact sets in Euclidean space[END_REF][START_REF] Kim | Homotopy Reconstruction via the Cech Complex and the Vietoris-Rips Complex[END_REF]. Their proofs are, however, different from ours, and (necessarily) more involved. While our proof is short and elementary, we provide optimal bounds on the distance between the manifold (or the set) of positive reach and its sample. The best previously known bounds (for some particular cases) can be found in [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF] and [START_REF] Attali | Vietoris-rips complexes also provide topologically correct reconstructions of sampled shapes[END_REF], respectively. Let M denote a manifold of positive reach, S a set of positive reach and let P be a sample. All sets are assumed to be compact unless stated otherwise. We denote the reach of a set X by rch(X ) and let R be a non-negative real number such that R≤rch(S) (resp. R≤rch(M)).

We denote the one-sided Hausdorff distance1 from P to M (resp. S) by δ, and the one-sided Hausdorff distance from M (resp. S) to P by ε. Distinguishing the two one-sided Hausdorff distances seems natural to the authors, because in measurements one would expect the measurement error δ (with the exception of some small number of outliers) to be often smaller than the sampling density ε. Similar assumptions seem to be common in the learning community, see e.g. [START_REF] Levin | The approximation power of moving least-squares[END_REF]. Niyogi, Smale, and Weinberger made similar assumptions on the support of the measure from which they sampled. We achieved the following conditions on ε and δ which, if satisfied, guarantee the existence of a radius r > 0 such that the union of balls p∈P B(p, r) deformation-retracts onto M (resp. S):

If S has positive reach, we achieved the condition ε

+ √ 2 δ ≤ ( √ 2 -1)R. If moreover S = M is a manifold, we achieved the condition (R -δ) 2 -ε 2 ≥ 4 √ 2 -5 R 2 (in this case we assume that δ ≤ ε ≤ R).
The set of pairs (ε, δ) that satisfy the conditions depicted in Figure 1. In this article we only consider samples for which we have precise bounds on ε and δ. In [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF], the authors also consider a setting where the point sample is drawn from a distribution centred on the manifold. They still recover the homotopy type of the underlying manifold with high probability. Our results can be applied to improve the bounds also in this context. However, we have not discussed this in detail, since combining both results is straightforward. We stress that in [START_REF] Boissonnat | Geometric and Topological Inference[END_REF][START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF][START_REF] Wang | Topological inference of manifolds with boundary[END_REF], the authors use ε/2 instead of our ε. As shown in Section 4, our conditions on ε and δ are optimal for sets of dimension at least 2 in the following sense: If the conditions are not satisfied we can construct a set of positive reach S (resp. manifold M) and a sample P such that there is no r ≥ 0 for which p∈P B(p, r) has the same homology as S (resp. M). In the case of 1-dimensional manifolds, we expect the optimal condition to be weaker. We would like to emphasize that, for noiseless samples, that is when δ = 0, both the constant √ 2 -1 (for general sets of positive reach), and the constant 2 -√ 2 (for manifolds) compare favorably with the previously best known constant 1 2 3 5 from [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF] for manifolds. In Proposition 7.1 of [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF] one encounters the condition ε < (3 -√ 8)R, for a particular case of the family of settings we consider, namely when the sampling condition is expressed through an upper bound ε on the Hausdorff distance, in other words δ = ε in our setting. The same constant 3 -√ 8 was proven independently in [2, Theorem 4] for general sets of positive reach. While the optimality was not previously proven, this constant is in fact optimal when δ = ε both for general sets with positive reach and for manifolds. To contrast the two related results in [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF], for δ = 0 and δ = ε respectively, with our bounds, we portray them as black dots in Figure 1, that depicts our sufficient and necessary conditions on the pair (ε, δ). The case of manifolds with boundary has been studied in [START_REF] Wang | Topological inference of manifolds with boundary[END_REF]Theorem 3.2], assuming lower bounds on both the reach of the manifold and the reach of its boundary. We improve also on this result by considering a manifold with boundary as a particular case of a set with positive reach, while our bounds only depend on the reach of the set itself and not the one of its boundary.

Outline, setting, and notation

In Section 2, we reprove -and extend -the homotopy reconstruction result of Niyogi, Smale, and Weinberger. In Section 3, we determine the bounds on the quality parameters ε and δ of the sample. Finally, in Section 4, we show that these bounds are optimal. We denote the closed ball in a Euclidean space centred at a point p with radius r by B(p, r). Open balls are denoted by B(p, r) • . If no centre is indicated we tacitly assume that the centre is at 0. The Minkowski sum of two sets A and B is denoted by A ⊕ B. The normal cone of a closed set S at a point q ∈ S is denoted by Nor(q, S) (see also Definition 21). The closest point projection on a set S is denoted by π S .

▶ Universal Assumption 1. Throughout this article, we work with a closed set S ⊆ R d with positive reach rch(S), and let R > 0 be a constant satisfying R ≤ rch(S). Furthermore, we consider a set P ⊆ R d , such that the one-sided Hausdorff distance from P to S is δ, and the one-sided Hausdorff distance from S to P by ε. That is,

S ⊆ P ⊕ B(ε) and P ⊆ S ⊕ B(δ).

We assume that δ, ε < R. If the set S is a submanifold of the Euclidean space, we denote it by M.

Figure 1 In blue-gray we depict the region in (ε, δ)-space for which there exists a radius r such that the union of balls P ⊕ B(r) captures the homotopy type of a set of positive reach. We do the same in yellow for a manifold. We stress that for δ ≥ ε the fact that a set of positive reach is also a manifold does not lead to better bounds. We have rescaled such that the reach is 1. The black points indicate the bounds that were known to Niyogi, Smale and Weinberger.

For most applications the assumption δ ≤ ε seems to be natural, but we do not need this for our analysis.

2

The geometric argument

In this section we show that if the union of balls P ⊕ B(r) covers a sufficiently large neighbourhood of S, and the parameter r is not too large, then P ⊕ B(r) deformation-retracts to S.

▶ Theorem 2. Assume that a parameter α > 0 is small enough, so that the α-neighbourhood S ⊕ B(α) of the set S is contained in the union of balls P ⊕ B(r). In other words,

S ⊕ B(α) ⊆ P ⊕ B(r). ( 1
)
If, moreover,

r 2 ≤ (R -δ) 2 -(R -α) 2 , ( 2
)
then, for any point q ∈ S, the intersection (q + Nor(q, S)) ∩ B(q, R) ∩ (P ⊕ B(r)) of the normal cone q + Nor(q, S), the ball B(q, R), and union of balls P ⊕ B(r), is star shaped, with the point q as its 'centre'. Furthermore, P ⊕ B(r) deformation-retracts onto S along the closest point projection.

▶ Remark 3. The statement of Theorem 2 does not use the hypothesis S ⊆ P ⊕ B(ε) from the universal assumption.

We refer to Figure 3 for a pictorial overview of the proof of Theorem 2.

(a) Any point in q + Nor(q, S) a distance less than α from S is covered by P ⊕ B(r).

(b) If (q +Nor(q, S))∩(P ⊕B(r)) is not star shaped there exists a point x where the segment L reenters a ball B(p ′ , r) (in blue) in P ⊕ B(r) after having left P ⊕ B(r) closer to q.

(c) The center p ′ of the ball B(p ′ , r) lies inside the half space H. The half space H lies at least a distance α from q.

(d) The ball of radius R is 'tangent' to the set S, thus it cannot contain any point of S in its interior. Since the distance between p ′ and S is bounded by δ, p ′ has to lie outside of the ball of radius R -δ. This contradicts the fact that p ′ lies in the halfspace H and is not too far from the normal space q + Nor(q, S). Proof of Theorem 2. We prove the claim by contradiction. For any point q ∈ S, the set (q + Nor(q, S)) ∩ (S ⊕ B(α)) is contained in the union of balls P ⊕ B(r), as illustrated in Figure 3a. Assume that there exists a point q ∈ S and a vector v ∈ Nor(q, S), with ∥v∥ = 1, such that the intersection of P ⊕ B(r) with the segment

L = def. {q + λv | λ ∈ [0, R)}
consists of several connected components (as illustrated in Figure 3b). Thanks to Equation ( 1), the connected component that contains q has length at least α. Let x be first point along L, seen from q, lying inside a connected component of (P ⊕ B(r)) ∩ L that does not contain q.

Then x lies at the intersection of the segment L and a ball B(p ′ , r), with p ′ ∈ P . Furthermore, the point p ′ is contained in the open half space H orthogonal to the vector v, that does not contain q, and whose boundary contains x. We stress that if p ′ lies on the boundary of H then the line L is tangent to the sphere ∂B(p ′ , r), which is still compatible with star-shapedness.

The situation is illustrated in Figure 3c.

Let z = q + Rv be the open endpoint of L. Since, by [14, Theorem 4.8 [START_REF] Edelsbrunner | Computational topology: an introduction[END_REF]], the intersection S ∩ B(z, R) • is empty and the distance between p ′ and S is bounded by δ, we know that

p ′ / ∈ B(z, R -δ) • . Thus, p ′ ∈ A = def. H ∩ (R d \ B(z, R -δ) • ).
Figure 4 The center of the ball creating a new connected component along the normal cone q + (q, S) is constrained to belong to the set A.

The sphere ∂B(z, R -δ) has a non-empty intersection with the plane ∂H. Indeed, the sphere passes through point q + δv which does not belong to H while its center z belongs to H; see Figure 3d. We can thus pick a point y in the intersection ∂H ∩ ∂B(z, R -δ). By the Pythagoras theorem, the minimal squared distance between A and L is:

inf a∈A ℓ∈L ∥a -ℓ∥ 2 = ∥x -y∥ 2 = ∥z -y∥ 2 -(∥z -q∥ -∥x -q∥) 2 ≥ (R -δ) 2 -(R -α) 2 ,
as illustrated in Figure 4. Hence, if

r 2 ≤ (R -δ) 2 -(R -α) 2 , ( 2
)
the ball B(p ′ , r) does not intersect L. Therefore, L ∩ (P ⊕ B(r)) cannot have more than one connected component. The set (q + Nor(q, S)) ∩ B(q, R) ∩ (P ⊕ B(r)) is thus star-shaped with centre q. Since r satisfies Equation ( 2), we deduce that δ + r < R, and thus

P ⊕ B(r) ⊆ S ⊕ B(R) • .
Thanks to this, the fact that the set (q + Nor(q, S)) ∩ B(q, R) ∩ (P ⊕ B(r)) is star-shaped with centre q, and [14, Theorem 4.8 [START_REF] Edelsbrunner | Computational topology: an introduction[END_REF]], the map

H :[0, 1] × (P ⊕ B(r)) → P ⊕ B(r), (t, x) → (1 -t)x + tπ S (x),
is well-defined. Furthermore, thanks to [14, Theorem 4.8 [START_REF] Boissonnat | Geometric and Topological Inference[END_REF]], the projection π S is continuous, and even Lipschitz. Thus, the map H is a deformation retract from the union of balls P ⊕ B(r) to the set S. ◀

In Appendix B, we provide an alternative proof of Theorem 2, similar to an argument presented in [START_REF] Chazal | Smooth manifold reconstruction from noisy and nonuniform approximation with guarantees[END_REF].

Bounds

In this section we first compute the bounds on the size α of the neighbourhood S ⊕ B(α) covered by the union of balls P ⊕ B(r) in terms of ε, δ, and r. We then combine these bounds with Equation (2) to infer (optimal) upper bounds on ε and δ, for which there exists a radius r such that the deformation retract from P ⊕ B(r) to S is possible. We do so first for sets of positive reach and then for manifolds. Somewhat counter-intuitively, it turns out to be easier to determine the bounds for sets of positive reach.

Sets of positive reach

For sets of positive reach, the bound on α is almost trivial. Nevertheless, it is tight, as we will see in Section 4.

▶ Lemma 4. Suppose that S ⊆ P ⊕ B(ε) for some ε ≥ 0. Then, for all α ≤ r -ε, the α-neigbourhood S ⊕ B(α) of S is contained in the finite union of balls P ⊕ B(r). That is,

S ⊕ B(α) ⊆ P ⊕ B(r). Proof. Indeed, S ⊕ B(α) ⊆ (P ⊕ B(ε)) ⊕ B(α) ⊆ P ⊕ B(ε + α) ⊆ P ⊕ B(r). ◀
From this, we derive the bounds on ε and δ (in terms of R).

▶ Proposition 5. If ε and δ satisfy ε + √ 2 δ ≤ ( √ 2 -1)R, (3)
there exists a radius r > 0 such that the union of balls P ⊕ B(r) deformation-retracts onto S along the closest point projection. In particular, r can be chosen as:

r ∈ 1 2 R + ε - √ ∆ , 1 2 R + ε + √ ∆ ,
where

∆ = 2(R -δ) 2 -(R + ε) 2 .
Proof. We combine the bound from Lemma 4 with the conditions of Theorem 2. More precisely, inserting α = r -ε in Equation (2) yields that

r 2 + (R -r + ε) 2 ≤ (R -δ) 2 . ( 4
)
Using the abc-formula for quadratic equations, this is equivalent to

r ∈ 1 2 R + ε - √ ∆ , 1 2 R + ε + √ ∆ ,
where

∆ = 2δ 2 + R 2 -4δR -2Rε -ε 2 = 2(R -δ) 2 -(R + ε) 2 is the discriminant. This interval is non-empty if the discriminant is non-negative, that is, if ε + √ 2 δ ≤ ( √ 2 -1)R. ◀ ▶ Remark 6.
The parameter δ is not necessarily smaller than ε, even if this would be natural in most applications.

Manifolds with positive reach

In this section, we show that the bounds from Proposition 5 can be improved further if the set of positive reach is a manifold. In Lemma 4, we used the triangle inequality to set α = r -ε.

If S is a manifold, however, the parameter α can be increased using more subtle arguments than the triangle inequality: Manifolds with positive reach are C 1,1 smooth 2 , i.e., differentiable with Lipschitz derivative. Moreover, Federer's normal cone Nor(q, M) (Definition 21) coincides with the 'classical' normal space N q M of an n-dimensional submanifold M of R d at every point q ∈ M. The tangent and normal cones of manifolds of positive reach are n-and (d -n)-dimensional linear spaces, respectively, that are not only dual, but also orthogonal.

In Lemma 7, we establish a lower bound for the parameter α in the case that S = M is a manifold. This bound is tight, as we will see in Section 4.

▶ Lemma 7. Suppose that M ⊆ P ⊕ B(ε) for some ε ≥ 0. Then, for any r ≥ α ≥ 0 satisfying

r 2 ≥ α 2 + α R R 2 + ε 2 -(R -δ) 2 + ε 2 , ( 5
) the α-neighbourhood M ⊕ B(α) of M is contained in the finite union of balls P ⊕ B(r). That is, M ⊕ B(α) ⊆ P ⊕ B(r).
Proof. Given a point q ∈ M, the tangent cone T q M and the normal cone N q M are orthogonal vector spaces satisfying T q M ⊕ N q M = R d . Since M ⊆ P ⊕ B(ε), the intersection P ∩ B(q, ε) is non-empty. Let p ∈ P ∩ B(q, ε). The vector p -q decomposes uniquely as 2 Topologically embedded manifolds with positive reach are C 1,1 embedded [START_REF] Federer | Curvature measures[END_REF][START_REF] Lieutier | Embedded topological manifolds have positive reach if and only if they are C 1,1[END_REF][START_REF] Lytchak | On the geometry of subsets of positive reach[END_REF][START_REF] Lytchak | Almost convex subsets[END_REF][START_REF] Rataj | Curvature measures of singular sets[END_REF][START_REF] Rataj | On the structure of sets with positive reach[END_REF]. with u T ∈ T q M, u N ∈ N q M, ∥u T ∥ = ∥u N ∥ = 1, and λ T , λ N ≥ 0 (see Figure 5). Since ∥p -q∥ ≤ ε,

p -q = λ T u T + λ N u N ,
λ 2 T + λ 2 N ≤ ε 2 . ( 6
)
Thanks to [14, Theorem 4.8 ( 12)] (Theorem 24 in Appendix A), the sets M and B(q + Ru N , R) • do not intersect, and thus:

p / ∈ B(q + Ru N , R -δ) • .
Hence, (p -q -Ru N ) 2 ≥ (R -δ) 2 . Applying the decomposition of p -q we obtain

(λ T u T + (λ N -R)u N ) 2 ≥ (R -δ) 2 ,
which implies that

λ 2 T + (R -λ N ) 2 ≥ (R -δ) 2 .
Combining this result with Equation ( 6) implies that

ε 2 -λ 2 N + (R -λ N ) 2 ≥ (R -δ) 2 ,
which can be rewritten as

2Rλ N ≤ R 2 + ε 2 -(R -δ) 2 . ( 7
)
Choose a vector v ∈ Nor(q, M) with ∥v∥ = 1, and let α ≥ 0. Then,

(p -(q + αv)) 2 = ((λ N u N -αv) + λ T u T ) 2 = (λ N u N -αv) 2 + (λ T u T ) 2 ≤ (λ N + α) 2 + λ 2 T ≤ (λ N + α) 2 + ε 2 -λ 2 N (by (6)) = α 2 + 2αλ N + ε 2 .
Using inequality [START_REF] Boissonnat | Shape reconstruction from planar cross sections[END_REF] to substitute 2λ N , we further obtain:

(p -(q + αv)) 2 ≤ α 2 + α R R 2 + ε 2 -(R -δ) 2 + ε 2 .
Thus, if

r 2 ≥ α 2 + α R R 2 + ε 2 -(R -δ) 2 + ε 2 ,
then the point q + αv lies in B(p, r) ⊆ P ⊕ B(r). Since this inclusion holds for any q ∈ M and v ∈ N q M, M ⊕ B(α) ⊆ P ⊕ B(r). ◀

As in Proposition 5, we now derive a bound on ε.

▶ Proposition 8. If ε and δ satisfy (R -δ) 2 -ε 2 ≥ 4 √ 2 -5 R 2 (8)
and δ ≤ ε, there exists a radius r > 0 such that the union of balls P ⊕ B(r) deformationretracts onto M along the closest point projection. In particular, r can be chosen as as in [START_REF] Chazal | Smooth manifold reconstruction from noisy and nonuniform approximation with guarantees[END_REF].

The bound is illustrated in Figure 1.

Proof. We combine the bound from Lemma 7 with the conditions of Theorem 2. More precisely, combining Equations ( 2) and ( 5) yields the following sufficient condition for L ∩ (P ⊕ B(r)) to be connected:

α 2 + α R R 2 + ε 2 -(R -δ) 2 + ε 2 ≤ r 2 ≤ (R -δ) 2 -(R -α) 2 . ( 9
)
The inequality between leftmost and rightmost members of [START_REF] Chazal | A sampling theory for compact sets in Euclidean space[END_REF], which needs to be satisfied for a non-empty range of values for r to exist, can be rearranged as:

0 ≥ ε 2 -(R -δ) 2 + R 2 + α 1 R ε 2 -R 2 -(R -δ) 2 + 2α 2 .
Using the abc-formula for quadratic equations, the above inequality is satisfied if α ∈ [α min , α max ], with

α min = 1 4 (R -δ) 2 + R 2 -ε 2 R - √ ∆ , α max = 1 4 (R -δ) 2 + R 2 -ε 2 R + √ ∆ , ( 10
)
where the discriminant ∆ is

∆ = 1 R 2 ε 2 -(R -δ) 2 2 -10 ε 2 -(R -δ) 2 -7R 2 .
The discriminant can be viewed as a polynomial in

y = ε 2 -(R -δ) 2 . Solving ∆(y) = 0 with respect to y yields y = R 2 5 ± 4 √ 2 . This in turn implies that ∆ is non-negative if either ε 2 -(R -δ) 2 ≤ R 2 5 -4 √ 2 or ε 2 -(R -δ) 2 ≥ R 2 5 + 4 √ 2 .
Thanks to Assumption 1, we are only interested in the case where ε, δ < R, and thus we can ignore the second inequality. Hence the interval [α min , α max ] is non-empty if

ε 2 -(R -δ) 2 ≤ 5 -4 √ 2 R 2 .
Plugging the bounds on α as given in [START_REF] Chazal | Geometric inference for measures based on distance functions[END_REF] into Equations ( 2) and ( 5) yields

1 + α min R ε 2 + α 2 min + α min R R 2 -(R -δ) 2 ≤ r 2 ≤ (R -δ) 2 -(R -α max ) 2 . ( 11
)
◀ ▶ Remark 9. We restricted ourselves to the case where δ ≤ ε, because if δ > ε, the fact that the set is a set of positive reach is a manifold [instead of an arbitrary set ]no longer helps. The geometric reason for this is that p in Figure 5 may lie in N q M.

Optimality

In this section we prove that the bounds provided in Section 3 are optimal in the following sense: To prove Propositions 10 and 11, we construct the set S, the manifold M, and the corresponding samples in Examples 15 and 16, respectively. Due to rescaling it suffices to construct sets of reach R = 1.

▶
▶ Remark 12. The dimension of the set in both examples will be two-dimensional. As mentioned in the introduction, we expect that better bounds can be obtained for onedimensional sets, i.e. curves, possibly with boundary.

▶ Remark 13. When δ ≥ ε, which in Figure 1 corresponds to the area above the diagonal δ = ε, the same bound is optimal whether the set is assumed manifold or not. Indeed, in this case the union of annuli S in Example 15 can be replaced by a union of circles, namely the inner boundaries of the annuli. Thus, the bound is tight for manifolds, including one-dimensional submanifolds in R 2 .

▶ Remark 14. For simplicity, the examples 15 and 16 are not connected. However, one could easily cut up our examples and glue the pieces together in a way that preserves the reach so that resulting examples still yield Propositions 10 and 11.

Sets of positive reach

The construction of the set proving Proposition 10 goes as follows:

▶ Example 15. We define S to be the disjoined union of annuli A i in R 2 , each of which has inner radius 1 and outer radius 1 + 2ε, and lies a distance at least 2 away from the other. We number the annuli from i = 1. A bound on their total number will be discussed later. The sample P consists of circles C i of radius 1 + ε lying in the middle of the annuli (C i ⊆ A i ), and a discrete set of points.

We describe the positioning of the discrete set of points in an inductive manner. Each pair of points, p i and pi , lies inside the annulus A i , at a distance δ from A i , and a distance 2r i-1 from each other. We define r 0 = δ+ε 2 . The distance r i is defined as follows: Let q i be the intersection point of the bisector of p i and pi , and the circle C i . Then r i is the circumradius of the triangle p i pi q i . The situation is illustrated in Figure 6. Denoting the angle at q i in the triangle p i pi q i by α i , we can express the inductive step as

r i+1 = r i sin α i ,
as a consequence of the law of sines. It is important to observe that, since r i+1 ≥ r i for each i, the sequence of angles satisfies α i+1 ≤ α i . The construction thus necessarily stops if either r i+1 = r i (i.e., α i = π/2) or r i+1 > 1 -δ. In the former case, p i+1 and pi+1 coincide with p i and pi , respectively. In the latter, it is not possible to place p i+1 and pi+1 in the annulus according to our conditions. The latter case is reached after finitely many iterations if r i+1 -r i is lower bounded, or equivalently if α i is bounded away from π/2. The proof of Proposition 10 focuses on establishing this bound.

Proof of Proposition 10. Assume that δ and ε fail to satisfy the bound (3), that is,

√ 2(1 -δ) -1 < ε. ( 12
)
We first show that each triangle p i pi q i is acute, that is, that α i < π/2. This is the case precisely when the height of the triangle, 1 + ε -(1 -δ) 2 -r 2 i , is larger than r i . This is a quadratic inequality in r i ,

r i < 1 + ε -(1 -δ) 2 -r 2 i ⇐⇒ 0 < (1 + ε) 2 -(1 -δ) 2 -2(1 + ε)r i + 2r 2 i .
The inequality holds for all

r i if 0 > 4 2(1 -δ) 2 -(1 + ε) 2 ,
which is equivalent to bound [START_REF] Edelsbrunner | Computational topology: an introduction[END_REF]. In particular, S consists of finitely many annuli. For any r ∈ [0, 1 -δ), the union of balls P ⊕ B(r) has different homology than the set S. Indeed, for r ∈ [0, r 0 ), each set (C i ∪ {p i , pi }) ⊕ B(r) has three connected components, for r ∈ [r i-1 , r i ), the set (C i ∪ {p i , pi }) ⊕ B(r) has the homotopy type of two circles that share a point, that is, a wedge of two circles, also known as a bouquet. The situation is illustrated in Figure 7. At the same time, each annulus A i ⊆ S has the homotopy type of a circle.

Manifolds

The construction of the manifold proving Proposition 11 goes as follows:

▶ Example 16. We define M to be a finite union of tori of revolution T i in R 3 . Each of these tori is the 1-offset of a circle of radius 2 in R 3 . Put differently, each T i is -up to Euclidean transformations -the surface of revolution around the y-axis of a circle of radius 1 in the xz-plane, centred at (2, 0, 0). In order for the reach of M to be equal to 1, we assume that the tori lie at a distance at least 2 apart from each other. The set T i is illustrated in Figure 8 (in blue). The sample P consists of tori with parts cut out, denoted by C i , and a discrete set of points. We construct each set C i by taking a δ-offset of the T i on its interior. In other words, we consider a surface of revolution around the y-axis of a circle of radius 1 -δ in the xz-plane centred at (2, 0, 0), with an ε-neighbourhood of the circle found by revolving the point (1, 0, 0) around the y-axis cut out. For the positioning of the discrete set of points, we let C i,δ be the circle found by revolving the point (1 -δ, 0, 0) around the y-axis. As in Example 15, we describe the positioning of the discrete points in an inductive manner. Each pair of points, p i and pi , lies on C i,δ at a distance 2r i-1 from each other. We define

r 0 = δ 2 + ε 2 + δ(ε 2 -δ 2 + 2δ) 2 = ε 2 (1 + δ) + δ 2 (3 -δ) 2 .
The distance r 0 is chosen such that the balls B(p 1 , r) and B(p 1 , r) start to intersect at the same value of r as these balls start to intersect C 1 ⊕ B(r). The computation of this distance is based on Figure 11.

Figure 10 We depict the points p and p in red as well as their offset {p, p} ⊕ B(r) in purple, the points q and q as well as their offset {q, q} ⊕ B(r) are indicated in magenta.

The consecutive r i 's are defined as follows: Let q i and qi be the two points in the intersection of the bisector of p i and pi and the set C i that lie closest to p i and pi . Then we define r i to be the radius of the tetrahedron with vertices p i pi q i qi . Note that {q 1 , q1 } = π C1 p1+ p1 2

. We continue this construction inductively until r i > 1 -δ.

In the proof of Proposition 10 acuteness played an essential role. We argued that if

√ 2(1 -δ) -1 < ε, (12)
then any triangle p i pi q i (see Example 15) is acute. A triangle is acute if and only if it contains its circumcentre. For the manifold case we therefore generalize acuteness to tetrahedra in the following manner:

▶ Definition 17 (Self-centred simplices). A simplex is called self-centred if it contains its circumcentre.

▷ Claim 18. Let C = C i and C δ = C i,δ be as defined in Example 16. Let p, p ∈ C δ , and define q, q to be the two points in the intersection of the bisector of p, p and C δ that are closest to p (and thus p). If

ε 2 > 6 -4 √ 2 -2δ + δ 2 , ( 13
)
then the simplex ppq q is self-centred.

In the proofs of Claim 18 and Proposition 11 we need a number of distances, namely r 0 = w 2 , r -1 = h, and w ′ , as indicated in Figure 11. It is convenient to calculate these in one go. We first focus on the first two. The sample P contains sets C i which are tori with a band cut out, r -1 denotes half the size of the gap. As we will discuss in the proof of Proposition 11 there is a topological transition at this radius. The distance w ′ is important in establishing self-centredness. From Figure 11 we immediately see that

1 = ℓ 1 + ℓ 2 (14)
ℓ 2 1 + h 2 = ε 2 (15) ℓ 2 2 + h 2 = (1 -δ) 2 (16
)

w 2 = (δ + ℓ 1 ) 2 + h 2 . ( 17
)
Combining these yields ℓ 1 = ε 2 -δ 2 + 2δ 2 (by using [START_REF] Federer | Curvature measures[END_REF] to substitute ℓ 1 for ℓ 2 in ( 16) and then subtracting the result from ( 15)) ( 18)

h 2 = ε 2 - ε 2 -δ 2 + 2δ 2
2
(plugging ( 18) into ( 15)) ( 19)

w 2 = δ 2 + ε 2 + δ(ε 2 -δ 2 + 2δ).
(by combining ( 17), [START_REF] Lytchak | On the geometry of subsets of positive reach[END_REF], and ( 19))

Figure 11 also yields

w ′2 = (δ ′ + δ + ℓ 1 ) 2 + h 2 ,
from which one derives that

w ′2 = ε 2 + (1 -(1 -δ) cos ψ) 2 + (1 -(1 -δ) cos ψ)(ε 2 -δ 2 + 2δ), (20)
using ( 18), [START_REF] Lytchak | Almost convex subsets[END_REF], and the formula for δ ′ (see Figure 11 Proof of Claim 18. By definition, the circumcentre of the simplex ppq q lies on the bisector of p and p and the bisector of q and q. Thus, if the position of p and p on C δ varies, then the simplex ppq q transitions from being non-self-centred to self-centred exactly when its

≤ (R -α) 2 + r 2 -(R -δ) 2 ≤ 0,
showing that cos φ ≤ 0, or equivalently φ ≥ π 2 . Hence, p ′ lies on the segment qx. Let us show that ∥q -p ′ ∥ ≤ α. Because p ′ belongs to the segment qx ⊆ qz, we have

∥q -p ′ ∥ = ∥q -z∥ -∥p ′ -z∥ = R -∥z -p∥ 2 -∥p ′ -p∥ 2 ≤ R -∥z -p∥ 2 -∥x -p∥ 2 ≤ R -(R -δ) 2 -r 2 ≤ α = ∥q -y∥.
hence, y lies between x and p ′ . This shows that the distance to p decreases as we move on the segment xy, starting from x and going toward y. It follows that that xy ⊆ B(p, r) ⊆ P ⊕B(r). Since yq ⊆ B(q, α) ⊆ S ⊕ B(α) ⊆ P ⊕ B(r), we deduce that the whole segment xq belongs to P ⊕ B(r). The proof is completed by using the same deformation retract argument as in the first version of the proof. ◀

Figure 2 A

 2 Figure 2 A set S of positive reach with a sample P . The set p∈P B(p, r) = P ⊕ B(r) is also indicated.

Figure 3 A

 3 Figure 3 A pictorial overview of the proof.

Figure 5

 5 Figure 5 Overview of the notation used in the proof of Lemma 7.

Proposition 10 .

 10 Suppose that the dimension d of the ambient space R d satisfies d ≥ 2, and the one-sided Hausdorff distances ε and δ do not satisfy the bound[START_REF] Attali | Vietoris-rips complexes also provide topologically correct reconstructions of sampled shapes[END_REF]. Then there exists a set S of positive reach and a sample P that satisfy Universal Assumption 1, while the homology of the union of balls P ⊕ B(r) does not equal the homology of S for any r.▶ Proposition 11. Suppose that the dimension d of the ambient space R d satisfies d ≥ 3, the one-sided Hausdorff distances ε and δ do not satisfy the bound[START_REF] Boissonnat | Geometric and Topological Inference[END_REF], and δ ≤ ε. Then there exists a manifold M of positive reach and a sample P that satisfy Universal Assumption 1, while the homology of the union of balls P ⊕ B(r) does not equal the homology of M for any r.

Figure 6

 6 Figure 6 Situation in one annulus of Example 15.

( a )

 a For all r ≤ r0 the offset of the sample (Ci ∪ {pi, pi}) ⊕ B(r) has three connected components for each of the three annuli.(b) At the radius r1 the cycle in the thickening of the sample of the first annulus dies, while at the same radius a cycle is created for the second annulus.(c) At the radius r2 the cycle in the thickening of the sample of the second annulus dies, while at the same radius a cycle is created for the third annulus.

Figure 7

 7 Figure 7 In the figure we only depict three annuli Ai sorted from left to right.

◀

Figure 8

 8 Figure 8 The (half of the) torus Ti depicted in blue and the sample the set Ci together with the points pi and pi in red. In black we indicate the circle C i,δ . The closest point projection of this circle onto M is indicated in blue.

Figure 9

 9 Figure[START_REF] Chazal | A sampling theory for compact sets in Euclidean space[END_REF] We depict (half of) the sample P in red and its offset P ⊕ B(r) in purple.

Figure 11

 11 Figure 11 This figure is used for the computation of r0 = w 2 , r-1 = h, and w ′ .

  ∥z -x∥ = ∥z -q∥ -∥q -x∥ ≤ R -α. It follows that 2∥z -x∥∥x -p∥ cos φ = ∥z -x∥ 2 + ∥x -p∥ 2 -∥z -p∥ 2(reshuffling (22))

We recall that the one sided Hausdorff distance from X to Y , denoted by d o H (X; Y ), is the smallest ρ such that the union of balls of radius ρ centred at X covers Y .

Funding This research has been supported by the European Research Council (ERC), grant No. 788183, by the Wittgenstein Prize, Austrian Science Fund (FWF), grant No. Z 342-N31, and by the DFG Collaborative Research Center TRR 109, Austrian Science Fund (FWF), grant No. I 02979-N35. Mathijs Wintraecken: Supported by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 754411, and the Austrian science fund (FWF) grant No. M-3073 Acknowledgements We thank Jean-Daniel Boissonnat, Herbert Edelsbrunner, and Mariette Yvinec for discussion.

circumcentre lies on the edge pp. At this point, the circumcentre of the face ppq also lies on the edge pp. The same also holds for the face ppq. This is again equivalent to the fact that the distance w ′ between the midpoint p+ p 2 and q (and the midpoint p+ p 2 and q) equals ∥p-p∥ 2

. Because the points p and p lie on the circle C δ this distance can be written in terms of the angle ψ = ∠p0 p+ p 2 , that is, ∥p-p∥ 2 = (1 -δ) sin ψ. The (one-sided Hausdorff) distance from p+ p 2 to C δ then equals

The distance w ′ is given by [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF]. Hence the transition between self-centred and non-selfcentred takes place when

Writing x = cos ψ, this translates into

The solution space of this quadratic equation is non-empty if the discriminant

is non-negative. Introducing the variable y = ε 2 -δ 2 + 2δ and applying that δ < 1, the inequality ∆ ≥ 0 reduces to (2 + y) 2 -16y ≥ 0, and thus

The above expression is the same as Equation [START_REF] Edelsbrunner | Triangulating topological spaces[END_REF]. In summary, this implies that there is no transition from self-centred to non-self centred simplices as the positions of p and p vary on C δ if the inequality (13) is satisfied.

with c > 0, then the circumcentre of the simplex ppq q lies in its interior. This guarantees that the construction in Example 16 consists of a finite number of tori.

We can now finally prove the main statement of this section.

Proof of Proposition 11. We note that for r ∈ [0, r 1 ) the set (C 1 ∪ {p 1 p1 }) ⊕ B(r) does not have the same homology as the torus T 1 , because the tetrahedron p 1 p1 q 1 q1 is self-centred. The reason for this is that For r ∈ [0, r 0 ), the set (C 1 ∪ {p 1 , p1 }) ⊕ B(r) has three connected components.

For r ∈ [r 0 , r 1 ), the set (C 1 ∪ {p 1 , p1 }) ⊕ B(r) has homology type of a torus with either a circle or a 2-sphere attached (depending on if r is smaller or larger than the cirumradius of the triangle p 1 p1 q 1 -which is in turn equal to the circumradius of the triangle p 1 p1 q1 -).

We note there is also a topological transition at r

2 2 (r -1 is h in Figure 11) . The homotopy type of each C i ⊕ B(r) changes from that of a circle to that of a torus. ◀

A Federer's sets of positive reach

We recall the results from [START_REF] Federer | Curvature measures[END_REF] that we used in the main text. Here we just assume that S is a closed set.

▶ Theorem 24 ([14, Theorem 4.8 (12)]). If p ∈ S and lfs(p) > ℓ > 0, then

Tan(p, S) is the convex cone dual to Nor(p, S), and

for u ∈ Tan(p, S).

B Alternative proof

Alternative proof of Theorem 2. Let us prove that the set (q + Nor(q, S)) ∩ B(q, R) ∩ (P ⊕ B(r)) is star-shaped with respect to q. For this, consider a point x ∈ (q+Nor(q, S))∩B(q, R)∩ (P ⊕ B(r)) and let us prove that the segment xq is also contained in P ⊕ B(r). We consider two cases. First, suppose that ∥x-q∥ ≤ α. In that case, xq ⊆ B(q, α) ⊆ S ⊕B(α) ⊆ P ⊕B(r) and we are done. Second, that ∥x -q∥ > α; see Figure 12. In that case, x ̸ = q and the half-line with origin at q and passing through x is well-defined. Let y be the point on this half-line whose distance to q is α. Let z be the point on this half-line whose distance to q is R. Because x ∈ B(q, R), we have that x lies on the segment qz. Let p be any point of P whose distance to x is smaller than or equal to r. It is this assumption that later gives ∥x -p∥ ≤ r. Let p ′ be the projection of p onto the straight-line passing through q and x. We have that the five points x, y, z, q and p ′ are aligned and y lies between x and q. We claim that y also lies between x and p ′ . The claim is clearly true if q lies between x and p ′ . Let us assume that q does not lie between x and p ′ , in other words, let us assume that p ′ is on the half-line with origin at q and passing through x. Let φ be the internal angle of triangle xpz at x. The cosine rule gives: ∥z -p∥ 2 = ∥z -x∥ 2 + ∥x -p∥ 2 -2∥z -x∥∥x -p∥ cos φ. [START_REF] Rataj | Curvature measures of singular sets[END_REF] By Theorem 24, the interior of B(z, R) does not intersect S and because P ⊆ S ⊕ B(δ), we have ∥z -p∥ ≥ R -δ. By construction, we have that ∥x -p∥ ≤ r. Furthermore,


Van-Giang Trinh
email: trinh.van-giang@lis-lab.fr

Belaid -6581-998x]

Benhamou

Kunihiko Hiraishi

Sylvain Soliman
email: sylvain.soliman@inria.fr

Minimal trap spaces of Logical models are maximal siphons of their Petri net encoding

Keywords: Logical models, Boolean models, Trap spaces, Attractor computation, Petri nets, Siphons

Boolean modelling of gene regulation but also of post-transcriptomic systems has proven over the years that it can bring powerful analyses and corresponding insight to the many cases where precise biological data is not sufficiently available to build a detailed quantitative model. This is even more true for very large models where such data is frequently missing and led to a constant increase in size of logical models à la Thomas. Besides simulation, the analysis of such models is mostly based on attractor computation, since those correspond roughly to observable biological phenotypes. The recent use of trap spaces made a real breakthrough in that field allowing to consider medium-sized models that used to be out of reach. However, with the continuing increase in model-size, the state-of-the-art computation of minimal trap spaces based on prime-implicants shows its limits as there can be a huge number of implicants.

In this article we present an alternative method to compute minimal trap spaces, and hence complex attractors, of a Boolean model. It replaces the need for prime-implicants by a completely different technique, namely the enumeration of maximal siphons in the Petri net encoding of the original model. After some technical preliminaries, we expose the concrete need for such a method and detail its implementation using Answer Set Programming. We then demonstrate its efficiency and compare it to implicant-based methods on some large Boolean models from the literature.

Introduction

From the observation that the transcriptional regulation behaved in a sigmoid step-like way, came the original idea to represent models of gene regulation as discrete event systems. Those Gene Regulation Networks (GRN) use thresholds or equivalently logical functions to represent the different regulations [START_REF] Glass | The logical analysis of continuous, non-linear biochemical control networks[END_REF][START_REF] Thomas | Boolean formalisation of genetic control circuits[END_REF][START_REF] Thomas | Biological feedback[END_REF][START_REF] Thomas | Regulatory networks seen as asynchronous automata: a logical description[END_REF].

Boolean modelling has proven over the years that it can bring powerful analyses and corresponding insight to the many cases where precise biological data is not sufficiently available to build a detailed quantitative model [START_REF] Wang | Boolean modeling in systems biology: an overview of methodology and applications[END_REF], even for modelling post-transcriptional mechanisms. This is even more true for very large models where such data is frequently missing and led to a constant increase in size of logical models à la Thomas [START_REF] Aghamiri | Automated inference of Boolean models from molecular interaction maps using CaSQ[END_REF]. Besides simulation, the analysis of such models is mostly based on attractor computation, since those correspond roughly to observable biological phenotypes. The recent use of trap spaces [START_REF] Klarner | Computing maximal and minimal trap spaces of Boolean networks[END_REF] made a real breakthrough in that field allowing to consider medium-sized models that used to be out of reach. However, with the continuing increase in model-size, the state-of-the-art computation of minimal trap spaces based on prime-implicants shows its limits as there can be a huge number of implicants.

It is worth noting that the recent method presented in [START_REF] Chevalier | Synthesis of Boolean networks from biological dynamical constraints using answer-set programming[END_REF] for computing minimal trap spaces avoids the prime-implicants computation by relying on the most-permissive semantics of Boolean models. This method has been implemented in the tool mpbn 4 demonstrated in [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF] for handling medium-sized models from the literature and very large synthetic models (up to 100,000 nodes). However, this method is only applicable for locally-monotonic Boolean models, whereas the prime-implicants based method [START_REF] Klarner | Computing maximal and minimal trap spaces of Boolean networks[END_REF] is applicable for general Boolean models (i.e., including both locally-monotonic and non-locally-monotonic ones).

The study [START_REF] Noual | About non-monotony in Boolean automata networks[END_REF] highlights the need for non-locally-monotonic Boolean models in both biological and theoretical aspects. Hence, it is still necessary to develop efficient methods for computing minimal trap spaces of large-scale general Boolean models.

Petri nets were introduced in the 60s as 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][START_REF] Murata | Petri nets: Properties, analysis and applications[END_REF]. The use of Petri nets for representing biochemical reaction systems, by mapping molecular species to places and reactions to transitions, hinted at already in [START_REF] Peterson | Petri net theory and the modeling of systems[END_REF][START_REF] Murata | Petri nets: Properties, analysis and applications[END_REF] was used more thoroughly 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. Siphons are such a concept, but they have not been used a lot for the study of biochemical systems [START_REF] Zevedei-Oancea | Topological analysis of metabolic networks based on Petri net theory[END_REF][START_REF] Blätke | Biomodel engineering with Petri nets[END_REF] even if the practical cost of computing their minimal/maximal elements appear much more manageable than the theoretical complexity would indicate [START_REF] Oanea | New algorithms for deciding the siphon-trap property[END_REF][START_REF] Nabli | On enumerating minimal siphons in Petri nets using CLP and SAT solvers: theoretical and practical complexity[END_REF].

In this article we present an alternative method to compute minimal trap spaces, and hence complex attractors, of a general Boolean model. It replaces the need for prime-implicants by a completely different technique, namely the enumeration of maximal siphons in the Petri net encoding of the original model. Af-ter some technical preliminaries, we expose the concrete need for such a method and detail its implementation using Answer Set Programming. We then demonstrate its efficiency and compare it to implicant-based methods on several large Boolean models from the literature.

All models used for evaluation and the implementation of the presented method are available at https://github.com/soli/trap-spaces-as-siphons and executable as a CoLoMoTo docker image.

Preliminaries

We will briefly recall here some preliminaries on Boolean models related to trap spaces and Petri nets. In the case of multi-level Logical models, an encoding into a Boolean model is always possible [START_REF] Didier | Mapping multivalued onto Boolean dynamics[END_REF].

Traps spaces

We recall here some definitions from [START_REF] Klarner | Computing maximal and minimal trap spaces of Boolean networks[END_REF] for the introduction of trap spaces. Minimal trap spaces prove to be a very good approximation of the attractors of a Boolean model under asynchronous update schemes and have become the de facto standard way to analyze models of a few tens of genes [START_REF] Klarner | PyBoolNet: a python package for the generation, analysis and visualization of Boolean networks[END_REF][START_REF] Fontanals | Control strategy identification via trap spaces in Boolean networks[END_REF].

Given a Boolean model

M = (V, F ) with nodes V = (v 1 , . . . , v n ) and Boolean functions F = (f 1 , . . . , f n ), its state-space is S M = B n with B = {0, 1}. A state s ∈ B n
is a mapping s: V → B that assigns either 0 (inactive) or 1 (active) to each node. At each time step t, node v i can update its state by

v i (t + 1) ∈ {v i (t), f i (v(t))}
where v(t) is the state of M at time t, v i (t+1) is the state of node v i at time t+1 and an update scheme specifies which value is chosen [START_REF] Thomas | Regulatory networks seen as asynchronous automata: a logical description[END_REF]. There are two main types of update schemes: synchronous, where for all i, v i (t + 1) = f i (v(t)), and asynchronous, where an i is chosen among the possible updates (f i (v(t)) ̸ = v i (t)), and all the other nodes j remain fixed at v j (t).

A non-empty set T ⊆ S M is a trap set of M if and only if for every x ∈ T and y ∈ S M with y is reachable from x it holds that y ∈ T .

A subspace m of S M is characterized by its fixed nodes (denoted by D m ) and free nodes. This subspace can be specified by an assignment m : D m → B where D m ⊆ V and m(u) is the value of node u ∈ D m . The remaining nodes of M, V \ D m , are said to be free, i.e., they can receive any Boolean value. We write subspaces like states but use in addition the symbol ⋆ to indicate that a node is free. A subspace m thus corresponds to the set of states M [START_REF] Klarner | Computing maximal and minimal trap spaces of Boolean networks[END_REF]. A trap space is defined as a subspace that is also a trap set. It is noted that trap spaces of a Boolean model are independent of the update scheme of this model [START_REF] Klarner | Computing maximal and minimal trap spaces of Boolean networks[END_REF]. Then, we define a partial order < on S ⋆ M as:

S M [m] := {s ∈ S M | ∀v ∈ D m : s(v) = m(v)}. For example, m = ⋆ ⋆ 1 means that D m = {v 3 }, m(v 3 ) = 1,
m < m ′ if and only if S M [m] ⊆ S M [m ′ ] and S M [m] ̸ = S M [m ′ ]. Consequently,
= (V, F ), where V = (x 1 , x 2 ) and F = (f 1 , f 2 ) with f 1 = (x 1 ∧ x 2 ) ∨ (¬x 1 ∧ ¬x 2 ), f 2 = (x 1 ∧ x 2 ) ∨ (¬x 1 ∧ ¬x 2 ).
Herein, ∧, ∨, and ¬ denote the conjunction, disjunction, and negation logical operators, respectively. A marking for a Petri net is a mapping m : P → N that assigns a number of tokens to each place. A place p is marked by a marking m if and only if m(p) > 0. We shall write pred(x) (resp. succ(x)) to represent the set of vertices that have a (non-zero weighted) arc leading to (resp. coming from) x.

px 1 px 1 px 2 px 2 t 1 x 1 t 2 x 1 t 1 x 2 t 2
The link between Boolean models à la Thomas and Petri nets was originally established in [START_REF] Chaouiya | Qualitative modelling of genetic networks: From logical regulatory graphs to standard Petri nets[END_REF] in order to make available formal methods like model-checking for the analysis of such systems. The basic encoding into 1-safe (i.e., never more than one token in each place) nets only holds for purely Boolean models but was later extended to multivalued Logical models in two ways, either in [START_REF] Chaouiya | Petri net representation of multivalued logical regulatory graphs[END_REF] with non 1-safe Petri nets or more recently in [START_REF] Chatain | Characterization of reachable attractors using Petri net unfoldings[END_REF] with 1-safe nets but many more places.

Since our study is focused on Boolean models, we briefly recall the original encoding here. Its basis is that every node (gene) v of the original model M = (V, F ) is represented by two separate places (p v and p v ), corresponding to its two states, active, and inactive, respectively. Each conjunct of the logical function that activates the gene will lead to a transition t, consuming the inactive place (i.e., a directional arc from p v to t), producing the active place (i.e., a directional arc from t to p v ), and with all other literals both consumed and produced (i.e., a bidirectional arc). And conversely for the inactivation. Let s be a state of the Boolean model and m s be its corresponding marking in the encoded Petri net. It holds that ∀v ∈ V , s(v) = 0 if and only if m s (p v ) = 1 and s(v) = 1 if and only if m s (p v ) = 1. Note also that at any marking m of the Petri net encoding a Boolean model, it always holds that m(p v ) + m(p v ) = 1.

The main property of this encoding is that it is completely faithful with respect to the update scheme of the original Boolean model. For each node v of M, only transitions corresponding to v can change the current marking of p v or p v . In addition, at any marking at most one of such transitions is enabled because m(p v ) + m(p v ) = 1 holds. Hence, for any update scheme in M, we have a corresponding firing scheme in P, which preserves the equivalence between the dynamics of M and P [START_REF] Chatain | Concurrency in Boolean networks[END_REF].

For illustration, let us reconsider the Boolean model shown in Example 1. Figure 1 cannot be both enabled. Then, the fully asynchronous update scheme in M corresponds to the classical firing scheme in P where only one of the enabled transitions for a given marking will be fired [START_REF] Murata | Petri nets: Properties, analysis and applications[END_REF].

Note that given a Boolean model in the standard SBML-Qual format [START_REF] Chaouiya | SBML qualitative models: a model representation format and infrastructure to foster interactions between qualitative modelling formalisms and tools[END_REF], i.e., the package of SBML v3 [START_REF] Keating | SBML Level 3: an extensible format for the exchange and reuse of biological models[END_REF] for such models, one can easily obtain its Petri net encoding in the Petri Net Markup Language (PNML)5 standard using the BioLQM6 library. This piece of software extracted from GINsim [START_REF] Chaouiya | Logical modelling of gene regulatory networks with GINsim[END_REF] and part of the CoLoMoTo7 [START_REF] Naldi | Cooperative development of logical modelling standards and tools with CoLoMoTo[END_REF] software suite allows for easy conversion between standard formats. It also accepts many other common formats for Boolean models, notably the .bnet files of the BoolNet [START_REF] Müssel | BoolNet -an R package for generation, reconstruction and analysis of Boolean networks[END_REF][START_REF] Klarner | PyBoolNet: a python package for the generation, analysis and visualization of Boolean networks[END_REF] tools. The conversion is executed as follows:

java -jar GINsim.jar -lqm <input.{sbml,bnet,zginml,...}> <output.pnml> Note that transforming a Boolean model defined by its functions into its Petri net encoding roughly relies on obtaining conditions for the activation and inactivation of the states. In [START_REF] Chaouiya | Qualitative modelling of genetic networks: From logical regulatory graphs to standard Petri nets[END_REF] this took the form of the whole truth table of the Boolean functions, but as shown in Appendix 1 of [START_REF] Chatain | Characterization of reachable attractors using Petri net unfoldings[END_REF] computing Disjunctive Normal Forms (DNF) of each Boolean function is enough. Though this might appear quite computationally intensive it is important to remark first that contrary to the prime-implicants case, there is no need to find minimal DNFs. One way to look at this is to consider that this amounts to a similar approach as that used in [START_REF] Chevalier | Synthesis of Boolean networks from biological dynamical constraints using answer-set programming[END_REF] but with the encoding of both activation and inhibition functions as DNFs in order to take into account possible non-local-monotonicity. This does not change the worst-case-complexity (obtaining a single DNF being exponential) but might matter a lot in practice. As such, we will explore how this transformation, here using BDDs in BioLQM, and the one based on the most-permissive semantics compare in the Section 6 on evaluation.

Siphons

Siphons are a static and classical property of Petri nets [START_REF] Peterson | Petri net theory and the modeling of systems[END_REF]. Note however that the use of siphons for the analysis of biological models, though it is not new, has been mostly relevant to the ODE-based continuous semantics of Chemical Reaction Networks [START_REF] Angeli | A Petri net approach to persistence analysis in chemical reaction networks[END_REF][START_REF] Angeli | Persistence results for chemical reaction networks with time-dependent kinetics and no global conservation laws[END_REF][START_REF] Degrand | Graphical conditions for rate independence in chemical reaction networks[END_REF].

We recall here the basic definition establishing that to produce something in a siphon you must consume something from the siphon. This corresponds to the idea that a siphon is a set of places that once unmarked remains unmarked. Definition 2. A siphon of a Petri net (P, T, W ) is a set of places S such that:

∀t ∈ T, S ∩ succ(t) ̸ = ∅ ⇒ S ∩ pred(t) ̸ = ∅.
Note that ∅ is trivially a siphon.

Minimal trap spaces as maximal conflict-free siphons

First, we add a definition related to any set of places of a Petri net encoding a Boolean model, and notably a siphon of such a net. Definition 3. A set of places of Petri net P encoding Boolean model M is conflict-free if it does not contain any two places corresponding to the active and inactive states of the same gene of M. Then, a conflict-free siphon S is said to be maximal if and only if there is no other conflict-free siphon S ′ such that S ⊂ S ′ .

Intuitively, a siphon is a set of places that once unmarked remains so. If it is conflict-free then its dual corresponds to a partial-state of the model such that whatever update, the fixed values remain so (since the unmarked places remain unmarked). This is precisely the definition of a trap space and maximality of the siphon is equivalent to as many fixed values as possible, hence minimality of the trap space. For example, the Boolean model given in Example 1 has two trap spaces, m 1 = 11 and m 2 = ⋆⋆. The Petri net encoding of this Boolean model has five generic siphons,

S 1 = ∅, S 2 = {p x1 , p x1 }, S 3 = {p x2 , p x2 }, S 4 = {p x1 , p x2 },
and S 5 = {p x1 , p x1 , p x2 , p x2 }. However, only S 1 and S 4 are conflict-free siphons and correspond to m 2 and m 1 , respectively. Since S 1 ⊂ S 4 , S 4 is a maximal siphon corresponding to the minimal trap space m 1 . Hereafter, we formally prove that a maximal conflict-free siphon is equivalent to a minimal trap space. Definition 4. Let m be a subspace of Boolean model M = (V, F ). A mirror of m is a set of places S in the Petri net encoding P of M such that:

∀v ∈ D m , m(v) = 0 ⇔ p v ∈ S, m(v) = 1 ⇔ p v ∈ S and ∀v ∈ V \ D m , p v ̸ ∈ S, p v ̸ ∈ S.
Theorem 1. Let M = (V, F ) be a Boolean model and P be its Petri net encoding. A subspace m is a trap space of M if and only if its mirror S is a conflict-free siphon of P.

Proof. First, we show that if m is a trap space of M, then S is a conflict-free siphon of P (*). If D m = ∅, then S = ∅ is trivially a conflict-free siphon of P. Thus, we consider the case that D m ̸ = ∅ (resp. S ̸ = ∅). Assume that S is not a siphon of P. Then, there is a transition t ∈ T such that S ∩ succ(t) ̸ = ∅ but S ∩ pred(t) = ∅. In other words, there is a place p ∈ S such that p ∈ succ(t) but p ̸ ∈ pred(t). Let v be the corresponding node in M of p. By the characterization of the encoding [START_REF] Chaouiya | Qualitative modelling of genetic networks: From logical regulatory graphs to standard Petri nets[END_REF], there is a directional arc from t to p and a directional arc from the complementary place of p to t. Without loss of generality, we assume that p = p v , then there is a directional arc from t to p v and a directional arc from p v to t. In addition, there is also no arc or a bidirectional arc between t and another place rather than p v and p v . Thus, there is no connecting arc between t and any place in

S \ {p v } because S ∩ pred(t) ̸ = ∅. In S M [m], a node in V \ D m can receive any Boolean value. Hence, there is a state s ∈ S M [m] such that m s (p ′ ) = 1, ∀p ′ ∈ pred(t) \ {p v }
where m s is the corresponding marking in P of s. We also have m s (p v ) = 0, leading to m s (p v ) = 1 by the characterization of the encoding [START_REF] Chaouiya | Qualitative modelling of genetic networks: From logical regulatory graphs to standard Petri nets[END_REF]. Now, t is enabled at marking m s . Its firing leads to a new marking

m ′ s such that m ′ s (p v ) = 1 and m ′ s (p v ) = 0. Let s ′ be the corresponding state in M of m ′ s . Since m is a trap space of M, s ′ ∈ S M [m]. Then, s ′ (v) = m(v), leading to m ′ s (p v ) = 0
, which is a contradiction. Hence, S is a siphon of P. By the definition of a mirror, S is also a conflict-free one.

Second, we show that if S is a conflict-free siphon of P, then m is a trap space of M (**). By the definition of a mirror, m is a subspace of M. Let s be an arbitrary state in S M [m] and m s be its corresponding marking in P. By the characterization of the encoding [START_REF] Chaouiya | Qualitative modelling of genetic networks: From logical regulatory graphs to standard Petri nets[END_REF], m s (p) = 0, ∀p ∈ S. In any marking m ′ s reachable from m s regardless of the firing scheme of P, we have m ′ s (p) = 0, ∀p ∈ S by the dynamical property on markings of a siphon [START_REF] Liu | A survey of siphons in Petri nets[END_REF]. Equivalently, in any state s ′ reachable from s regardless of the update scheme of M, we have

s ′ (v) = s(v) = m(v), ∀v ∈ D m . Then, s ′ ∈ S M [m]
. By the definition of a trap space and the arbitrariness of s, m is a trap space of M.

From (*) and (**), we can conclude the proof. ⊓ ⊔ Theorem 2. Let M be a Boolean model and P be its Petri net encoding. A subspace m is a minimal trap space of M if and only if its mirror S is a maximal conflict-free siphon of P.

Proof. First, we show that if m is a minimal trap space of M, then S is a maximal conflict-free siphon of P (*). Since m is a trap space of M, S is a conflict-free siphon of P by Theorem 1. Assume that S is not maximal. Then, there is another conflict-free siphon S ′ such that S ⊂ S ′ . By Theorem 1, there is a trap space m ′ corresponding to S ′ . Following the definition of a mirror,

S M [m ′ ] ⊂ S M [m]
, thus m ′ < m. This is a contradiction because m is a minimal trap space. Hence, S is a maximal conflict-free siphon of P.

Second, we show that if S is a maximal conflict-free siphon of P, then m is a minimal trap space of M (**). Since S is a conflict-free siphon of P, m is a trap space of M by Theorem 1. Assume that m is not minimal. Then, there is another trap space m ′ such that m ′ < m. In other words, S M [m ′ ] ⊂ S M [m]. Let S ′ be the mirror of m ′ . S ′ is a conflict-free siphon by Theorem 1. Following the definition of a mirror, S ⊂ S ′ , which is a contradiction because S is a maximal conflict-free siphon. Hence, m is a minimal trap space of M.

From (*) and (**), we can conclude the proof.

⊓ ⊔

By Theorem 2, we can reduce the problem of computing all minimal trap spaces of a Boolean model to the problem of computing all maximal conflict-free siphons of its Petri net encoding. Note that in the case of stable states, this can be put in regard to the classical relationship between siphons and deadlocks in Petri nets. It might actually be possible to generalize our result to any 1-safe place-complementary Petri net to define a notion of trap space that might be useful for the analysis of Petri nets, but this is out of the scope of this article.

It is noted that there are no existing methods specifically designed for computing maximal conflict-free siphons (even maximal siphons) of a Petri net. The reason might be that researchers mainly focus on minimal siphons [START_REF] Liu | A survey of siphons in Petri nets[END_REF]. Hence, we here propose a new method based on Answer Set Programming (ASP) [START_REF] Gebser | Potassco: The Potsdam answer set solving collection[END_REF] for computing maximal conflict-free siphons of a Petri net. The details of the proposed method shall be given in the next section.

Answer set programming-based method

First, we show the characterization of all conflict-free siphons of the encoded Petri net P = (P, T, W ). Suppose that S is a generic siphon of P. If a place p should belong to S, then by definition all the transitions in pred(p) must belong to succ(S). Note that succ(S) = p∈S succ(p). A transition t belongs to succ(S) if and only if there is at least one place p ′ in S such that p ′ ∈ pred(t). Hence, for each transition t ∈ pred(p), we can state that

p ∈ S ⇒ p ′ ∈pred(t) p ′ ∈ S. (1)
The system of all the rules of the above form with respect to all pairs (p, t) where p ∈ P, t ∈ T, t ∈ pred(p) fully characterizes all generic siphons of a Petri net and has been used with SAT solvers in [START_REF] Oanea | New algorithms for deciding the siphon-trap property[END_REF][START_REF] Nabli | On enumerating minimal siphons in Petri nets using CLP and SAT solvers: theoretical and practical complexity[END_REF]. To make S to be a conflict-free siphon, we need to add to the system the rule

p v ∈ S ⇒ p v ̸ ∈ S ∧ p v ∈ S ⇒ p v ̸ ∈ S (2)
for each node v ∈ V . By definition, the final system fully characterizes all conflictfree siphons of the encoded Petri net. Then, we translate the above characterization into the ASP L as follows. We introduce atom p-v (resp. n-v) to denote place p v (resp. p v ), ∀v ∈ V . The set of all atoms in L is given as A = v∈V {p-v, n-v}. For each pair (p, t) where p ∈ P, t ∈ T, t ∈ pred(p), we translate the rule (1) into the ASP rule a_1; ... ; a_k :-a.

where a ∈ A is the atom representing place p and {a_1, . . . , a_k} ⊆ A is the set of atoms representing places in pred(t). The rule ( 2) is translated into the ASP rule :-p-v, n-v. for each v ∈ V . This ASP rule guarantees that two places representing the same node in M never belong to the same siphon of P, representing the conflictfreeness. Naturally, a Herbrand model (see, e.g, [START_REF] Gebser | Potassco: The Potsdam answer set solving collection[END_REF]) of L is equivalent to a conflict-free siphon of P. To guarantee that a Herbrand model is also a stable model (an answer set), we need to add to L the two choice rules {p-v}. {n-v}.

for each v ∈ V . Note that the number of atoms of L is only 2n, whereas the ASP encoding shown in [START_REF] Klarner | Computing maximal and minimal trap spaces of Boolean networks[END_REF] has as many atoms as the number of prime-implicants of the Boolean model and that number might be exponential in n. In [START_REF] Chevalier | Synthesis of Boolean networks from biological dynamical constraints using answer-set programming[END_REF], there is an ASP characterization of trap spaces that does not rely on minimal DNFs either and thus seems very similar to our ASP encoding. Remarkably it only requires the DNF for the activation part, using the information that it will only be used for locally-monotonic Boolean models. We would therefore expect that, when available, it will have comparable performance on the ASP part (the ASP program would be approximately twice smaller, though redundancy is not always bad in that field), but can also avoid combinatorial explosion of the Petri net encoding for some formula where the activation DNF is simple but the inhibition is not. Since mpbn is included in our benchmark this will be evaluated in our experiments. Now, a solution (simply an answer set) A ⊆ A of L is equivalent to a conflictfree siphon S of P, thus a trap space m of M. The conversion from

A to m is straightforward. If p-v ∈ A then v ∈ D m and m(v) = 0. Conversely, if n-v ∈ A then v ∈ D m and m(v) = 1. Otherwise, v ̸ ∈ D m .
Computing multiple answer sets is built into ASP solvers and the solving collection POTASSCO [START_REF] Gebser | Potassco: The Potsdam answer set solving collection[END_REF] also features the option to find set-inclusion maximal answer sets with respect to the set of atoms. Naturally, a set-inclusion maximal answer set of L is equivalent to a maximal conflict-free siphon of P, thus a minimal trap space of M. By using this built-in option, we can compute all the set-inclusion maximal answer sets of L (resp. all the minimal trap spaces of M) in one execution.

Motivating example

For a few years now we have been collaborating with biologists who build very large detailed and annotated maps and now wish to analyze the dynamics of the corresponding models. One of the main maps studied this way represents knowledge about the Rheumatoïd Arthritis [START_REF] Singh | Computational systems biology approach for the study of rheumatoid arthritis: from a molecular map to a dynamical model[END_REF], and was the main motivation for the development of a tool to automatically transform it into an executable Boolean model [START_REF] Aghamiri | Automated inference of Boolean models from molecular interaction maps using CaSQ[END_REF]. In the supplementary material of the paper, an excerpt of the map, focused around the apoptosis (cell death) module is transformed into a model of reasonable size, namely 180 Boolean variables (model F5_RA_apoptosis_executable_module.sbml of supplementary material S3, and model "RA-apoptosis" of Section 6). The study of such model, though, is a big hurdle. Indeed, as stated in the article about another model of the same size: "The size of the CaSQ-inferred MAPK model (181 nodes) made the calculation of stable states a non-realistic endeavour."

In practice, even if there is a huge number of attractors in such a model, obtaining a sample of those can reveal very useful to invalidate the model and lead to further refinement. In particular, it provides a feature-rich alternative to random simulations for this type of very non-deterministic model. Being able to detect that there are inconsistencies with published experimental data in some of the first 1000 attractors, for instance, can lead to a much quicker Systems Biology loop: model, invalidate, refine.

However, using a state-of-the-art tool like PyBoolNet [START_REF] Klarner | Computing maximal and minimal trap spaces of Boolean networks[END_REF] on that model actually fails at the phase of prime-implicant generation. mpbn [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF] does not give any answer either because it recognizes that model as non-locally-monotonic. And hence, it is not possible to extract any (complex) attractor at all. This is also true for the Alzheimer model also mentioned in that same article and originally from [START_REF] Ogishima | AlzPathway, an updated map of curated signaling pathways: towards deciphering Alzheimer's disease pathogenesis[END_REF] (F4 file in the original supplementary material, and "Alzheimer" in Table 2), but actually not for the MAPK model for which the first trap spaces can be obtained in reasonable time. The current practice usually revolves then around fixing some inputs to plausible values and reducing the model accordingly. While this approach makes sense, it relies on potentially arbitrary decisions, and hides away critical modelling choices that were actually not part of the original Boolean model or even of the starting map.

Using the method presented above, it is possible to convert the model to PNML in about one second and to obtain the first 1000 minimal trap spaces (including ones that contain more than one state) in a few milliseconds. Unfortunately since this was not available at the time, the analysis of the model remained very high-level and qualitative, instead of being able to use the rich information of computed minimal trap spaces.

Evaluation

To assess the efficiency of the proposed method, implemented as a Python package named Trappist, we compare it with the state-of-the-art method implemented in the tool PyBoolNet [START_REF] Klarner | Computing maximal and minimal trap spaces of Boolean networks[END_REF][START_REF] Klarner | PyBoolNet: a python package for the generation, analysis and visualization of Boolean networks[END_REF] on both its own repository of models and large models from the literature. In addition, we also include the tool mpbn [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF] to the benchmarks although it only handles locally-monotonic models, whereas both Trappist and PyBoolNet can handle general models.

To solve the ASP problems, we used the same ASP solver CLINGO [START_REF] Gebser | Potassco: The Potsdam answer set solving collection[END_REF] and the same configuration as that used in PyBoolNet [START_REF] Klarner | Computing maximal and minimal trap spaces of Boolean networks[END_REF][START_REF] Klarner | PyBoolNet: a python package for the generation, analysis and visualization of Boolean networks[END_REF]. Specifically, we used the configuration -heuristic=Domain -enum-mod=domRec -dom-mod=3 (subset maximality, equivalent to the deprecated -dom-pref=32 -heuristic=domain -dom-mod=7 used by PyBoolNet). We ran all the benchmarks on an apple laptop whose environment is CPU: Intel® Core™ i7 1.20GHz x 4, 16 GB DDR4 RAM, MacOS 12.3.1. Finally, we set a time limit of two minutes for each model. Note that we did not get the opportunity to run this benchmark on a proper dedicated computing workstation, the results are therefore only indicative of global trends and should not be interpreted as precise performance. This is also why in some rare cases finding all minimal trap spaces was faster than with the limit to 1000. For all the above tools, namely PyBoolNet, mpbn and CLINGO, we used the version available in the latest CoLoMoTo docker image tagged 2022-05-01. All the models and a CoLoMoTo notebook realizing the benchmarks can be found at https://github.com/soli/trap-spaces-as-siphons. These can be run on a Docker image in the cloud by clicking the "Binder" button.

PyBoolNet repository

As shown in Table 1, for most of the models of the official PyBoolNet repository8 , the results are comparable with all minimal trap spaces found very fast. For 5 of the 29 models, mpbn did not give any answer because it recognized these models as not locally-monotonic. Note that on some very small models, Trappist is sometimes slower than PyBoolNet and/or mpbn, but still significantly under one second. Moreover, we believe that the result on the arellano_rootstem model is caused by the cold start of the JVM for BioLQM. On the contrary, on every model that was a bit challenging for PyBoolNet or mpbn, the new method is far more efficient with speedups between one and two orders of magnitude.

Selected models

We used a set of real-world Boolean models lying in various scales collected from numerous bibliographic sources. These models are quite big (in size), complex (i.e., having high average in-degree, which is related to the number of primeimplicants) and most of them have never been fully analyzed. We then applied PyBoolNet, mpbn, and Trappist to computing minimal trap spaces of these realworld models. It is notable that unlike existing analysis shown in the literature, we did not fix specific values for source nodes (i.e., some node v such that f v = v) in these models. Table 2 shows the experimental results on those models. Hereafter, we analyze in detail the results with respect to minimal trap space computation. The first observation is that for 26 of the 33 models (more than 78%), mpbn did not give any answer because it recognized that these models as not locallymonotonic. For 6 of the 33 models where mpbn returned the answers, mpbn and Trappist are comparable in computation time, though surprisingly mpbn appears a bit slower on average. Note however that mpbn was the only tool to provide a solution for the SN-5 model, thus confirming that if the activation function is in the right form, not having to compute the inactivation function's disjunctive normal form can render a difficult problem tractable. However, since Table 2: Timing comparisons between PyBoolNet (PBN), mpbn and Trappist on selected models from the literature. Column n (resp. s) denotes the number of nodes (resp. source nodes) of each model. Column |M | denotes the number of minimal trap spaces and for each method is given the computation time in seconds. "DNF" means that the method did not finish the computation (stopping at the first 1000 minimal trap spaces, or all for the last column) within the timeout of two minutes. "NM" indicates a non-locally-monotonic model. [START_REF] Kim | Reduction of complex signaling networks to a representative kernel[END_REF] 1953 669 >1000 DNF 23.62 5.93 DNF 32 SN-5 [START_REF] Kim | RMOD: A tool for regulatory motif detection in signaling network[END_REF] 2746 829 >1000 DNF 28.94 DNF DNF 33 turei-2016 [START_REF] Lee | Signal flow control of complex signaling networks[END_REF] 4691 1257 ??? DNF DNF DNF DNF mbpn can handle only locally-monotonic models and Trappist can handle general models, it is difficult to further compare between them. Hence, we focus on only comparisons between PyBoolNet and Trappist in the following observations. The second observation is that the proposed method vastly outperforms Py-BoolNet in computational time, on each and every model, and sometimes with orders of magnitude of difference (e.g., for most models in the 100-1000 nodes size range). Note that for all the cases where PyBoolNet did not manage to finish before the timeout, as marked by "DNF" in Table 2, the timeout occurred during the computation of the prime-implicants. Hence, not even a single minimal trap space was output by that method. The computational advantage is therefore immediately a practical advantage since on the one hand the state-of-the-art method did not allow any analysis whatsoever of the models, and on the other hand the proposed method could provide, very often under one second, the first thousand minimal trap spaces. For modellers having a critical look at a model and in a model, invalidate, refine loop this means a huge difference in the models that are amenable to study.

PBN mpbn

Note that even with a very restricted time-limit of two minutes, it was possible with the proposed technique to find all minimal trap spaces of small models (roughly under 130 nodes, i.e., considered as quite big up to now). Though it might seem impractical to handle tens of thousands of such possible complex attractors in a manual way, i.e., to compare them to specific experimental conditions and corresponding data, we hope that an automatic analysis of such attractors might become possible with systematic verification methods, not unlike that described in [START_REF] Hernandez | Computational verification of large logical models-application to the prediction of T cell response to checkpoint inhibitors[END_REF]. Since the ASP code is declarative by nature, it is also possible to add to it supplementary constraints coming from the modeler in case one is looking for specific attractors. Finally, sampling from the ASP-generated solutions as is done in [START_REF] Chevalier | Synthesis and simulation of ensembles of Boolean networks for cell fate decision[END_REF] would allow for a different type of exploration.

The third observation is that for all the models where PyBoolNet finished before the timeout, once PyBoolNet went through the prime-implicant phase, its ASP solving phase quickly returned the first 1000 minimal trap spaces, all under one second. For these models, the ASP solving phase of the proposed method also took very short time, all under one second. Hence, with the experimental results shown in this paper, the practical differences between our ASP encoding and that of PyBoolNet are not distinctly exposed. The fact that our new ASP encoding is guaranteed to be linear in the number of nodes of the original model does not seem to be crucial here, however a much deeper analysis of those cases remains to be done.

The last observation is that for very large models (i.e., more than two thousand nodes) the proposed method did not manage to finish the Petri net conversion before the timeout, as marked by "???" in Table 2. This points to the fact that our current choice of using a BDD-based translation to obtain that Petri net encoding, though it provides a small/efficient ASP might be too costly to handle the largest models. In such a case, a more naive encoding might provide a much larger ASP program, with many redundant rules, but easier/faster to obtain. The evaluation of the feasibility of such strategy, and of its impact on smaller instances, remains to be done and is out of the scope of this first article.

Recognizing that a model is locally-monotonic and applying in that specific case dedicated strategies as those of mpbn might also be a partial solution as shown in the SN-5 case, but not for the turei-2016 model.

Note that though enumerating the extremal siphons of a Petri net is exponential (see [START_REF] Nabli | On enumerating minimal siphons in Petri nets using CLP and SAT solvers: theoretical and practical complexity[END_REF] for instance) this is apparently not the bottleneck of the proposed method, showing once again that networks obtained from biochemical models do have a specific structure.

Conclusion

In this article we proposed a new method for the computation of minimal trap spaces of Boolean models, based on a new concept called maximal conflict-free siphons. This method is evaluated on large models from the literature and shows that it can scale up much better than the state-of-the-art prime-implicants based techniques. We believe that this opens up the way to a much better analysis of large Boolean models, which is needed with the advent of automatic modelgeneration pipelines [START_REF] Ostaszewski | COVID19 Disease Map, a computational knowledge repository of virus-host interaction mechanisms[END_REF].

Though we benchmarked this new approach against state-of-the-art tools, there are many more evaluations that we plan to do in the future. First, the Bi-oLQM platform that we are using is providing another implicant-based method using BDDs in http://colomoto.org/biolqm/doc/tools-trapspace.html and though we expect it to behave mostly like PyBoolNet, that remains to be checked. Note that this also raises the question of replacing altogether the BioLQM preprocessing step we use to obtain the Petri net encoding, since its use of BDDs might not be optimal for that step. The trade-off between a small Petri net, and hence small ASP, and the time it takes to compute it (in other words, the trade-off of allowing redundant constraints) has to be evaluated in depth.

Second, the experimental results shown in this paper mostly expose the differences caused by the prime-implicant phase of PyBoolNet. There is much more information required to distinctly study the practical differences between our ASP encoding and that of PyBoolNet, and notably their size/efficiency ratio. Hence, we plan to conduct experiments on more real-world models or maybe random models that can be randomly generated by using BoolNet [START_REF] Müssel | BoolNet -an R package for generation, reconstruction and analysis of Boolean networks[END_REF]. Such experiments should include a time-course of the number of ASP solutions found for proper comparison of the ASP encodings.

In addition, there are possibly other methods for computing maximal conflictfree siphons in Petri nets, like SAT/MaxSAT approaches [START_REF] Nabli | On enumerating minimal siphons in Petri nets using CLP and SAT solvers: theoretical and practical complexity[END_REF]. Although these approaches do not directly support the maximal conflict-free siphon computation now, we plan to investigate them in the future. They could replace our ASP program if they outperform it. However, the current method appears to already perform very well even on the biggest models we have considered.

Finally, we think that the links between Petri nets and Boolean models that we stumbled upon in this method might have deeper roots. Exploring those connections might lead both to interesting topics of research for Petri nets, like a notion of trap-spaces, and for Boolean models.

  and corresponds to the set of states {001, 011, 101, 111}. Let S ⋆ M denote the set of all possible subspaces of M. Note that |S ⋆ M | = 3 n and S M ⊂ S ⋆

Example 1 .

 1 a trap space m is minimal if and only if there is no trap space m ′ ∈ S ⋆ M such that m ′ < m. For example, let us consider the Boolean model shown in Example 1. Figure 1(a) shows the dynamics of this model under the fully asynchronous update (i.e., only one node is nondeterministically selected in order to be updated at each time step). The model has all two trap spaces, m 1 = 11 and m 2 = ⋆⋆. Since m 1 < m 2 , m 1 is a minimal trap space of the Boolean model. We give a Boolean model M

x 2 (

 2 b) Petri net encoding of the model. Circles denote places, whereas rectangles denote transitions.

Fig. 1 :

 1 Fig. 1: Dynamics and encoding of the Boolean model of Example 1

  (b) shows the Petri net encoding of this Boolean model. Place p x1 (resp. p x1 ) in P represents the activation (resp. the inactivation) of node x 1 in M. Marking {p x1 , p x2 } in P represents state 10 in M. Transitions t 1 x1 and t 2 x1 represent the update of node x 1 . Of course, in any marking t 1 x1 and t 2 x1

Table 1 :

 1 Timing comparisons between PyBoolNet, mpbn and Trappist on the PyBoolNet repository. Column n denotes the number of nodes of each model. Column |M | denotes the number of minimal trap spaces and for each method is given the computation time in seconds, asking only for the first 1000 trap spaces. A number in bold indicates a ratio greater than three compared to the best result. "NM" indicates a non-locally-monotonic model.

	model	n	|M | PyBoolNet mpbn Trappist
	1 arellano_rootstem	9	4	0.05	0.01	0.20
	2 calzone_cellfate	28	27	0.03 NM	0.03
	3 dahlhaus_neuroplastoma 23	32	0.05	0.02	0.03
	4 davidich_yeast	10	12	0.04	0.01	0.02
	5 dinwoodie_life	15	7	0.03	0.01	0.01
	6 dinwoodie_stomatal	13	1	0.03	0.01	0.01
	7 faure_cellcycle	10	2	0.04	0.01	0.02
	8 grieco_mapk	53	18	0.04	0.02	0.02
	9 irons_yeast	18	1	0.05	0.01	0.02
	10 jaoude_thdiff	103 >1000	1.44 0.90	0.09
	11 klamt_tcr	40	8	0.04	0.01	0.03
	12 krumsiek_myeloid	11	6	0.03	0.01	0.01
	13 multivalued	13	4	0.03	0.01	0.01
	14 n12c5	11	5	35.16	0.01	0.02
	15 n3s1c1a	2	2	0.02	0.01	0.01
	16 n3s1c1b	2	2	0.02	0.01	0.01
	17 n5s3	4	3	0.03 NM	0.01
	18 n6s1c2	5	3	0.03	0.01	0.01
	19 n7s3	6	3	0.02	0.01	0.01
	20 raf	3	2	0.02	0.01	0.01
	21 randomnet_n15k3	15	3	0.03 NM	0.01
	22 randomnet_n7k3	7	10	0.03 NM	0.01
	23 remy_tumorigenesis	34	25	2.14	0.02	0.02
	24 saadatpour_guardcell	13	1	0.03	0.01	0.01
	25 selvaggio_emt	56 >1000	1.02 0.52	0.09
	26 tournier_apoptosis	12	3	0.04	0.01	0.01
	27 xiao_wnt5a	7	4	0.03	0.01	0.01
	28 zhang_tlgl	60	156	0.22 NM	0.05
	29 zhang_tlgl_v2	60	258	0.09	0.15	0.02

https://www.pnml.org/

http://www.colomoto.org/biolqm/

http://colomoto.org/

https://github.com/hklarner/pyboolnet/tree/master/pyboolnet/repository


Julie Baltenneck

Ben Moussa

Loïc Couderc

Areski Flissi

Jacques Pédron

Cindy E Morris

Marie-Anne Barny

Odile Berge

Hélène Touzet

Erwan Gueguen

Tn-phyto: essential genomes of nine bacterial phytopathogens

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Dickeya dadantii 3937 Dickeya zeae A586 Dickeya fangzhongdai B16 Pectobacterium atrosepticum SCRI1043 Pectobacterium versatile A73 Pseudomonas syringae CC440 Pseudomonas syringae CC1557 chromosome

					Results
	SPECIES	REPLICATE	INITIAL NUMBER OF READS	NUMBER OF UNIQUELY MAPPED READS	TAs	DENSITY	MEAN READ	# ESSENTIAL GENES /
		1	11 021 397	10 148 785	171 791	50.1	57.3	435/4687
		2	11 657 925	10 699 150	171 791	50.3	60.4
		1	11 420 214	10 996 640	194 008	52.6	56.6	444/4475
		2	11 814 804	11 276 490	194 008	53.0	58.0
		1	14 652 216	13 811 935	166 816	38.1	82.6	447/4530
		2	11 834 554	11 157 191	166 816	37.9	66.7
		1	9 100 581	8 638 648	232 639	45.9	35.4	441/4733
		2	10 454 137	9 843 266	232 639	46.5	40.3
		1	11 354 590	10 579 502	219 900	55.5	55.1	407/5468
		2	10 215 789	9 512 400	219 900	55.3	49.6
	Pectobacterium aquaticum A212	1 2	12 091 457 11 184 408	10 283 920 9 460 081	199 301 199 301	62.7 62.6	49.0 45.1	480/4344
		1	15 818 299	14 860 409	132 365	64.1	112.2	429/5206
		2	12 946 819	12 215 317	132 365	63.8	92.2
		1	10 445 631	9 718 582	131 221	65.5	74.5	456/5244
		2	10 820 446	10 146 505	131 221	65.6	77.8
	Pseudomonas syringae TA043	1 2	9 425 329 7 475 968	8 883 550 6 963 428	121 303 121 303	70.2 69.7	73.1 57.3	431/4756

TOTAL GENES Funding: ANR-19-CE35-0016

Conclusions and perspectives

The aim of our project is to identify virulence factors of phytopathogenic bacteria isolated from the environment as well as those allowing their survival in water. For this purpose, we use Tn-seq (2), an innovative mutagenesis technique coupled with high throughput sequencing chosen to rapidly identify, on a large panel of strains, genes necessary for growth in a given condition, the plant or water. This study is carried out on strains of Dickeya, Pectobacterium and on the Pseudomonas syringae complex, all of which are capable of growing on a wide range of hosts as well as in water, thus allowing the most complete identification of genes of interest. Another axis of this ANR is the development of a complete bioinformatics pipeline, allowing to process data obtained under multiple conditions and for strains of different species. The results obtained will help us to predict the virulence of strains and/or to develop antibacterial strategies. Overall, this project will pave the way for studies of other environmental pathogens by proposing an original conceptual and technological framework that is extremely promising.

TN-PHYTO: ESSENTIAL GENOMES OF NINE BACTERIAL PHYTOPATHOGENS

Baltenneck J. 1 , Ben Moussa H. 2 , Couderc L. 3 , Flissi A. 5 , Pédron J. 2 , Morris C.E 4 , Barny MA. 2 , Berge O. 4 , Touzet H. 5 , Gueguen E 1* , Condemine G* 1 .

Conditions to apply

Materials and methods

Mutant libraries are constructed by conjugation between an Escherichia coli MDFpir strain and each of the strains of interest. During this conjugation, the plasmid pSAM-Ec, which contains the mariner Himar9 transposon that confers kanamycin resistance (3), a transposase and an ampicillin resistance gene, is transferred to the strain of interest. It is a suicide plasmid that cannot replicate in the recipient strain. Thus, the transposon will integrate uniquely and randomly into the genome of the strain of interest. The mutants obtained can then be selected on medium containing kanamycin.

The quality of the libraries created for each strain was analyzed by Novaseq sequencing (MGX, Montpellier).

It allowed to determine the percentage of transposon insertion in the genome of each strain. A good quality library is considered to have an insertion density of more than 30%. The correlation coefficient shows a high reproducibility of the transposition. These two elements attest to the good quality of the mutant libraries constructed.

(A) Now that the quality of the library is checked, the objective is to apply different conditions to the mutant libraries in order to determine the essential genes for bean infection or bacterial survival in a specific medium (TSB 50%, sterilized river water).