Louis Fippo Fitime

Olivier Roux

Carito Guziolowski
email: carito.guziolowski@irccyn.ec-nantes.fr

Loïc Paulevé
email: loic.pauleve@lri.fr

Identification of Bifurcations in Biological Regulatory Networks using Answer-Set Programming

Aiming at assessing differentiation processes in complex dynamical systems, this paper focuses on the identification of states and transitions that are crucial for preserving or pre-empting the reachability of a given behaviour. In the context of non-deterministic automata networks, we propose a static identification of so-called bifurcations, i.e., transitions after which a given goal is no longer reachable. Such transitions are naturally good candidates for controlling the occurrence of the goal, notably by modulating their propensity. Our method combines Answer-Set Programming with static analysis of reachability properties to provide an under-approximation of all the existing bifurcations. We illustrate our discrete bifurcation analysis on several models of biological systems, for which we identify transitions which impact the reachability of given long-term behaviour. In particular, we apply our implementation on a regulatory network among hundreds of biological species, supporting the scalability of our approach.

Introduction

The emerging complexity of dynamics of biological networks, and in particular of signalling and gene regulatory networks, is mainly driven by the interactions between the species, and the numerous feedback circuits it generates [START_REF] Thomas | Biological Feedback[END_REF][START_REF] Plathe | Feedback loops, stability and multistationarity in dynamical systems[END_REF][START_REF] Richard | Negative circuits and sustained oscillations in asynchronous automata networks[END_REF][START_REF] Paulevé | Static analysis of boolean networks based on interaction graphs: a survey[END_REF]. One of the prominent and fascinating features of cells is their capability to differentiate: starting from a multi-potent state (for instance, a stem cell), cellular processes progressively confine the cell dynamics in a narrow state space, an attractor. Deciphering those decisional processes is a tremendous challenge, with important applications in cell reprogramming and regenerative medicine.

Discrete models of network dynamics, such as Boolean and multi-valued networks [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF][START_REF] Bernot | Semantics of biological regulatory networks[END_REF], have been designed with such an ambition. These frameworks model nodes of the network by variables with small discrete domains, typically Boolean. Their value changes over time according to the state of their parent nodes. Exploring the dynamical properties of those computational models, such as reachability (the ability to evolve to a particular state) or attractors (long-run behaviors), allows to understand part of important cellular processes [START_REF] Sahin | Modeling ERBB receptor-regulated G1/S transition to find novel targets for de novo trastuzumab resistance[END_REF][START_REF] Abou-Jaoudé | Model checking to assess t-helper cell plasticity[END_REF][START_REF] Calzone | Predicting genetic interactions from boolean models of biological networks[END_REF].

In classical control theory the reachability and safety control of discrete systems is an important topic which is both critical and challenging. Safety verification or reachability analysis aims at showing that starting from some initial conditions a system cannot evolve towards to some unsafe region in the state space. In a stochastic setting the different trajectories originating from one initial state have a different likelihood and one can then evaluate what is the probability that the system reaches the assigned set of states starting from a given initial distribution of the set of initial states. In safety problems, where the evolution of the system can be influenced by some control input, one should select it appropriately so as to minimize the probability that the state of the system will move to an unsafe set of states. In this work we address the computation of sets of states from which the system can evolve into an undesired set of states given a model represented as a discrete finite-state of interacting components, such as an Automata Network. This topic is in particular highly relevant for systems biology and systems medicine, since it can suggest intervention sites or therapeutic targets in Biological Regulatory Networks (BRNs) that may counteract pathological behavior.

Contributions. In this work we introduce the notion of bifurcation in Automata Networks (ANs) and provide a scalable method for their identification relying on declarative programming with Answer-Set Programming (ASP) [START_REF]Knowledge Representation, Reasoning and Declarative Problem Solving[END_REF]. A bifurcation corresponds to a transition after which the system looses the capability to reach a given goal state. Identifying bifurcations extensively relies on the reachability problem, which is PSPACE-complete in ANs and related frameworks [START_REF] Cheng | Complexity results for 1-safe nets[END_REF]. In order to obtain an approach tractable on large biological networks, we show how to combine techniques from static analysis of ANs dynamics, from concurrency, and from constraint programming in order to relax efficiently the bifurcation problem. Our method identifies correct bifurcations only (no false positives), but, due to the embedded approximations, is incomplete (false negatives may exist). To our knowledge, this is the first integrated method to extract such decisive transitions from models of interaction networks. Related work. Many works have been devoted to the control theory and system verification; in particular in reachability analysis and safety control for system verification. We limit here to those developed for discrete systems.

In the case of deterministic systems, a large number of methods have been proposed to aid in the verification of safety-critical computer software and hardware. Due to the inherently discrete nature of these systems, early attempts in this area have concentrated on purely discrete systems [START_REF] Cassandras | Discrete event systems: modeling and performance analysis[END_REF]. A verification procedure has been proposed for discrete event dynamic systems in [START_REF] Moon | Automatic verification of discrete chemical process control systems[END_REF]. They use a finite state machine representation to model discrete event dynamic systems, temporal logic statements (to represent specifications imposed on their operations) and a model-checking verification method introduced by Clarke et al. [START_REF] Clarke | Automatic verification of finitestate concurrent systems using temporal logic specifications[END_REF] to test if the given relationship between the model and the specifications holds. The same method has been applied to the verification of programmable logic controllers [START_REF] Moon | Automatic verification of sequential control systems using temporal logic[END_REF] and control procedures for batch plants [START_REF] Moon | Formal verification of batch processing control procedures[END_REF].

For systems with complex behavior, the development of modeling frameworks and verification methods has been done among others by [START_REF] Grossman | Hybrid systems[END_REF]. Because of the significant increase of the state space, many approximation methods were proposed for reachability computation and can be grouped in two main approaches. The first one is based on an approximate simulation relation to obtain an abstraction of the original system [START_REF] Girard | Approximate simulation relations for hybrid systems[END_REF]. In the second approach, static analysis methods inspired from Cousot et al. [START_REF] Cousot | Abstract interpretation: A unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF] analyze the system dynamic and cope with the state-space explosion. In this group we can cite the development of static analysis of reachability properties based on an over-and under-approximation using the Process Hitting [START_REF] Paulevé | Static analysis of boolean networks based on interaction graphs: a survey[END_REF] and extended for the ANs. Besides methods that approximated the computation of reachability, authors in [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF][START_REF] Prajna | A framework for worst-case and stochastic safety verification using barrier certificates[END_REF] proposed methods for safety verification that do not require the computation of reachable sets, but instead relied on the notion of barrier certificates [START_REF] Prajna | Barrier certificates for nonlinear model validation[END_REF].

A key notion underlying network behavior is that state-space is organized into a number of basins of attraction, connecting states according to their transitions, and summing up the network's global dynamics. For systems having many attractors, some of them may correspond to desired behaviors, while others, to unsafe regions. In the global network dynamic a system can move towards an unsafe region. Therefore, beyond ensuring the properties reachability and safety verification, proposing strategies to interfere on a system and force a desired behavior is important. Previous works introduced the concept of Minimal Intervention Sets (MISs) [START_REF] Klamt | Minimal cut sets in biochemical reaction networks[END_REF] for biological networks and later generalized in [START_REF] Samaga | Computing combinatorial intervention strategies and failure modes in signaling networks[END_REF] for addressing Boolean models of signaling networks. This concept allowed one to choose an intervention strategy to provoke a desired/observed response in certain target nodes. Compared to the approach present of this paper, they do not take into account the transient dynamics of the system, which limits considerably the domain of predictions. [START_REF] Paulevé | Under-approximating cut sets for reachability in large scale automata networks[END_REF]2] proposed approaches based on cut sets to identify nodes/reactions whose perturbation would prevent the occurrence of a given state/reaction. Whereas those prediction can help to control the reachability of an attractor, they do not allow to capture the differentiations processes, as do bifurcations.

Outline Section 2 introduces the Automata Network formalisms on which our method relies. The notion of bifurcation in ANs is introduced in section 3. In section 4, we give a brief overview of methods for checking reachability properties which is a core task for the characterization of bifurcations. In section 5, we combine static analysis, dynamics unfolding, and ASP, aiming at providing a scalable identification of the bifurcation transitions for a given goal. In section 6, we evaluate a prototype implementation on several large biological models in order to support the scalability of our approach. Section 7 concludes the paper.

We consider an Automata Network (AN) as a finite set of finite-state machines having transitions between their local states conditioned by the state of other automata in the network. An AN is defined by a triple (Σ, S, T ) (definition 1) where Σ is the set of automata identifiers; S associates to each automaton a finite set of local states: if a ∈ Σ, S(a) refers to the set of local states of a; and T associates to each automaton the list of its local transitions. Each local state is written of the form a i , where a ∈ Σ is the automaton in which the state belongs to, and i is a unique identifier; therefore given a i , a j ∈ S(a), a i = a j if and only if a i and a j refer to the same local state of the automaton a. For each automaton a ∈ Σ, T (a) refers to the set of transitions of the form t = a i -→ a j with a i , a j ∈ S(a), a i = a j , and the enabling condition of t, formed by a (possibly empty) set of local states of automata different than a and containing at most one local state of each automaton.

Definition 1 (Automata Network (Σ, S, T )). An Automata Network (AN) is defined by a tuple (Σ, S, T ) where -Σ is the finite set of automata identifiers; -For each a ∈ Σ, S(a) = {a i , . . . , a j } is the finite set of local states of automaton a; S We write a i -→ a j ∈ T ∆ ⇔ (a i , , a j ) ∈ T (a).

At any time, each automaton is in one and only one local state, forming the global state of the network. Assuming an arbitrary ordering between automata identifiers, the set of global states of the network is referred to as S as a shortcut for a∈Σ S(a). Given a global state s ∈ S, s(a) is the local state of automaton a in s, i.e., the a-th coordinate of s.

A local transition t = a i -→ a j ∈ T is applicable in a global state s ∈ S when a i and all the local states in are in s. The application of the local transition, noted s • t, replaces the local state of a with a j (definition 2). It results in a (global) transition s t -→ s where s = s • t. In this paper, we consider the asynchronous semantics of ANs: only one local transition can be applied at a time, meaning only one automaton changes its local state by the transitions between two global states. In this semantics, different local transitions may be applicable to a same state, which may potentially lead to very different dynamics. The choice of the transition is non-deterministic. A global state s is reachable from s, noted s → * s , if and only if there exists a (possibly empty) sequence of transitions leading from s to s . Definition 2 (Transition, reachabilility). Given a state s ∈ S and a local transition t = a i -→ a j ∈ T such that s(a) = a i and ∀b k ∈ , s(b) = b j , s • t is the state s where a i has been replaced by a j :

∀b ∈ Σ, (s • t)(b) = a j if b = a s(b) otherwise We then write s t - → s where s = s•t. The reachability binary relation → * ⊆ S×S satisfies s → * s ∆ ⇔ s = s ∨ ∃t ∈ T : s t - → s ∧ s → * s
T (a) = {t 1 = a 1 ∅ -→ a 0 , t 2 = a 0 b0 -→ a 1 , t 3 = a 0 b0,c0 ---→ a 2 } T (b) = {t 4 = b 0 ∅ -→ b 1 , t 5 = b 1 a0 -→ b 0 } T (c) = {t 6 = c 0 a1 -→ c 1 , t 7 = c 1 b1 -→ c 0 , t 8 = c 1 b0 -→ c 2 }
From the given initial state s 0 = a 0 , b 0 , c 0 , 3 transitions can be applied: t 2 , t 3 , and t 4 ; the application of the latter results in s 0 • t 4 = a 0 , b 1 , c 0 (automaton b is now in state b 1 ).

Bifurcations

From an initial state s 0 and a goal local state, we call a bifurcation a transition from a state where the goal is reachable to a state where the goal is not reachable, i.e., there exists no sequence of transition that lead to a state containing the goal local state.

Let us consider the AN of figure 1, with s 0 = a 0 , b 0 , c 0 and the goal a 2 . Figure 2 shows all the possible transitions from s 0 . The states with a gray background are connected to a state containing a 2 (in thick/blue). The transitions in thick/red are bifurcations: once in a white state, there exist no sequence of transitions leading to a 2 . In other words, bifurcations are the transitions from a gray state to a white state. In this example, t 8 is the (unique) local transitions responsible for bifurcations from s 0 to a 2 .

In this paper, given an AN (Σ, S, T ), we are interested in identifying the local transitions t b ∈ T that trigger a bifurcation from a state reached from s 0 ∈ S for a given goal, describing a set of states S g ⊆ S. We call s b a global state where a bifurcation occurs, and s u the global state after the bifurcation:

s u = s b • t b .
The goal is reachable from s b but not from s u . This is illustrated by figure 3. Note that, as illustrated, s b is not inevitably reached: we allow the existence of alternative paths of transitions to the goal. Definition 3 formalizes the notion of bifurcation, where the goal is specified by a local state g 1 (hence S g = {s ∈ S | s(g) = 1}). Note that this goal specification does not loose generality, as one can build an automaton g with local states g 0 and g 1 , and with a local transitions from g 0 to g 1 conditioned by each desired goal state.

a 2 , b 1 , c 0 a 2 , b 0 , c 0 a 1 , b 0 , c 0 a 0 , b 0 , c 1 s 0 = a 0 , b 0 , c 0 a 1 , b 1 , c 0 a 0 , b 1 , c 0 a 1 , b 1 , c 1 a 0 , b 1 , c 1 a 1 , b 0 , c 1 a 1 , b 0 , c 2 a 0 , b 0 , c 2 a 1 , b 1 , c 2 a 0 , b 1 , c 2 t 8 t 8

Definition 3 (Bifurcation).

Given an AN (Σ, S, T ), a global state s 0 ∈ S and a goal local state g 1 with g ∈ Σ and g 1 ∈ S(g), a bifurcation is a transition

s b t b -→ s u of the AN with s b , s u ∈ S and t b ∈ T , such that s 0 → * s b and ∀s ∈ S where s u → * s , s (g) = g 1 .

Reachability and formal approximations

The identification of a bifurcation as presented in the previous section can be decomposed in three steps: we want to find s b ∈ S and t b ∈ T such that [START_REF] Abou-Jaoudé | Model checking to assess t-helper cell plasticity[END_REF]

s b is reachable from s 0 (s 0 → * s b ); (2) g 1 is reachable from s b (∃s g ∈ S : s g (g) = g 1 ∧ s b → * s g ); (2) g 1 is not reachable from s u = s b •t b (∀s ∈ S, s u → * s ⇒ s(g) = g 1 ).
Therefore, alongside the enumeration of candidate s b and t b , reachability is at the core of our bifurcation identification.

In this section, we give a brief overview of the basics of reachability checking, stressing the methods we use in this paper.

State graph and partial order reductions

Given two states s, s of an AN (or an equivalent Petri net), verifying s → * s is a PSPACE-complete problem [START_REF] Cheng | Complexity results for 1-safe nets[END_REF].

The common approach for reachability checking is to build the (finite) set of all the states reachable from s until finding s , by exploring all the possible transitions. However, such a set can be rapidly intractable with large models. Techniques relying on symbolic representations, notably using Binary Decision Diagrams (BDDs) and variants [START_REF] Hamez | Building efficient model checkers using hierarchical set decision diagrams and automatic saturation[END_REF] can improve the scalability of this approach by several orders of magnitude [START_REF] Cimatti | NuSMV 2: An opensource tool for symbolic model checking[END_REF].

In many cases, numerous transitions modelled by ANs are concurrent: their application is independent from each other. For instance, if t 1 and t 2 are concurrent in a state s, one can apply indifferently s • t 1 • t 2 and s • t 2 • t 1 . Such features can be exploited to provide very compact representations of the reachable states in a concurrent system, taking into account the partial order of transition applications. Unfoldings, and more precisely their complete finite prefixes [START_REF] Esparza | Unfoldings -A Partial-Order Approach to Model Checking[END_REF], allow to compute efficiently such compact representations.

In this paper, part of our method uses complete finite prefixes of unfoldings to compute the states that are reachable from the fixed initial state s 0 . Indeed, because biological networks are typically very large, but also very sparse (each node/automaton interacts with a few others, compared to the size of the network), they exhibit a high degree of concurrency for their transitions, making unfolding approaches very effective.

Formal approximations

When facing a large AN, it may turn out that the reachable state space is too large for the aforementioned exact verification of reachability. Moreover, the complexity of the reachability problem can be prohibitive when numerous verifications have to be done, for instance when enumerating candidate initial states (such as s b in our case for checking goal reachability from it).

In this paper, we rely on the reachability approximations for ANs introduced in [START_REF] Paulevé | Static analysis of biological regulatory networks dynamics using abstract interpretation[END_REF][START_REF] Folschette | Sufficient conditions for reachability in automata networks with priorities[END_REF]. We will use both over-approximations (OA) and under-approximations (UA) of the reachability problem: s → * s is true only if OA(s → * s ) is true and s → * s is true if UA(s → * s ) is true; but the converses do not hold in general:

UA(s → * s ) ⇒ s → * s ⇒ OA(s → * s )
The approximations rely on static analysis by abstract interpretation of AN dynamics. We give here the basic explanations for the over-and underapproximation. The analysis rely on the decomposition of the systems dynamics in automata, in order to derive necessary or sufficient conditions for a reachability property of the form s → * s .

The core objects are the local paths within two local states a i , a j of a same automaton a. We call a i a j an objective and define local-paths(a i a j ) the set of the acyclic paths of local transitions between a i and a j . Definition 4 gives the formalization of local-paths where we use the following notations: for a local transition t = a i -→ a j ∈ T , orig(t)

∆ = a i , dest(t) ∆ = a j , enab (t)
∆ = ; ε denotes the empty sequence, and |η| is the length of sequence η.

Definition 4 (local-paths).

Given an objective a i a j ,

-if i = j, local-paths(a i a i ) ∆ = {ε}; -if i = j, a sequence η of transitions in T (a) is in local-paths(a i a j ) if and only if orig(η 1 ) = a i , dest(η |η| ) = a j , ∀n, 1 ≤ n < |η|, dest(η n ) = orig(η n+1 ), and ∀n, m, |η| ≥ n > m ≥ 1, dest(η n ) = orig(η m ). We write t ∈ η ∆ ⇔ ∃n, 1 ≤ n ≤ |η| : η n = t.
Given a local path η, η denotes the union of the conditions of all the local transitions composing it:

η ∆ = |η| n=1 enab(η n ) In the AN of figure 1, local-paths(a 0 a 2 ) = {a 0 b0,c0 ---→ a 2 }; local-paths(c 2 c 1 ) = ∅.
Focusing on the reachability of a single local state g 1 from a state s where s(g) = 0, the analyzes essentially start with the local paths in local-paths(g 0 g 1 ): if g 1 is reachable, then at least of the local path η has to be realizable, meaning that all the local states of its conditions ( η) should be reachable. This lead to a recursive reasoning by repeating the procedure with the objectives from s to the local states in η.

The dependence relationships between the local paths of the different automata can be represented as a graph, where the nodes are all the local states, all the possible objectives, and all their local paths. Such a graph is called a Local Causality Graph (LCG), and abstracts all the executions of the AN.

From a complexity point of view, local paths are computed for each pair of local states within every automata; the length of a local path being at most the number of local states within the automaton, the number of local paths is at most polynomial in the number of local transitions and exponential in the size of the single automaton. In practice, the automata are small, typically between 2 and 4 states for biological models. Therefore, LCGs turn out to be very small compared to the reachable state space of biological networks. They have been successfully applied for analyzing dynamics of ANs with hundreds or thousands of automata, which were intractable with standard model-checking approaches [START_REF] Paulevé | Static analysis of biological regulatory networks dynamics using abstract interpretation[END_REF][START_REF] Paulevé | Under-approximating cut sets for reachability in large scale automata networks[END_REF].

The over-approximation and under-approximation reduce to finding subgraphs of LCGs that satisfy some particular structural properties, which have been proven to be necessary or sufficient for the reachability property, respectively. Whereas the over-approximation can be verified in a time linear with the LCG [START_REF] Paulevé | Static analysis of biological regulatory networks dynamics using abstract interpretation[END_REF], the under-approximation we consider in this paper requires to enumerate over many possible sub-LCGs, but checking if a sub-LCG satisfies the sufficient condition is linear in its size.

Note that further refinements of local-paths have been considered for the mentioned approximations, but for the sake of simplicity, we stick to this coarsegrained presentation in the scope of this paper.

Appendix A gives examples of LCGs satisfying necessary or sufficient conditions for reachability properties in the AN of figure 1.

Identification of bifurcations using ASP

Among the states reachable from s 0 , we want to find a state s b from which (1) the goal is reachable and (2) there exists a transition to a state from which the goal is not reachable. Putting aside the complexity of reachabilities checking, the enumeration of candidate states s b is a clear bottleneck for the identification of bifurcations in an AN.

Our approach combines the formal approximations and (optionally) unfoldings introduced in the previous section with a constraint programming approach to efficiently identify bifurcations. As discussed in the previous section, checking the over-/under-approximations from candidate states and sub-LCGs is easy. For the case of unfolding, checking if a state s belongs to the state space represented by a complete finite prefix is NP-complete [START_REF] Esparza | Unfolding based algorithms for the reachability problem[END_REF]. Therefore, a declarative approach such as Answer-Set Programming (ASP) [START_REF]Knowledge Representation, Reasoning and Declarative Problem Solving[END_REF] is very well suited for specifying admissible s b and t b , and obtaining efficient enumerations of solutions by a solver.

We first present the general scheme of our method, and then given details on its implementation with ASP.

General scheme

A sound and complete characterization of the local transitions t b ∈ T triggering a bifurcation from state s 0 to the goal g 1 would be the following: t b is a bifurcation transition if and only if there exists a state s b ∈ S such that

(C1) s u → * g 1 (C2) s b → * g 1 (C3) s 0 → * s b where s u = s b • t b , s → * g 1 ∆ ⇔ ∀s ∈ S, s → * s ⇒ s (g) = g 1 and s → * g 1 ∆ ⇔ ∃s g ∈ S : s g (g) = g 1 ∧ s → * s g .
However, in an enumeration scheme for s b candidates, checking reachability and non-reachability of the goal from each s b candidate ((C1) and (C2)) is prohibitive. Instead, we relax the above constraints as follows:

(I1 # ) ¬OA(s u → * g 1 ) (I2 # ) UA(s b → * g 1 ) (I3) s b ∈ unf-prefix(s 0 ) (I3 # ) UA(s 0 → * s b )
where unf-prefix(s 0 ) is the set of all reachable states from s 0 represented as the prefix of the unfolding of the AN which has to be pre-computed (section 4.1). Either (I3) or (I3 # ) can be used, at discretion. From OA and UA properties (section 4.2), we directly obtain:

(I1 # ) ⇒ (C1) (I2 # ) ⇒ (C2) (I3) ⇔ (C3) (I3 # ) ⇒ (C3)
Therefore, our characterization is sound (no false positive) but incomplete: some t b might be missed (false negatives). Using (I3) instead of (I3 # ) potentially reduces the false negatives, at the condition that the prefix of the unfolding is tractable. When facing a model too large for the unfolding approach, we should rely on (I3 # ) which is much more scalable but may lead to more false negatives.

Relying on the unfolding from s b (unf-prefix(s b )) is not considered here, as it would require to compute a prefix from each s b candidate, whereas unf-prefix(s 0 ) is computed only once before the bifurcation identification.

ASP syntax and semantics

We give a very brief overview of ASP syntax and semantics that we use in the next section. Please refer to [START_REF] Gelfond | The stable model semantics for logic programming[END_REF][START_REF]Knowledge Representation, Reasoning and Declarative Problem Solving[END_REF][START_REF] Gebser | Answer Set Solving in Practice[END_REF] for an in-depth introduction to ASP.

An ASP program consists in a set of predicates, and logic rules of the form:

1 a0 ← a1, . . ., an, not an+1, . . ., not a n+k .

where a i are atoms, terms or predicates in first order logic. Essentially, such a logical rules states that when all a 1 , . . . , a n are true and all a n+1 , . . . , a n+k are false, then a 0 has to be true as well. a 0 can be ⊥ (or simply omitted): in such a case, the rule is satisfied only if the right hand side of the rule is false (at least one of a 1 , . . . , a n is false or at least one of a n+1 , . . . , a n+k is true). A solution consists in a set of true atoms/terms/predicates with which all the logical rules are satisfied.

ASP allows to use variables (starting with an uppercase) instead of terms/predicates: these pattern declarations will be expanded to the corresponding propositional logic rules prior to the solving. For instance, the following ASP program has as unique (minimal) solution b(1) b( 2) c(1) c(2).

2 c(X) ← b(X). 3 b(1). 4 b(2).
In the following, we also use the notation n { a(X) : b(X) } m which is true when at least n and at most m a(X) are true where X ranges over the true b(X).

ASP implementation

We present here the main rules for implementing the identification of bifurcation transitions with ASP. A significant part of ASP declarations used by (I1 # ), (I2 # ), (I3), and (I3 # ) are generated from the prior computation of local-paths and, in the case of (I3), of the prefix of the unfolding. Applied on figure 1, our implementation correctly uncovers t 8 as a bifurcation for a 2 .

Declaration of local states, transitions, and states Every local state a i ∈ S(a) of each automaton a ∈ Σ are declared with the predicate ls(a,i). We declare the local transitions of the AN and their associated conditions by the predicates tr(id,a,i,j) and trcond(id,b,k), which correspond to the local transition a i {b k }∪ ----→ a j ∈ T . States are declared with the predicate s(ID,A,I) where ID is the state identifier, and A, I, the automaton and local state present in that state. Finally, the goal g 1 is declared with goal(g,1).

For instance, the following instructions declare the automaton a of figure 1 with its local transitions, the state s 0 = a 0 , b 0 , c 0 , and the goal being a 2 : Declaration of s b , t b , and s u The bifurcation transition t b , declared as btr(b), is selected among the declared transitions identifiers (line 6). If a i -→ a j is the selected transition, the global state s u (recall that s (I1 # ) declaration of ¬OA(s u → * g 1 ) This part aims at finding the state s u from which g 1 is not reachable. For that, we designed an ASP implementation of the reachability over-approximation presented in section 4.2. It consists in building a Local Causality Graph (LCG) from pre-computed local-paths. A predicate oa_valid(G,ls(A,I)) is then defined upon G to be true when the local state a i is reachable from the initial state s G . The full implementation is given in appendix B.1.

We instantiate a LCG u with the initial state s u by declaring that the goal is not reachable (oa_valid) (line 11). (I2 # ) declaration of UA(s b → * g 1 ) This part aims at finding the state s b from which g 1 is reachable. Our designed ASP implementation of the reachability under-approximation (section 4.2) consists in finding a sub-LCG G with the satisfying properties for proving the sufficient condition. The edges of this sub-LCG are declared with the predicate ua_lcg(G,Parent,Child). The graph is parameterized by (1) a context which specifies a set of possible initial states and (2) an edge from the node root to the local state(s) for which the simultaneous reachability has to be decided from the supplied context. The full implementation is given in appendix B.2.

We instantiate the under-approximation for building a state s b from which the goal g 1 is reachable: g 1 is a child of the root node of graph b (line 12). The context is subject to the same constraints as s b from s u (lines 13 and 14 reflect lines 9 and 10). Then, s b defines one local state per automaton among the context from which the reachability of g 1 is ensured (line 15), and according to lines 9 and 10.

(I3) declaration of s b ∈ unf-prefix(s 0 ) Given a prefix of an unfolding from s 0 (section 4.1), checking if s b is reachable from s 0 is an NP-complete problem [START_REF] Esparza | Unfolding based algorithms for the reachability problem[END_REF] which can be efficiently encoded in SAT [START_REF] Heljanko | Using logic programs with stable model semantics to solve deadlock and reachability problems for 1-Safe petri nets[END_REF] (and hence in ASP). A synthetic description of the ASP implementation of reachability in unfoldings is given in appendix C. The interested reader should refer to [START_REF] Esparza | Unfoldings -A Partial-Order Approach to Model Checking[END_REF]. Our encoding provides a predicate reach(a,i) which is true if a reachable state contains a i . Declaring s b reachable from s 0 is done simply as follows:

16 reach(A,I) ← s(b,A,I).

(I3 # ) declaration of UA(s 0 → * s b ) An alternative to (I3) which does not require to compute a complete prefix of the unfolding is to rely on the underapproximation of reachability similarly to (I2 # ). The under-approximation is instantiated for the reachability of s b from s 0 with the following statements:

17 ua_lcg(0,root,ls(A,I)) ← s(b,A,I). 18 ctx(0,A,I) ← s(0,A,I).

Experiments

We evaluate our method for the computation of bifurcation transitions in models of actual biological networks showing differentiation capabilities. We have selected three networks showing at least two attractors reachable from a same initial state. In these cases, by supplying an adequate goal state representing one of the attractor, we expect to identify transitions causing a bifurcation from this attractor. We start by introducing our three case studies.

Immunity control in bacteriophage lambda (Lambda phage) . A number of bacterial and viral genes take part in the decision between lysis and lysogenization in temperate bacteriophages. In the lambda case, at least five viral genes (refered to as cI, cro, cII, N and cIII) and several bacterial genes are involved. We apply our method on an automata network equivalent to the model introduced in [START_REF] Thieffry | Dynamical behaviour of biological regulatory networks-ii. immunity control in bacteriophage lambda[END_REF] and with two different goals, corresponding to lysis or lysogenization.

Epidermal Growth Factor & Tumor Necrosis Factor α (EGF/TNF) is a model that combines two important mammalian signaling pathways induced by the Epidermal Growth Factor (EGF) and Tumor Necrosis Factor alpha (TNF α ) [START_REF] Macnamara | State-time spectrum of signal transduction logic models[END_REF]7]. EGF and TNF α ligands stimulate ERK, JNK and p38 MAPK cascades, the PI3K/AKT pathways, and the NFkB cascade. This network encompasses cross-talks between these pathways, as well as two negative feedback loops.

T-helper cell plasticity (Th) has been studied in [START_REF] Abou-Jaoudé | Model checking to assess t-helper cell plasticity[END_REF] in order to investigate switches between attractors subsequent to changes of input conditions. It is a cellular network regulating the differentiation of T-helper (Th) cells, which orchestrate many physiological and pathological immune responses. T-helper (CD4+) lymphocytes play a key role in the regulation of the immune response. By APC activation, native CD4 T cells (Th0) differentiate into specific Th subtypes producing different cytokines which influence the activity of immune effector cell types. Several subtypes (Th1, Th2, Th17, Treg, Tfh, Th9, and Th22) have been well established. We report in this paper four (Pro Th1, Pro Th2, Pro Th17, Pro HTG(Th1 + Th17)) different initial state from which different attractors are reachable. The network is composed of 101 nodes and 221 interactions.

Automata Network Goal (I3) (I3 # ) M-C |pfx| |t b | Time |t b | Time |t b |
Method. For each selected model with initial state and goal, we performed the bifurcation identification following either (I1 # ), (I2 # ), (I3) (unfolding from s 0 ); or (I1 # ), (I2 # ), (I3 # ) (reachability under-approximation). We use clingo 4.5.3 [START_REF] Gebser | Clingo = ASP + control: Preliminary report[END_REF] as ASP solver, and Mole [START_REF] Schwoon | Mole[END_REF] for the computation of the unfolding ((I3)).

The computation times correspond to the total toolchain duration, and includes the local-paths computation, unfolding, ASP program generation, ASP program loading and grounding, and solving. Note that the local-paths computation (and ASP program generation) is almost instantaneous for each case. Source code and models are provided in the supplementary material file.

Results. Table 1 summarizes the results of the experiments. The last experiment shows the limit of the exact analysis of the reachable state space: the computation of the prefix is not tractable on this model. However, the alternative approach (I3 # ) allows to identify bifurcation transitions in this large model. Following section 5.1, (I3 # ) always results in less bifurcations transitions than (I3) with our models. It can be explained with the additional approximation for the reachability of s b from s 0 , using the notations of section 3. One can finally remark that when (I3) is tractable, (I3 # ) shows a slightly slower solving time, albeit of the same order of magnitude. This suggest that checking if a state belongs to an unfolding is more efficient than checking its under-approximation. Finally, because there exists to our knowledge no other method for identifying bifurcations, we cannot compare our results with different methods, and in particular with an exact method in order to appreciate the false negative rate obtained by the (I1 # )-(I2 # )-(I3) scheme.

Conclusion

This paper presents an original combination of computational techniques to identify transitions of a dynamical system that can remove its capability to reach a (set of) states of interest. Our methodology combines static analysis of Automata Networks (ANs) dynamics, partial order representations of the state space, and constraint programming to efficiently enumerate those bifurcations.

To our knowledge, this is the first intregated approach for deriving bifurcation transitions from concurrent models, and ANs in particular. Bifurcations are key features of biological networks, as they model decisive transitions which control the differentiation of the cell: the bifurcations decide the portions of the state space (no longer) reachable in the long-run dynamics. Providing automatic for capturing those differentiations steps is of great interest for biological challenges such as cell reprogramming [START_REF] Crespo | A general strategy for cellular reprogramming: The importance of transcription factor cross-repression[END_REF][START_REF] Abou-Jaoudé | Model checking to assess t-helper cell plasticity[END_REF], as they suggest targets for modulating undergoing cellular processes.

Given an initial state of the AN and a goal state, our method first computes static abstractions of the AN dynamics and (optionnaly) a symbolic representation of the reachable state space with so-called unfoldings. From those prior computations, a set of constraints are issued to identify bifurcation transitions. We used Answer-Set Programming to declare the admissible solutions and the solver clingo to obtain their efficient enumerations. For large models, the unfolding may be intractable: in such a case, the methods relies only on reachability over-and under-approximations. By relying on those relaxations which can be efficiently encoded in ASP, our approach avoids costly exact checking, and is tractable on large models, as supported by the experiments.

Further work will consider the complete identification of bifurcation transitions, by allowing false positives (but no false negatives). In combination with the under-approximation of the bifurcations presented in this paper, it will provide an efficient way to delineate all the transitions that control the reachability of the goal attractor. Future work will also focus on exploiting the identified bifurcations for driving estimations of the probability of reaching the goal at steady state, in the scope of hybrid models of biological networks [START_REF] Shmulevich | Probabilistic boolean networks: a rule-based uncertainty model for gene regulatory networks[END_REF][START_REF] Stoll | Continuous time boolean modeling for biological signaling: application of gillespie algorithm[END_REF]. We propose here a possible encoding of the necessary condition for reachability in ANs outlined in section 4.2 and introduced in [START_REF] Paulevé | Static analysis of biological regulatory networks dynamics using abstract interpretation[END_REF]. Starting from s u (g) = g 0 , the analysis starts with the local paths of the objective g 0 g 1 : g 1 is reachable only if all the conditions of the transitions of at least one local path η ∈ local-paths(g 0 g 1 ) are reachable. This recursive reasoning can be modelled with a graph relating dependencies between objectives, local paths, and local states.

A Examples of Local Causality Graphs

The local paths computed a priori are used to generate the template declaration of the directed edges of the LCG oa_lcg(G,Parent,Child) from each possible objective a i a j . If local-paths(a i a j ) = ∅, the objective a i a j is linked to a node bottom:

1 oa_lcg(G,obj(a,i,j),bottom) ← oa_lcg(G,_,obj(a,i,j)).

otherwise, for local-paths(a i a j ) = {η 1 , . . . , η n }, we declare a node lpath for each different local path m ∈ {1, . . . , n} as a child of a i a j : 2 oa_lcg(G,obj(a,i,j),lpath(obj(a,i,j),m)) ← oa_lcg(G,_,obj(a,i,j)).

then, for each different local state b k ∈ η m in the conditions of the local transitions of η m , we add an edge from the lpath node to ls(b,k):

3 oa_lcg(G,lpath(obj(a,i,j),m),ls(b,k)) ← oa_lcg(G,_,obj(a,i,j)).

In the case when the local path requires no condition ( η m = ∅, this can happen when the objective is trivial, i.e., a i a i , or when the local transitions do not dependent on the other automata), we link the lpath to a node top:

4 oa_lcg(G,lpath(obj(a,i,j),m),top) ← oa_lcg(G,_,obj(a,i,j)).

A LCG G for over-approximation is parameterized with a state s G : if a local path has a local state a j in its transition conditions, the node ls(a,j) is linked, in G, to the node for the objective a i a j (line 5), with a i = s G (a). It is therefore required that state s G defines a (single) local state for each automaton referenced in G (line 6).

5 oa_lcg(G,ls(A,I),obj(A,J,I)) ← oa_lcg(G,_,ls(A,I)), s(G,A,J). 6 1 { s(G,A,J) : ls(A,J) } 1 ← oa_lcg(G, _, ls(A, _)).

The necessary condition for reachability is then declared using the predicate oa_valid(G,N) which is true if the node N satisfies the following condition: it is not bottom (line 7); and, in the case of a local state or objective node, one of its children is oa_valid (lines 8 and 9; or in the case of a local path, either top is its child, or all its children (local states) are oa_valid (lines 10 and 11).

B.2 UA(s → * s ): sufficient condition for reachability

We give here a declarative implementation of the sufficient condition for reachability in ANs outlined in section 4.2 and introduced in [START_REF] Folschette | Sufficient conditions for reachability in automata networks with priorities[END_REF]. The under-approximation consists in building a graph relating objectives, local paths, and local states which satisfies several constraints. If such a graph exists, then the related reachability property is true. Similarly to (I1 # ), we give template declarations for the edges with the predicate ua_lcg(G,Parent,Child). We assume that the reachability property is specified by adding an edge from root to ls(a,i) for each local state to reach.

The graph ua_lcg is parameterized with a context which is a set of local states, declared with the predicate ctx(G,A,J). Every local states a i of the graph that are not part of the reachability specification belong to that context (line 12); and are linked to the objective a j a i for each a j in the context (line 13).

A first constraint is that each objective in the graph is linked to one and only one of its local path. Therefore, objectives without local paths (local-paths(a i a j ) = ∅) cannot be included (line 14), for the others, a choice has to be made among local-paths(a i a j ) = {η 1 , . . . , η n } (line 15).

As for oa_lcg, each local path is linked to all the local states composing its transition conditions: for each m ∈ {1, . . . , n}, for each b k ∈ η m , 16 ua_lcg(G,lpath(obj(a,i,j),m),ls(b,k)) ← ua_lcg(G,_,obj(a,i,j)).

The graph has to be acyclic. This is declared using a predicate conn(G,X,Y) which is true if the node X is connected (there is a directed path) to Y (line 17). A graph is cyclic when conn(G,X,X) (line 18). 17 conn(G,X,Y) ← ua_lcg(G,X,Y). conn(G,X,Y) ← ua_lcg(G,X,Z), conn(G,Z,Y). 18 ← conn(G,X,X).

Then, if the node for an objective a i a j is connected to a local state a k , the under-approximation requires a i a j to be connected with a k a j (assuming that a has at least 3 local states, definition not shown):

19 ua_lcg(G,obj(A,I,J),obj(A,K,J)) ← not boolean(A), conn(G,obj(A,I,J),ls(A,K)).

When a local transition is conditioned by at least two other automata (for instance c • ai,bj ---→ c • ), the under-approximation requests that reaching b j does not involve other local states from a others that a i . This is stated by the indep(G,Y,a,i,ls(b,j)) which cannot be true if b j is connected to a local state a k with k = i line 20. Then, the under-approximation requires that at most one indep predicate is false, for a given LCG G and a given local path Y (line 21). Such an independence should also hold between the local states of the reachability specification (line 22).

C ASP implementation of reachability in unfoldings

A (prefix of an) unfolding is an acyclic bipartite digraph where nodes are either events (application of a transition) or conditions (change of local state) [START_REF] Esparza | Unfoldings -A Partial-Order Approach to Model Checking[END_REF]. We use the predicate post(X,Y) to denote an edge from X to Y ; and h(C,ls(A,I)) to denote that the condition C corresponds to the local state a i . Figure 5 shows an example of unfolding.

A state s belongs to the prefix if it is possible to build a configuration such that all the local states in s have a unique corresponding condition on the cut of the configuration (line 1).

A configuration is a set of events, and we use e(E) to denote that the event E belongs to the configuration. By definition, if E is in a configuration, all its parent events are in the configuration (line 2). There should be no conflicts between two events of a configuration: two events are in conflict if they share a common parent condition (line 3).

A condition is on the cut if its parent event is in the configuration (line 4), and none of its children event is in the configuration (line 5).

∆=

  a∈Σ S(a) is the finite set of global states; LS ∆ = a∈Σ S(a) denotes the set of all the local states. -T = {a → T a | a ∈ Σ}, where ∀a ∈ Σ, T a ⊆ S(a) × 2 LS\S(a) × S(a) with (a i , , a j ) ∈ T a ⇒ a i = a j and ∀b ∈ Σ, | ∩ S(b)| ≤ 1, is the mapping from automata to their finite set of local transitions.

Figure 1

 1 Figure 1 represents an AN (Σ, S, T ) of 3 automata (Σ = {a, b, c}), with S(a) = {a 0 , a 1 , a 2 }, S(b) = {b 0 , b 1 }, S(c) = {c 0 , c 1 , c 2 }, and 8 local transitions defined as follows:

Figure 1 .

 1 Figure 1. An example of Automata Network (AN).' Automata are represented by labelled boxes, and local states by circles where ticks are their identifier within the automaton -for instance, the local state a0 is the circle ticked 0 in the box a. A transition is a directed edge between two local states within the same automaton. It can be labelled with a set of local states of other automata. Grayed local states stand for the global state a0, b0, c0 .

Figure 2 .

 2 Figure 2. Transition graph of the AN in figure 1 from the initial state s0 = a0, b0, c0 . The goal a2 is in thick/blue; the states connected to the goal are in gray; the bifurcations to the goal in thick/red, labelled with the local transitions in the AN definition.

Figure 3 .

 3 Figure 3. General illustration of a bifurcation. s0 is the initial state, Sg is a set of states in which the goal local state is present. The dashed arrows represent a sequence (possibly empty) of transitions. The plain red arrow is a bifurcation from a global state s b to su, and t b is the associated local transition.

  u = s b • t b ) should satisfy s u (a) = a j (line 7) and, ∀b k ∈ , s u (b) = b k (line 8). The state s b should then match s u , except for the automaton a, as s b (a) = a i (lines 9 and 10).

6 1 {

 1 btr(ID) : tr(ID,_,_,_) } 1. 7 s(u,A,J) ← btr(ID),tr(ID,A,_,J). 8 s(u,B,K) ← btr(ID),trcond(ID,B,K). 9 s(b,A,I) ← btr(ID),tr(ID,A,I,_). 10 s(b,B,K) ← s(u,B,K),btr(ID),tr(ID,A,_,_),B!=A.

11 ←

 11 oa_valid(u,ls(G,I)),goal(G,I).

Figure 4

 4 Figure4gives examples of Local Causality Graphs (section 4.2) for approximation reachability of a 2 in the AN of figure1. The left LCG does not satisfy the necessary condition (no local paths from c 2 to c 0 ), hence a 2 is not reachable from the given initial state a 1 , b 0 , c 2 . The middle LCG does satisfy the necessary condition. And the right LCG is a valid sub-LCG for the sufficient condition for a 2 reachability.

a 2 b 0 b 0 c 2 c 0 c 0 a 1 a 2 b 0 c 1 c 0 a 2 b 0 b 1 b 1 b 0 b 1 b 0 b 0 a 0 a 0 b 1 b 1 c 0 c 0 a 0 a 2 rootFigure 4 .

 2020211101024 Figure 4. Exemples of Local Causality Graphs for (left) over-approximation of a2 reachability from a1, b0, c2 (middle) over-approximation of a2 reachability from a1, b0, c1 (right) under-approximation of a2 reachability from a0, b1, c1 . The small circles represent the local paths.

  20 indepfailure(Y,ls(A,I)) ← indep(G,Y,A,I,N), conn(G,N,ls(A,K)), K!=I. 21 ← indepfailure(Y,N),indepfailure(Y,M),M!=N. 22 indep(G,root,A,I,ls(B,J)) ← ua_lcg(G,root,ls(A,I)),ua_lcg(G,root,ls(B,J)),B != A. For η m ∈ local-paths(a i a j ), for each local transition a • -→ a • ∈ η m , for each couple of different local states in its condition b k , c l ∈ , b k = c l :23 indep(G,lpath(obj(a,i,j),m),b,k,ls(c,l)) ← ua_lcg(G,_,lpath(obj(a,i,j),m)).

Figure 5 .

 5 Figure 5. Unfolding of the AN of figure 1. Events are boxed nodes, conditions have no borders and indicate both the automata local state and the condition identifier.

Table 1 .

 1 Experimental results for the identification of bifurcations depending if (I3) or (I3 # ) is used, compared to a exact model-checking (M-C) using NuSMV[START_REF] Cimatti | NuSMV 2: An opensource tool for symbolic model checking[END_REF]. Models Th_th1, Th_th2, Th_th17, Th_HTG are the same automata network but have different initial state. For each model, two different goals have been tested. |Σ| is the number of automata, and |T | the number of transitions; |pfx| is the size (number of events in the partial order structure) of the prefix of the unfolding from the initial state of the model; |t b | is the number of identified bifurcation transitions. Computation times have been obtained on an Intel R Core TM i7-4770 3.40GHz CPU with 16GiB of RAM. OT indicates an out-of-time execution (more than one hour).

	Time

ua_lcg(b,root,ls(G,I)) ← goal(G,I).

ctx(b,A,I) ← btr(ID),tr(ID,A,I,_).

ctx(b,B,K) ← s(u,B,K),btr(ID),tr(ID,A,_,_),B!=A.

1 { s(b,A,I) : ctx(b,A,I) } 1 ← ctx(b,A,_).

ctx(G,A,I) ← ua_lcg(G,N,ls(A,I)), N != root.

ua_lcg(G,ls(A,I),obj(A,J,I)) ← ua_lcg(G,_,ls(A,I)), ctx(G,A,J).

← ua_lcg(G,_,obj(a,i,j)).

1 { ua_lcg(G,obj(a,i,j),lpath(obj(a,i,j),1..n)) } 1 ← ua_lcg(G,_,obj(a,i,j)).