Maxime Folschette
email: maxime.folschette@irccyn.ec-nantes.fr

Loïc Paulevé

Katsumi Inoue

Morgan Magnin

Olivier Roux

Identification of Biological Regulatory Networks from Process Hitting models

Keywords: qualitative modeling, model abstraction, interaction graph, René Thomas's parameters, answer set programming

Qualitative formalisms offer a well-established alternative to the more traditionally used differential equation models of Biological Regulatory Networks (BRNs). These formalisms led to numerous theoretical works and practical tools to understand emerging behaviors. The analysis of the dynamics of very large models is however a rather hard problem, which led us to previously introduce the Process Hitting framework (PH), which is a particular class of nondeterministic asynchronous automata network (or safe Petri nets). Its major advantage lies in the efficiency of several static analyses recently designed to assess dynamical properties, making it possible to tackle very large models.

In this paper, we address the formal identification of qualitative models of BRNs from PH models. First, the inference of the Interaction Graph from a PH model summarizes the signed influences between the components that are effective for the dynamics. Second, we provide the inference of all René-Thomas models of BRNs that are compatible with a given PH. As the PH allows the specification of nondeterministic interactions between components, our inference emphasizes the ability of PH to deal with large BRNs with incomplete knowledge on interactions, where Thomas's approach fails because of the combinatorics of parameters.

The inference of corresponding Thomas models is implemented using Answer Set Programming, which allows in particular an efficient enumeration of (possibly numerous) compatible parametrizations.

Introduction

As regulatory phenomena play a crucial role in biological systems, they need to be studied accurately. Biological Regulatory Networks (BRNs) consist in sets of either positive or negative mutual effects between the components. With the purpose of analyzing these systems, they are often modeled as graphs which make it possible to determine the possible evolutions of all the interacting components of the system. Besides continuous models of physicists, generally designed through systems of ordinary differential equations, modeling regulatory networks by means of Boolean networks has become popular in the wake of Kauffman's work [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF]. We based our work on a logical formalism developed by René Thomas from 1973 [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF] that generalizes upon Boolean networks in the sense that it allows variables to have more than two values and transitions between states to occur asynchronously.

In this approach, the different levels of a component, such as concentration or expression levels, are abstractly represented by (positive) integer values and transitions between these levels may be considered as instantaneous. Hence, qualitative state graphs may be derived from which we are able to formally find out all the possible behaviors expressed as sequences of transitions between these states. Nevertheless, these dynamics can be precisely established only with regard to some discrete parameters, hereafter called "Thomas's parameters", which stand for kinds of "focal points", i.e., the evolutionary tendency from each state and depending on the set of the other currently interacting components.

Thomas's modeling has motivated numerous works around the link between the Interaction Graph (IG) (summarizing the global influences between components) and the possible dynamics (e.g., [START_REF] Richard | Necessary conditions for multistationarity in discrete dynamical systems[END_REF][START_REF] Remy | Graphic requirements for multistability and attractive cycles in a boolean dynamical framework[END_REF]), model reduction (e.g., [START_REF] Naldi | A reduction of logical regulatory graphs preserving essential dynamical properties[END_REF]), formal checking of dynamics (e.g., [START_REF] Richard | Modern Formal Methods and Applications[END_REF][START_REF] Naldi | Decision diagrams for the representation and analysis of logical models of genetic networks[END_REF]), and the incorporation of time (e.g., [START_REF] Siebert | Incorporating time delays into the logical analysis of gene regulatory networks[END_REF][START_REF] Ahmad | Analysing formal models of genetic regulatory networks with delays[END_REF]) and probability (e.g., [START_REF] Twardziok | Stochasticity in reactions: a probabilistic boolean modeling approach[END_REF]) dimensions, to name but a few.

While the formal checking of dynamical properties is often limited to small networks because of the state graph explosion, an other major drawback of this framework is the difficulty to specify Thomas's parameters, especially for large networks. Other works have also been carried out with the aim of dealing with large systems (e.g., [START_REF] Corblin | A declarative constraint-based method for analyzing discrete genetic regulatory networks[END_REF][START_REF] Bollman | Fixed points in discrete models for regulatory genetic networks[END_REF][START_REF] Delaplace | Analysis of modular organisation of interaction networks based on asymptotic dynamics[END_REF][START_REF] Ay | Scalable steady state analysis of boolean biological regulatory networks[END_REF]). Because of the complexity of such systems, these works generally focus on some restricted properties, such as finding steady states.

In order to address the formal analysis of additional dynamical properties (e.g., state reachabilities) within very large BRNs, we recently introduced in [START_REF] Paulevé | Refining dynamics of gene regulatory networks in a stochastic π-calculus framework[END_REF] a new formalism, named the Process Hitting (PH), to model concurrent systems having components with a few qualitative levels. A PH model describes, in an atomic manner, the possible evolutions of a process (representing one component at one level) triggered by the hit of at most one other process in the system. This framework can be seen as a special class of formalisms like Petri Nets or Communicating Finite State Machines, where the arity and the effect of synchronizations between components are restricted. Thanks to the particular structure of interactions within a PH, very efficient static analysis methods have been developed to over-and under-approximate reachability properties making tractable the formal analysis of the qualitative dynamics of BRNs with hundreds or thousands of components, which was impossible with other state-of-the-art 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][START_REF] Folschette | Under-approximation of reachability in multivalued asynchronous networks[END_REF].

The PH is suitable to model BRNs with different levels of abstraction in the specification of cooperations (associated influences) between components. The concept of cooperation refers to the way two (or more) components jointly influence a third one. In other words, it captures the logical functions stating how various elements coalesce and act together upon an other element among the network. In particular, PH allows to specify partial (nondeterministic) cooperations that superpose the dynamics of candidate fully-specified (deterministic) cooperations. This permits to model BRNs with a partial knowledge on precise evolution functions for components by capturing the largest (the most general) dynamics, without a combinatorics enumeration of compatible Thomas models.

The aim of this paper is to formally establish links between several layers suited for specifying the qualitative dynamics of BRNs, namely: the Interaction Graph (IG), referencing the sign of direct regulations between components; the Process Hitting (PH), which permits to specify the qualitative dynamics of the system with various degrees of knowledge for joint regulations; and Thomas's parameters, which completely specify the functions governing the evolution of components. Another motivation for the work depicted in this paper is that it constitutes an important step to go further in the study of large biological regulatory networks systems beyond simple analysis: it makes it possible to partly control some behaviors. Indeed, when one wants to modify some behaviors (e.g., so as to avoid the reachability of some states or to create stable states) we intend to be able to get upstream the results in accordance at the regulatory level. The outcome of the present work is that we could make such changes on the PH model and then recover afterwards the BRN corresponding to the transformed behavioral system. In a more general framework and in an interesting way, we should be able to synthesize families of BRNs having some featuring behavioral properties.

Firstly, we derive rules that, given a PH model, infer the IG corresponding to the dynamics of the encoded BRN. The obtained IG relates only components that have actual influence in the dynamics. This typically results into IGs that contain less interactions than an hypothetical starting prior IG, as the specification of component evolutions may reveal non-functional regulations. Based on the derived IG, various static analysis allows to conclude on global properties of the system dynamics (see [START_REF] Paulevé | Static analysis of boolean networks based on interaction graphs: a survey[END_REF] for a short survey). The most known are the Thomas's conjectures (that have been proved since, e.g., [START_REF] Remy | Graphic requirements for multistability and attractive cycles in a boolean dynamical framework[END_REF][START_REF] Richard | Necessary conditions for multistationarity in discrete dynamical systems[END_REF][START_REF] Richard | Negative circuits and sustained oscillations in asynchronous automata networks[END_REF] for Boolean and discrete networks), which relate the absence of positive cycle to the impossibility for multi-stationarity (distinct attractors), and the absence of negative cycle to the impossibility for oscillatory behaviors.

Second, we tackle the inference of Thomas's parameters that are compatible with a given PH model, i.e., that fully specify the evolution of components while respecting the cooperations specified in PH. More formally, the resulting BRN dynamics is ensured to be included in the PH dynamics: any transition in the (asynchronous) Thomas model exists in the PH model. In practice, the number of compatible Thomas's parameters can be huge for large networks where only a few cooperations are specified. This highlights that the PH is suitable to analyze dynamics of large-scale networks where the knowledge on cooperations is often incomplete: instead of dealing with a possibly non-tractable number of Thomas models, one can already capture precise dynamical features using a single PH model which gathers all the plausible dynamics.

Finally, we discuss on an implementation of these inference schemes using Answer Set Programming (ASP) [START_REF]Knowledge Representation, Reasoning and Declarative Problem Solving[END_REF], which turns out to be effective for these enumerative searches. The whole method is also applied to the model of the bacteriophage lambda immunity control with 4 components, in order to give a detailed application of our method, and to the biological model of the ERBB receptor-regulated G1/S transition, which contains 20 components and allows to show how joint actions can be defined inside a PH model in order to refine its dynamics and study a more precise class of underlying Thomas models.

Our work is related to the approach of [START_REF] Khalis | The SMBioNet method for discovering models of gene regulatory networks[END_REF] which relies on temporal logic, and to [START_REF] Corblin | Applications of a formal approach to decipher discrete genetic networks[END_REF][START_REF] Corblin | Automatic inference of regulatory and dynamical properties from incomplete gene interaction and expression data[END_REF] that use constraint programming. All three of them aim at determining a class of models which are consistent with available partial data on the regulatory structure and dynamical properties. Our method is based on a model rather than on constraints, which allows to define some properties on the system structure (such as cooperations). Furthermore, we claim that we are able to deal with larger biological networks. Finally, it must be noticed that we are not interested in this paper in the derivation of one PH from a BRN (which was previously described in [START_REF] Paulevé | Refining dynamics of gene regulatory networks in a stochastic π-calculus framework[END_REF]) but, on the contrary, in finding out a set of BRNs from one PH.

The contributions presented in this paper significantly extend and improve the preliminary results introduced in [START_REF] Folschette | Concretizing the process hitting into biological regulatory networks[END_REF]. In addition to the improvement of the efficiency and accuracy of the IG inference, we have added support for non-monotonous regulations, that are regulations being positive or negative depending on a particular context. This allows to apply our methodology to a wider class of models. Furthermore, some parts have been completely rewritten in order to explain more precisely the efficiency of the ASP implementation, and give a new detailed biological application of our work.

Outline. Sect. 2 recalls the PH and Thomas frameworks; Sect. 4 defines the IG inference from PH; Sect. 5 details the parameters inference and the enumeration of Thomas's parametrizations compatible with a PH; Sect. 6 discusses the implementation of the latter in ASP. Sect. 7 illustrates our method on a biological model and discusses its applicability on large models. Notations. i; j is the set of integers {i, i + 1, . . . , j}.

Frameworks

The Process Hitting framework

We give the definition and semantics of the Process Hitting (PH), and its usage to model cooperation between concurrent components. Two examples of PH modeling a BRN at different abstraction levels are given. They serve as running examples in the rest of this article.

A PH (Def. 1) gathers a finite number of concurrent processes grouped into a finite set of sorts. A process belongs to a unique sort and is noted a i where a is the sort and i the identifier of the process within the sort a. At any time, one and only one process of each sort is present; a state of the PH thus corresponds to the set of such processes.

The concurrent interactions between processes are defined by a set of actions. Actions describe the replacement of a process by another of the same sort conditioned by the presence of at most one other process in the current state of the PH. An action is denoted by

a i → b j b k where a i , b j , b k are pro- cesses of sorts a and b. It is required that b j = b k and that a = b ⇒ a i = b j . An action h = a i → b j b k is
read as "a i hits b j to make it bounce to b k ", and a i , b j , b k are called respectively hitter, target and bounce of the action, and can be referred to as hitter(h), target(h), bounce(h), respectively. Definition 1 (Process Hitting). A Process Hitting is a triple (Σ, L, H):

• Σ ∆ = {a, b, . . . } is the finite set of sorts; • L ∆ = a∈Σ L a is the set of states with L a = {a 0 , . . . , a la } the finite set of processes of sort a ∈ Σ and l a a positive integer with a = b ⇒ ∀(a i , b j ) ∈ L a × L b , a i = b j ; • H ∆ = {a i → b j b k , . . . | (a, b) ∈ Σ 2 ∧ (a i , b j , b k ) ∈ L a × L b × L b ∧b j = b k ∧ a = b ⇒ a i = b j }
is the finite set of actions.

P denotes the set of all processes (P

∆ = {a i | a ∈ Σ ∧ a i ∈ L a }).
The sort of a process a i is referred to as Σ(a i ) = a and the set of sorts present in an action h ∈ H as Σ(h) = {Σ(hitter(h)), Σ(target(h))}. Given a state s ∈ L, the process of sort a ∈ Σ present in s is denoted by s[a], that is the a-coordinate of the state s. If a i ∈ L a , we define the notation

a i ∈ s ∆ ⇔ s[a] = a i .
Given a set of sorts g ⊂ Σ, L(g) denotes the set of all the sub-states of the sorts in g:

L(g) ∆ = b∈g L b .
The set of sorts having an action on a given sort a is noted hitters(a) (Eq. ( 1)).

∀a ∈ Σ, hitters(a) ∆ = {b ∈ Σ | ∃b i → a j a k ∈ H} (1) An action h = a i → b j b k ∈ H is playable in s ∈ L if and only if s[a] = a i and s[b] = b j .
In such a case, (s • h) stands for the state resulting from the play of the action h in s, that is

(s • h)[b] = b k and ∀c ∈ Σ, c = b, (s • h)[c] = s[c]. For the sake of clarity, ((s • h) • h ), h ∈ H is abbreviated as (s • h • h ).
PH is a particular sub-class of safe Petri nets [START_REF] Bernardinello | A survey of basic net models and modular net classes[END_REF], i.e., Petri nets with at most one token per place, having groups of mutually exclusive places (processes) acting as the sorts, and where each transition has at most one read-arc (connected to the hitter ), and exactly one incoming arc (from the target) and one outgoing arc (to the bounce). The relationship between PH and Petri nets is detailed in [START_REF] Paulevé | From the Process Hitting to Petri Nets and Back[END_REF].

H = {a 2 → b 1 b 0 , b 0 → a 2 a 1 , c 0 → a 2 a 1 , b 0 → a 1 a 0 , c 0 → a 1 a 0 , b 1 → a 0 a 1 , c 1 → a 0 a 1 , b 1 → a 1 a 2 , c 1 → a 1 a 2 } . The action h = b 1 → a 1 a 2 is playable in the state s = b 1 , a 1 , c 0 ; and s • h = b 1 , a 2 , c 0 .
This PH example actually models a BRN where the component a has three qualitative levels and components b and c are Boolean. In this BRN, b and c activate a, while a inhibits b. The inhibition of b by a is only effective when a is at level 2 (noted a 2 ); in the other cases, b cannot evolve in any direction. The activation of a by b (c) is encoded by the actions making the level of a increase (resp. decrease) when b (c) is present (resp. absent). It is worth noticing that the activation of a by b (c) is independent from c (b). This may express a lack of knowledge on the cooperation between these two regulators: we thus model an over-approximation of the possible actions.

Modeling conjunctions. PH restricts the causality of transitions within sorts to the presence of at most one process in another sort (the hitter). One may want to model transitions that should be taken only when several processes are present, i.e., modeling conjunctions -or cooperations -between the presence of several processes. As described in [START_REF] Paulevé | Refining dynamics of gene regulatory networks in a stochastic π-calculus framework[END_REF], such cooperations can be encoded in PH with the use of some intermediate sorts, called cooperative sorts, that merge information on the state of cooperating sorts. Fig. 2 shows an example of cooperation between processes b 1 and c 1 to make a 1 bounce to a 2 : a cooperative sort bc is defined with 4 processes (one for each sub-state of the presence of processes b 1 and c 1 ). For the sake of clarity, the processes of bc are indexed using the substate they represent. Hence, bc 01 represents the sub-state b 0 , c 1 , and so on. Each process of sorts b and c hit bc to make it bounce to the process reflecting the state of the sorts b and c (e.g., b 1 → bc 00 bc 10 and b 1 → bc 01 bc 11 ). Then, it is process bc 11 which hits a 1 to make it bounce to a 2 instead of independent hits from b 1 and c 1 .

We note that cooperative sorts are standard PH sorts and do not involve any special treatment regarding the semantics of related actions. It is also worth noticing that the use of such intermediate sorts may trigger spurious transitions due to potential incoherences between the state of the cooperative sort and the actual state of cooperating sorts. However, such phenomena can be fixed by the use of priorities between actions [START_REF] Folschette | Under-approximation of reachability in multivalued asynchronous networks[END_REF].

When the number of cooperating processes is large, it is possible to chain several cooperative sorts to prevent the combinatoric explosion of the number of processes created within cooperative sorts. For instance, if b 1 , c 1 , and d 1 cooperate, one can create a cooperative sort bc with 4 processes reflecting the presence of b 1 and c 1 , and a cooperative sort bcd with 4 processes reflecting the presence of bc 11 and d 1 . Such constructions are helpful in PH as the static analysis of dynamics developed in [START_REF] Paulevé | Static analysis of biological regulatory networks dynamics using abstract interpretation[END_REF] does not suffer from the number of sorts, but on the number of processes within a single sort.

Example 2. The PH in Fig. 3 results from the refinement of the PH in Fig. 1 where several cooperations have been specified. In particular, the bounce to a 2 is the result of a cooperation between b 1 and c 1 ; and the bounce to a 0 of a cooperation between b 0 and c 0 . Hence, this PH expresses a BRN where a requires both b and c active to reach its highest level, and a does not become inactive unless both b and c are inactive.

Thomas's modeling

The principle of qualitative modeling was introduced as synchronous Boolean networks by Stuart Kauffman on the one hand [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF], and asynchronous René-Thomas networks on the other [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF]. Both models have been considered of interest and led to numerous works. In the following, we choose to focus on Thomas's modeling, with its asynchronous, unitary and multi-level semantics.

In this section, we concisely present Thomas's modeling of a Biological Regulatory Network (BRN) and its dynamics, merely inspired by [START_REF] Thomas | Biological feedback[END_REF][START_REF] Richard | Modern Formal Methods and Applications[END_REF][START_REF] Bernot | Semantics of biological regulatory networks[END_REF]].

Thomas's formalism lies on two complementary descriptions of the system. First, the Interaction Graph (IG) models the structure of the system by defining the components' mutual influences and the conditions of these influences. Then the parametrization removes the ambiguity when a component is targeted by (at least) two different influences. In other words, it specifies the levels towards which a component tends when a given configuration of its regulators applies.

The IG is composed of nodes (a, b, c, . . . ) that represent components, and edges (a s,t -→ b, . . . ) labeled with a sign (s) and a threshold (t) that stand for regulations between these components (Def. 2). The activity, concentration rate or presence of each component in a given state of the system is modeled by an abstract discrete value called expression level. The maximum expression level of a component a is denoted l a . The sign of an edge denotes the kind of regulation it models: it can be positive (+), negative (-) or non-monotonous (±), the latter meaning that the regulation have different sign depending on the context [START_REF] Noual | About non-monotony in boolean automata networks[END_REF]. Regarding the dynamics, an edge a s,t -→ b states that a influences the evolution of b in a certain way when its expression level is above or equal to the threshold t; when its expression level is strictly below this threshold, another kind of influence is expressed, which usually consists of the opposite influence.

The IG is not sufficient to unambiguously define the dynamics of the system. In particular, a component may be regulated by both a positive and negative interaction derived from two different components. To refine the dynamics, a parametrization has been added to the model [START_REF] Snoussi | Qualitative dynamics of piecewise-linear differential equations: a discrete mapping approach[END_REF]. This parametrization allows to define focal points, i.e., it clarifies the level towards which a component evolves depending on the expression of its regulators. Some of the information included in the IG and the parametrization may appear as slightly redundant. Especially, if the single parametrization were enriched with additional constraints (e.g., monotonicity requiring an ordering on the value of parameters), then it may directly capture the signs. However, distinguishing the IG and parametrization allows to decouple the graph (which is generally one of the prior knowledge from the biologists) from the constraints on its dynamics. It opened the way to numerous studies (e.g., necessary condition on sustained oscillations) solely based on the structure of the IG (e.g., the analysis of positive or negative circuits).

Definition 2 (Interaction Graph

). An Interaction Graph (IG) is a couple G = (Γ, E) with Γ the finite set of components, and E the finite set of regulations between two nodes, labeled with a sign and a threshold :

E ∆ = {a s,t -→ b, . . . | a, b ∈ Γ ∧ s ∈ {+, -, ±} ∧ t ∈ 1; l a } ,
where a regulation from a to b is uniquely referenced:

∀a s,t -→ b ∈ E, ∀a s ,t --→ b ∈ E, s = s ∧ t = t .
Given this definition, we denote as a shortcut:

E s ∆ = {a s,t -→ b ∈ E} for s ∈ {+, -, ±}.
Furthermore, for all component b ∈ Γ, we denote E -1 (b) the set of its regulators as defined in Eq. ( 2). Then, for all component a regulating b, i.e., if a s,t -→ b ∈ E, we denote levels(a → b) (resp. levels(a → b)) the interval of expression levels of a above (resp. below) the threshold t (Def. 3). For all expression levels of a that belong to levels(a → b), a is expected to have the influence corresponding to the sign s on b; for all expression levels belonging to levels(a → b), the opposite influence is expected. This definition holds because of the uniqueness of the edge a

E -1 (b) ∆ = {a ∈ Γ | ∃a s,t -→ b ∈ E} (2) a b c +1 +1 -2 +1 +1 K a,∅ = 0 K b,∅ = 0 K a,{b} = 1 K b,{a} = 0 K a,{c} = 1 K b,{b} = 1 K a,{b,c} = 2 K b,{a,b} = 0 K c,∅ = 0 K c,{c} = 1
E + = {b +1 --→ a, c +1 --→ a, b +1 --→ b, c +1 --→ c} E -= {a -2 --→ b} E ± = ∅
Hence:

E -1 (a) = {b, c} E -1 (b) = {a, b} E -1 (c) = {c}
We also have especially:

levels(a → b) = 2; 2 levels(a → b) = 0; 1 A state s of an IG (Γ, E) is an element in S ∆ = a∈Γ 0; l a . s[a]
refers to the level of component a in s. For each possible state, the set of resources of a given component is the set of regulators of this component whose expression level is above the threshold of the regulation (Def. 4). In other words, in every state s,

a regulator b of a component a is either a resource (if s[b] ∈ levels(b → a)) or not a resource (if s[b] ∈ levels(b → a)).
Relying on this observation, the specificity of Thomas's approach lies in the use of discrete parameters to represent the focal level interval towards which the component will evolve in every configuration of the resources (Def. 5).

Definition 4 (Resources (Res)).

For a given component a ∈ Γ and a state s ∈ S, the set of regulators of a whose level in s is above the related threshold is called the set of resources of a in s and is noted Res a (s):

Res a (s) ∆ = {b ∈ E -1 (a) | s[b] ∈ levels(b → a)}
Definition 5 (Parameter K a,ω and Parametrization K). For a given component a ∈ Γ, and ω ⊂ E -1 (a) a set of regulators of a, the parameter K a,ω ∈ 0; l a is a non-negative integer. The complete map K of parameters on G is called a parametrization on G.

An IG and a parametrization make up a complete BRN. Regarding the dynamics, a parameter K a,ω is the value towards which a will tend in the states where its resources are exactly the regulators in ω. Indeed, from a given state s, a transition to another state s is possible provided that exactly one component a evolves of one expression level towards K a,Resa(s) , as stated by the definition of the transition relation s → s (Def. 6). However, a cannot evolve if its expression level already equals the parameter K a,Resa(s) .

Definition 6 (Asynchronous dynamics (→)). The dynamics of a BRN using Thomas's parameters is given by the transition relation → ∈ S × S defined by:

∀s, s ∈ S, s → s ⇐⇒ ∃a ∈ Γ, s[a] = K a,Resa(s) ∧ s [a] = s[a] + δ a (s) ∧ ∀b ∈ Γ, b = a ⇒ s[b] = s [b] with: δ a (s) = +1 if s[a] < K a,Resa(s) -1 if s[a] > K a,Resa (s)
Example 4. Fig. 4(right) gives a Parametrization on the IG of Fig. 4(left). In this BRN, the following transitions are possible given the semantics defined in Def. 6:

a 0 , b 1 , c 1 → a 1 , b 1 , c 1 → a 2 , b 1 , c 1 → a 2 , b 0 , c 1 → a 1 , b 0 , c 1 ,
where a i denotes the component a at level i. This sequence of states ends in a steady state: no evolution is possible in a 1 , b 0 , c 1 .

Remark 1 (Parameters equivalence). In the scope of the asychronous dynamics (Def. 6), one can remark that some parameters values can be equivalent, i.e., the resulting dynamics is the same. It is notably the case with self-regulations: if a s,t -→ a, then for any ω ⊆ E -1 (a), two values v and v of the parameter K a,ω\{a} are equivalent if v ≥ t and v ≥ t. Similary, two values v and v of the parameter K a,ω∪{a} are equivalent if v < t and v < t.

Interaction Graph of the dynamics. From the dynamics specification of a BRN, one can infer back the IG that contain only the functional regulations. Let us define f a (s) = s[a] if s[a] = K a,Resa(s) and f a (s) = s[a] + δ a (s) otherwise. A positive (resp. negative) interaction from b to a is inferred if there exists a state s such that increasing b in s would increase (resp. decrease) f a ; in other words, if there exists a state s so that: f a (s{b i }) < (resp. >) f a (s{b i+1 }), where i < l b and s{b i } denotes the state s where the component b is assigned to i. If all the interactions from b to a are positive (resp. negative), there is a positive (resp. negative) edge from b to a in the IG, denoting a monotonous influence of b on a. In the case when there exists two differently signed interactions, the edge in the IG is said non-monotonous. This abstract representation of the dynamical system can then be used to derive global properties on the dynamics (e.g., [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]).

From Local Transitions to Global Functions

Process Hitting (PH) and Thomas modeling rely on two different paradigms for the specification of network dynamics: Thomas's parametrization defines a global function which associates to each component and each resource configuration a level towards which the component eventually evolves in the given configuration. Such a function is always deterministic (whereas the associated asynchronous dynamics can be nondeterministic due to concurrent evolution of components). In contrast, PH models specify a list of local transitions between the levels (processes) of the components (sorts) that are conditioned with the presence of other processes. In addition, PH allows the use of intermediate sorts, such as cooperative sorts, that do not refer to the components of the network.

Hence, identifying Thomas models from PH requires to lump a set of transitions that can be applied in a given configuration into a single process that correspond to the level towards which the component evolves. This step is detailed in Subsect. 3.1 in which is defined the focals function associating the farthest reachable processes, called focal processes, of a given sort for a given configuration. It is possible that several different focal processes are identified, indicating an nondeterministic evolution of the sort; and it is also possible that no focal process exist, when in presence of terminal cycles of transitions. In those cases, there is no possible Thomas specification for the PH. The Subsect. 3.2 discusses on constraints for obtaining deterministic (local) behaviors in PH.

Finally, Subsect. 3.3 addresses the separation of PH sorts into components and intermediate cooperative sorts, the latter having to be masked in Thomas models. It notably defines the predecessors, the (groups of) regulators of a sort, and the extension of a configuration to incorporate the state of the cooperative sorts.

Focal Processes

Given a sort a ∈ Σ and a configuration delimited by a sub-set of processes of sort a and a sub-state σ of the sorts having an action on a, focals(a, S a , σ) (Def. 7) is the set of processes of sort a towards which a will eventually converge in the scope of the configuration. The configuration delimits the set of states s ∈ L where s[a] ∈ S a , and for all sort b ∈ hitters(a

), b = a, s[b] = σ[b].
focals(a, S a , σ) is obtained from the digraph where a i ∈ S a is connected to a j ∈ L a only if there exists an action b k → a i a j where either b k = a i or b k is in σ. focals(a, S a , σ) is empty if there is any cyclic terminal connected components; otherwise, it is exactly the leafs of the digraph.

Definition 7 (focals(a, S a , σ)). Given a sort a ∈ Σ, a sub-set of its processes S a ⊂ L a and a sub-state σ ∈ b∈hitters(a),b =a L b , focals(a, S a , σ)

∆ = ∅ if tscc(V, E) = ∅ {a i ∈ V | (a i , a j ) ∈ E} otherwise, where E ∆ = {(a i , a j ) ∈ (S a × L a ) | ∃b k → a i a j ∈ H : (b k = a i ∨ σ[b] = k)} (3)
V ∆ = S a ∪ {a j ∈ L a | ∃(a i , a j ) ∈ E} (4)
and tscc(V, E) are the non-elementary terminal strongly connected components of the digraph (V, E):

tscc(V, E) = {W ∈ scc(V, E) | #W ≥ 2 ∧ ∀a i ∈ W, (a i , a j ) ∈ E ⇒ a j ∈ W }
with scc(V, E) the strongly connected components of the digraph (V, E).

From Def. 7 can be derived Property 1 -if focals(a, S a , σ) is empty, there exists a limit cycle in the evolution of a in the given configuration; and Property 2 -if focals(a, S a , σ) is non empty, all the evolutions of a in the scope of the given configuration eventually terminate, and the resulting processes are in focals(a, S a , σ). In other words, if focals(a, S a , σ) is empty, there exists a sequence of actions that may be played an unbound number of times (cycle) in the given configuration; if it is non-empty, it is ensured that any state in the configuration converges in a bounded number of steps to a process in focals(a, S a , σ). , either h ∈ H with target(h) ∈ S a playable in s and s[a] ∈ focals(a, S a , σ); or there exists a sequence of actions h 1 , . . . , h n in H sequentially playable in s with target(h 1 ) ∈ S a and ∀m ∈ 1; n -1 , bounce(h m ) = target(h m+1 ) where bounce(h n ) ∈ focals(a, S a , σ) and where

h n+1 ∈ H with target(h n+1 ) = bounce(h n ) that is playable in s • h 1 • • • • • h n .
Example 5. In the PH of Fig. 1, we obtain:

focals(a, L a , b 0 , c 0 ) = {a 0 } focals(a, L a , b 1 , c 1 ) = {a 2 } focals(a, L a , b 1 , c 0 ) = ∅ focals(a, {a 1 }, b 1 , c 0 ) = {a 0 , a 2 }

Restrictions for Determinism

As discussed at the beginning of this section, and as formally described in the above sub-section, in the scope of a given configuration, a PH sort may converge to one among several different processes (indeterminism) or may never converge (infinite loop).

We call a deterministic sort (Def. 8) a sort a that has a single focal process for each possible configuration of its predecessors in the scope of L a .

Definition 8 (Deterministic sort). A sort a ∈ Σ is a deterministic sort if and only if each configuration σ of its predecessors leads a to a unique focal process, denoted a(σ):

∀σ ∈ L(hitters(υ)), focals(a, L a , σ) = {a(σ)}
For the inference of the interaction graph and Thomas's parameters from a PH model, the following sections assume that all the cooperative sorts are deterministic.

Separating components from cooperative sorts

The identification of a BRN from a PH assumes that the PH defines two types of sorts: the sorts corresponding to BRN components -noted Γ -, and the cooperative (intermediate) sorts -noted ∆ -which should not appear in the BRN. In this subsection, we first give a criteria for identifying the sorts of a PH as either a component or cooperative sort. Then, we characterize the well-formed PH for IG inference.

The delimitation of sorts modeling components relies on the observation that their processes represent (ordered) qualitative levels. Hence an action on such a sort cannot make it bounce to a process at a distance more than one. This set of sorts in denoted by Γ (Eq. ( 5)), whereas the set of cooperative sorts is denoted by ∆ (Eq. ( 6)).

Γ ∆ = {a ∈ Σ | b i → a j a k ∈ H, |j -k| > 1} (5)
∆ ∆ = Σ \ Γ (6)
Given a PH and a partition of its sorts in components Γ and cooperative sorts ∆, Property 3 establishes conditions for BRN identification: in addition of having Γ compatible with Γ and cooperative sorts being deterministic (Def. 8), we also require that there is no cycle between cooperative sorts, and that sorts being never hit (i.e., serving as an invariant environment) are components.

Property 3 (Well-formed Process Hitting for BRN identification). A PH is well-formed for BRN identification with Σ = Γ ∪ ∆ only if the following conditions hold:

1. Γ ∩ ∆ = ∅, Γ ⊆ Γ and ∀υ ∈ ∆, υ is a deterministic sort (Def. 8); 2. there is no cycle between cooperative sorts (the digraph

(Σ, {(a, b) ∈ (Σ × Σ) | ∃a i → b j b k ∈ H ∧ a = b ∧ {a, b} ∩ Γ = ∅}) is acyclic); 3. sorts having no action hitting them belong to Γ ({a ∈ Σ | b i → a j a k ∈ H} ⊂ Γ).
Example 6. In the PH of Fig. 3, bc is a deterministic sort as defined in Def. 8:

focals(bc, L bc , b 0 , c 0 ) = {bc 00 } focals(bc, L bc , b 0 , c 1 ) = {bc 01 } focals(bc, L bc , b 1 , c 0 ) = {bc 10 } focals(bc, L bc , b 1 , c 1 ) = {bc 11 }
Hence, both Fig. 1 and Fig. 3 are well-formed PH for BRN identification with Γ = {a, b, c} and ∆ = {bc}.

Assuming a PH and a split Γ and ∆ of its sort that satisfy Property 3, Def. 9 characterizes the set of predecessors of a sort a as the sorts influencing a through direct actions or intermediate cooperative sorts. The predecessors of a that are components are the regulators of a, denoted reg(a) (Eq. ( 9)).

Definition 9 (pred(a) and reg(a)). Given a ∈ Σ, pred(a) ⊆ Σ is the smallest set satisfying the following conditions:

hitters(a) ⊆ pred(a) (7)
υ ∈ pred(a) ∩ ∆ ⇒ pred(υ) ⊆ pred(a) (8)
The components that are predecessors of a are referred to as reg(a):

reg(a) ∆ = pred(a) ∩ Γ (9)
As described in the introduction of this section, the identification of Thomas models from PH relies on the inference of focal processes in the different configurations of the regulators of each component. In order to reduce the scope of the configurations that need to be enumerated, we introduce the notion of group of regulators of sort a as a set of components that have a joint influence on a. More precisely, b and c are in the same group of regulators of a if there exists an intermediate cooperative sort υ that hits a such that b and c are regulators of υ. Those groups form a partition of reg(a). We denote X(a) the finest partition in groups of regulators (Def. 10), i.e., where each group is minimal. If the cooperative sorts are all deterministic (Def. 8), it is sufficient to specify the configuration σ of a group g of regulators of a sort a to obtain the configuration of all the sorts hitting a directly. Indeed, because of the absence of cycles between cooperative sorts (Property 3), one can recursively evaluate the focal process of each cooperative sort that hits a with respect to the configuration σ of the regulators g. Such an extended configuration is denoted by ς g a (σ), formalized in Def. 11

Definition 11 (Configuration extension ς g a (σ)). Given a ∈ Σ, a sub-set of regulators g ⊆ reg(a) and a configuration σ ∈ L(g), ς g a (σ) is the configuration σ with the corresponding processes of the cooperative sorts in hitters(a):

ς g a (σ) ∆ = σ υ(σ) | υ ∈ hitters(a) ∩ ∆ ∧ reg(υ) ⊆ g (10)
υ(σ

) ∆ = focals(υ, L υ , σ ) if hitters(υ) ⊆ dom(σ ) υ(σ υ (σ ) | υ ∈ hitters(υ) ∩ ∆ ) otherwise, (11)
where dom(σ ) is the set of sorts defined in σ and denotes the union of two (sub)states.

Example 7. In the PH of Fig. 3,

ς {b,c} a ( b 0 , c 1 ) = b 0 , c 1 , bc 01 and ς {b,c} a ( b 1 , c 1 ) = b 0 , c 1 , bc 11 .

Interaction Graph Inference from Process Hitting

The Interaction Graph (IG) is an abstract representation of the direct qualitative influences, positive and/or negative, between the components of the system. As discussed in Sect. 1, the IG allows to efficiently characterize global dynamical properties for the concrete system, such as the capability for multistationarity or oscillation.

The inference of the IG also allows to check for the consistency of the model with respect to prior knowledge on the influences. This is the case for Thomas models, where a prior IG is required to specify the parameters of the dynamics. As detailed in [START_REF] Richard | Modern Formal Methods and Applications[END_REF], the IG of a dynamical model is consistent with a prior IG if it is a sub-graph of the prior IG. On the one hand, having influences in the model not referenced in the prior IG indicates either mistakes in the encoding of the dynamics, or potential influences that might need to be verified in the biological system. One the other hand, some influences specified in the prior IG may not have been necessary to encode the desired dynamics, indicating some non-minimality of the prior IG.

We consider hereafter a global PH (Σ, L, H) and a split of sorts Γ ∪ ∆ = Σ satisfying Property 3 on which the IG inference is to be performed. The inference of the IG is described in Subsect. 4.1, and is illustrated on small examples in Subsect. 4.2.

Inference of Influences

We aim at inferring that b activates (inhibits) a if there exists a configuration where increasing the level of b makes possible the increase (decrease) of the level of a, following the standard IG inference from Boolean and discrete networks [START_REF] Richard | Negative circuits and sustained oscillations in asynchronous automata networks[END_REF]. This can be seen as looking for changes of derivatives in the dynamics of a when the configuration changes.

Proposition 1 details the inference of all existing influences between components occurring with a threshold t. The inference relies on finding configurations σ of groups g of regulators of a such that the increase of one of the regulators b ∈ g from b t to b t+1 changes the direction of the evolution of a. If b = a (Eq. ( 12)), this is achieved by finding a configuration σ of g and a such that there exists a bounce a j from a i in σ{b t } (i.e., the state σ where the process of sort b has been replaced with b t ) which is different from a bounce a k from a i in σ{b t+1 }. If so, b regulates a at threshold t + 1 with the sign of k -j. The case where b = a (Eq. ( 13)) is similar: we look for a t such that there exists a bounce in opposite direction in σ{a t } than in σ{a t+1 }. If there is no bounce in one of these two configurations, we also need to ensure that the configuration corresponds to a local fixed point between all the regulators of a (instead of in scope of g). The set of configurations where no bounce occur on a i is given by Φ(a i ) (Eq. ( 16)).

Proposition 1 (Influences inference). We define the set of positive (resp. negative) influences Ê+ (resp. Ê-) by

∀a ∈ Γ, ∀b ∈ reg(a), b = a,∀s ∈ {+, -}, ∀t < l b b t+1 --→ a ∈ Ês ∆ ⇔∃g ∈ X(a), b ∈ g, ∃σ ∈ L(g ∪ {a}), ∃a j ∈ B g a (σ{b t }), ∃a k ∈ B g a (σ{b t+1 }), k > j ∧ s = + ∨ k < j ∧ s = - (12) ∀a ∈ Γ,∀s ∈ {+, -}, ∀t < l a a t+1 --→ a ∈ Ês ∆ ⇔∃g ∈ X(a), ∃σ ∈ L(g ∪ {a}), ∃a j ∈ B g∪{a} a (σ{a t }), ∃a k ∈ B g∪{a} a (σ{a t+1 }), (j = t ⇒ ∃σ ∈ Φ(a t ) : σ{a t } ⊆ σ ) ∧ (k = t + 1 ⇒ ∃σ ∈ Φ(a t+1 ) : σ{a t+1 } ⊆ σ ) ∧ (k ≥ t + 1 ∧ j ≤ t ∧ s = + ∨ k < t + 1 ∧ j > t ∧ s = -) (13)
where B g a (σ) ⊆ L a is defined as follows:

Bg a (σ) ∆ = {a j | b k → a i a j ∈ H, b ∈ g ∧ b k ∈ ς g a (σ) ∧ a i ∈ σ} (14) B g a (σ) ∆ = Bg a (σ) if Bg a (σ) = ∅ {σ[a]} otherwise ( 15
)
and Φ(a i ) ⊆ L(reg(a) ∪ {a}) is defined as follows:

Φ(a i ) ∆ = {σ ∈ L(reg(a) ∪ {a}) | a i ∈ σ ∧ focals(a, {a i }, ς reg(a) a (σ)) = {a i }} (16)
We are now able to infer the edges of the final IG by considering positive and negative influences (Proposition 2). We infer a positive (resp. negative) edge if there only exist corresponding influences with the same sign. If an influence is both positive and negative, we infer a non-monotonous edge. In the end, the threshold of each edge is the minimum threshold for which an influence has been found. Note that as regulations match with sign changes in the evolution of the regulated components, if two positive (negative) regulations from b to a have been inferred, necessarily, a third regulation of opposite sign has been inferred in between. The only exception is when a is never hit (reg(a) = ∅); in such a case, we arbitrarily pick the minimum threshold (that is, 1) for the self-activation.

Proposition 2 (Interaction Graph inference). We infer G = (Γ, E) using Proposition 1 as follows:

E -= {a -,t --→ b ∈ Ê-| a t -→ b ∈ Ê+ } E + = {a +,t --→ b ∈ Ê+ | a t -→ b ∈ Ê-∧ t = min{l | a l - → b ∈ Ê+ }} E ± = {a ±,t --→ b | ∃a t -→ b ∈ Ê+ ∧ ∃a t -→ b ∈ Ê- ∧ t = min{l | a l - → b ∈ Ê-∪ Ê+ }}

Examples

The IG inference on the PH of Fig. 3 returns the IG in Fig. 4(left): it finds back all the expected edges alongside with their signs and thresholds. Furthermore, the IG inference on the PH of Fig. 1, which is not refined, gives the same result, despite the absence of the cooperative sort. This is due to the fact that the inference of Proposition 1 is able to analyze the influences from different regulators independently if they do not share a cooperative sort.

If an action a 2 → b 0 b 1 is added to the PH of Fig. 3, then two nonmonotonous edges towards b are inferred instead of the previous signed edges:

E + = {b +1 --→ a, c +1 --→ a, a +1 --→ a, c +1 --→ c} E -= ∅ E ± = {a ±2 --→ b, b ±1 --→ b}
This is due to the fact that the actions a 2 → b 1 b 0 and a 2 → b 0 b 1 introduce an oscillation only caused by a, which cannot be represented in Thomas's modeling.

Parametrization inference

Given a PH and an IG (either arbitrary, or inferred following the previous section), this section addresses the identification of the parametrization for the corresponding Thomas model.

As described in Subsect. 5.1, this identification relies on the computation of the focal processes for each configuration of the parameters. When there is a unique focal process of a component for a given configuration of its regulators, the focal process matches with the value of the associated Thomas's parameter.

However, in the general case, as described in Sect. 3 and even when the PH is well-formed for parameter inference (characterized in the next subsection by Property 4), there may be several focal processes (nondeterministic behavior) or none (cyclic behavior). Such a case occurs notably when the PH encodes the union of several Boolean or discrete networks, as described in [START_REF] Paulevé | Refining dynamics of gene regulatory networks in a stochastic π-calculus framework[END_REF][START_REF] Paulevé | Logical Modeling of Biological Systems[END_REF].

Hence, we propose in Subsect. 5.2 the notion of compatible parametrization with respect ot the PH dynamics: a Thomas model is compatible with a PH if its dynamics is included in the PH dynamics, i.e., all the transitions in the Thomas model are possible transitions in the PH. This relaxed notion of parameterization compatibility allows to enumerate all parametrizations compatible with a given PH, i.e., all Thomas model whose dynamics is included in the PH dynamics.

Parameters inference

This subsection addresses the inference of independent discrete parameters from a given PH. The inference is equivalent to the one in [START_REF] Paulevé | Refining dynamics of gene regulatory networks in a stochastic π-calculus framework[END_REF]. In addition, we characterize the well-formed PH for parameter inference property (Property 4), which implies that any process in levels(b → a) (resp. levels(b → a)) share the same behavior regarding a.

Property 4 (Well-formed PH for parameter inference). A PH is wellformed for parameter inference if and only if it is well-formed for BRN identification (Property 3) and that the associated IG (Γ, E) verifies that:

∀b → a ∈ E, b = a, ∀σ ∈ L(reg(a)), ∀(i, j ∈ levels(b → a) ∨ i, j ∈ levels(b → a)), ∀a k ∈ L a , focals(a, {a k }, ς reg(a) a (σ{b i })) = focals(a, {a k }, ς reg(a) a (σ{b j })) (17) ∀a ∈ Γ,∀b ∈ reg(a), b = a ∧ b → a / ∈ E, ∀σ ∈ L(reg(a)), ∀b i , b j ∈ L b , ∀a k ∈ L a , focals(a, {a k }, ς reg(a) a (σ{b i })) = focals(a, {a k }, ς reg(a) a (σ{b j })) (18)
Let K a,ω be a Thomas's parameter for a given component a ∈ Γ with ω ⊂ E -1 (a) a set of its regulators. As described in Subsect. 2.2, K a,ω specifies to which values a eventually evolves in the configuration matching with ω.

The configuration of the PH corresponding to ω is given as follows. For each component b ∈ E -1 (a), we define σ b a,ω (Eq. ( 19)) as the process of b with the level in levels

(b → a) if b ∈ ω, or in levels(b → a) if b / ∈ ω. Because of Property 4,
if several possible levels exist, this process can be chosen arbitrarily. If b does not regulate a in the IG, all processes of sort b should have the same action on a, so the process b 0 is arbitrarily selected (Eq. ( 20)). The configuration σ a,ω corresponding to ω (Eq. ( 21)) is then obtained by extending the configuration of the regulators of a to the (deterministic) cooperative sorts (Def. 11 in Subsect. 3.3). Finally, we denote by S a (σ) (Eq. ( 22)) the set of processes of sort a that are compatible with a configuration σ: if a is specified in σ, then S a (σ) = {a i } where a i ∈ σ; otherwise S a (σ) is the set of all processes of sort a, that is, L a .

∀b → a ∈ E, σ b a,ω ∆ = min(levels(b → a)) if b ∈ ω, min(levels(b → a)) if b / ∈ ω (19) ∀b ∈ reg(a), b → a / ∈ E, σ b a,ω ∆ = b 0 (20) σ a,ω ∆ = ς reg(a) a ( σ b a,ω | b ∈ reg(a) ) (21) ∀a ∈ Γ, S a (σ) ∆ = {a i } if a i ∈ σ L a otherwise. ( 22
)
Therefore, we obtain that K a,ω = focals(a, S a (σ a,ω ), σ a,ω ) if this latter is a singleton, denoting a deterministic behavior (Proposition 3).

Proposition 3 (Parameter inference). Let (Σ, L, H) be a Process Hitting with Σ = Γ ∪ ∆ and an associated IG G = (Γ, E) well-formed for parameter inference. For all a ∈ Γ, for any ω ⊆ E -1 (a), if focals(a, S a (σ a,ω ), σ a,ω ) = {a i }, then K a,ω = i.

Example 8. When applied to the refined PH of Fig. 3, and using the IG of Fig. 4(left) which is the result of the IG inference on this PH model (as explained in Subsect. 4.2), Proposition 3 infers the parametrization of Fig. 4(right). The inference of parameters thus finds back all parameters of the original model when the cooperations are fully defined.

Example 9. Regarding the non-refined PH of Fig. 1, and using the same IG of Fig. 4(left), the parameters of Fig. 4(right) are inferred, except the parameters K a,{b} and K a,{c} for which Proposition 3 is not conclusive. The obtained parametrization is therefore partial. This is due to the lack of precision in the cooperation between b and c, caused by the absence of the cooperative sort.

Given the Proposition 3, we see that in some cases, the inference of the targeted parameter is impossible. This can be due to a lack of cooperation between regulators: when two regulators independently hit a component, their actions can have opposite effects, leading to two possible evolutions. Such an indeterminism is not possible in Thomas modeling as in a given configuration of regulators, a component can only have an interval attractor, and can therefore evolve in only one direction. In order to avoid such inconclusive cases, one has to ensure that no such behavior is allowed by either removing undesired actions or using cooperative sorts to prevent opposite influences between concurrent regulators.

Admissible parametrizations

When building a BRN, one has to find the parametrization that best describes the desired behavior of the studied system. Complexity is inherent to this process as the number of possible parametrizations for a given IG is exponential w.r.t. the number of components. However, the method of parameters inference presented in this section gives some information about necessary parameters given a certain dynamics described by a PH. This information thus drops the number of possible parametrizations, allowing to find the desired behavior more easily.

We first delimit the validity of a parameter (Property 5) in order to ensure that any transition in the resulting BRN is allowed by the studied PH. This is verified by the existence of a hit making the concerned component bounce into the direction of the value of the parameter in the matching context. Thus, assuming Property 4 holds, any transition in the inferred BRN corresponds to at least one transition in the PH, proving the correctness of our inference. Any parameter that does not satisfy Property 5 is therefore excluded from the enumeration. We remark that all parameters inferred by Proposition 3 satisfy this property.

Property 5 (Parameter validity).

A parameter K a,ω is valid w.r.t. the PH if and only if the following equation is verified:

∀a i ∈ C a a,ω , a i = K a,ω =⇒ (∃c k → a i a j ∈ H, c k ∈ C c a,ω ∧ a i < K a,ω ⇒ j > i ∧ a i > K a,ω ⇒ j < i)
Then, we use some additional biological constraints on Thomas's parameters given in [START_REF] Bernot | Semantics of biological regulatory networks[END_REF], that we sum up in the following three properties: Property 6 (Extreme values assumption). Let G = (Γ, E) be an IG. A parametrization K on G satisfies the extreme values assumption if and only if:

∀b ∈ Γ, E -1 (b) = ∅ =⇒ ∃ω ⊂ E -1 (b), K b,ω = 0 ∧ ∃ω ⊂ E -1 (b), K b,ω = l b
Property 7 (Activity assumption). Let G = (Γ, E) be an IG. A parametrization K on G satisfies the activity assumption if and only if:

∀b ∈ Γ, ∀a ∈ E -1 (b), ∃ω ⊂ E -1 (b), K b,ω = K b,ω∪{a}
Property 8 (Monotonicity assumption). Let G = (Γ, E) be an IG. A parametrization K on G satisfies the monotonicity assumption if and only if:

∀b ∈ Γ, ∀A + ⊂ {a ∈ Γ | a +,t --→ b ∈ E + }, ∀A -⊂ {a ∈ Γ | a -,t --→ b ∈ E -}, K b,ω∪A -≤ K b,ω∪A +
Example 10. The parametrization inferred in Example 9 was partial because K a,{b} and K a,{c} could not be inferred. It is however possible to enumerate all admissible parametrizations compatible with both the inferred parameters, and the properties of this subsection. This enumeration gives 9 different parametrizations which correspond to the 3 possible values for both K a,{b} and K a,{c} :

K a,{b} ∈ {0, 1, 2} and K a,{c} ∈ {0, 1, 2}
We note that all these solutions satisfy Properties 5 to 8. All the parametrizations obtained from these combinations are thus admissible. Finally, we note that the value 1 belongs to the possible values for both parameters. Therefore this enumeration allows, from the model in Fig. 1, to find the behavior of the model refined with a cooperative sort described in Fig. 3.

Answer Set Programming implementation concepts

Answer Set Programming (ASP) is a logic programming paradigm [START_REF]Knowledge Representation, Reasoning and Declarative Problem Solving[END_REF][START_REF] Baral | Using answer set programming to model multi-agent scenarios involving agents' knowledge about other's knowledge[END_REF], which has been chosen to address the enumeration of all admissible parametrizations. ASP can efficiently enumerate a large set of possible answers, and it is easy to constrain the answers according to some properties. Applied to the models and constraints defined in this paper, ASP appears to be efficient for tackling the inherent complexity of the enumeration of parametrizations.

Hereafter, we synthesize some key points of our ASP implementation with the enumeration example.

Simple rules and answer sets

ASP is based on a set of rules of the form:

H head ← A 1 , A 2 , . . . , A n , ¬B 1 , ¬B 2 , . . . , ¬B m body .
where the body is a series of atoms (A i ) and negations of atoms (¬B i ). In the case of simple rules (as opposed to the cardinality rules of Subsect. 6.3), the head is also an atom (H). Such a rule, whose formal semantics is defined below, schematically states that if all atoms A 1 , A 2 , . . . , A n are true and all atoms B 1 , B 2 , . . . , B m are not true (negation by failure), then H has to be true.

An atom is composed of a predicate and a series of arguments (possibly empty). For example, the following atom:

p(x 1 , x 2 , . . . , x r )
is composed of the predicate p and r arguments: x 1 , x 2 , . . . , x r . Each argument is either a constant, which is a representation of a piece of data (component name, expression level, . . . ), or a variable, which is in fact a shorthand for any possible constants (variables are detailed below). In this paper, constants are either numerical or consist of a single lowercase letter (e.g., a, b, c, 1, 2, . . . ) while variables are always denoted by a single capital letter (e.g., A, P , Q, . . . ). We do not use function symbols as arguments in this work.

Example 11. Consider a PH model such as depicted in Fig. 1: in order to state the existence of each component, we use a predicate called component with two arguments:

component(x, n).
where x is the name of the component and l x = n is its maximum expression level.

An ASP program is a set of rules as described above. Solving an ASP program means finding an answer set, which is a minimal set of atoms that respect all the rules. In order to formally define this notion of answer set, let us define a definite rule as a rule with no negation of atoms (noted "¬" above), and a definite program as a program containing only definite rules. Let S be a set of atoms: a definite rule is satisfied by S if H ∈ S, or if ∃i ∈ 1; n , A i / ∈ S. Given this definition of satisfaction, we define the answer set of a definite program Π as the (unique) minimal set S of atoms that satisfies all the rules in Π.

In order to consider the general case, for all non-definite program Π and set S of atoms, we denote Π S the reduct of Π w.r.t. S, defined from Π by: 1. deleting all the rules that have a negation of atom ¬B i in the body where B i ∈ S, and 2. removing all negations of atoms in the bodies of the remaining rules.

Then, a set of atoms S is an answer set of a non-definite program Π if it is the answer set of Π S , which is a definite program. We note that several answer sets can be solution to the same non-definite program, and in practice a solver can be directed to enumerate them all.

Note that it is possible to define a simple rule with no body part. Such a rule is called a fact, and its head atom consequently has to belong to all answer sets. For instance, the information describing the studied model (the original PH model and the inferred IG and parameters) are expressed in ASP using facts.

Example 12. In order to define the 3 components of Fig. 1, we use the following program:

component(a, 2). component(b, 1).

component(c, 1).

This program contains only facts using the predicate component defined in Example 11, and a, b, c, 1 and 2 are constants.

Variables

To describe the sets of all expression levels of each component (i.e., the set 0; l a for each a ∈ Γ), one can use atoms of the form component_levels(a, k) to state that k ∈ 0; l a . Variables here come in handy to enumerate each possible constant k for each component a: before solving, any rule containing variables is grounded, that is, replaced by an equivalent set of rules with constants only. The following rule, for example, contains three variables (A, K and M ) and enumerates the set of possible expression levels of each component in the system:

component_levels(A, K) ← component(A, M ), 0 ≤ K ≤ M.
where the notation "≤" stands for a shortcut in ASP which has the same meaning as the mathematical operator.

Example 13. Regarding Fig. 1, the previous rule together with the facts of Example 12 will give the following answer set:

{ component(a, 2), component (b, 1), component(c, 1), component_levels(a, 0), component_levels(a, 1), component_levels(a, 2), component_levels(b, 0), component_levels(b, 1), component_levels(c, 0), component_levels(c, 1). } 6.3

. Cardinality rules

As an extension of simple rules, cardinality rules turn out to be convenient to enumerate a set of answer sets. The head of a cardinality rule specifies a set of atoms H and two integers min and max, and is denoted:

min { H } max ← A 1 , A 2 , . . . , A n , ¬B 1 , ¬B 2 , . . . , ¬B m .
Given such a rule, as many answer sets as possible are created, so that each answer set S verifies: min ≤ |S ∩ H| ≤ max and every atom H i ∈ S ∩ H respects the simple rule:

H i ← A 1 , A 2 , . . . , A n , ¬B 1 , ¬B 2 , . . . , ¬B m .
In other words, all answer sets contain a subset of H whose cardinality goes from min to max, and for which the condition in the body of the cardinality rule is met. The set of atoms H = {H 1 , H 2 , . . . , H p } is often defined as: H = {P | Q}, which is a shorthand for "the set of atoms of the form P for which Q is true".

Cardinality rules turn out to be convenient to enumerate all possible parametrizations by creating multiple answer sets. For functional purposes, a unique label is assigned to every possible set of resources of a given component. Thus, we denote ω p the set of resources of a given component a labeled by p, and naturally, K a,ωp is the related parameter. We note that labeling the sets of resources of a component is obviously equivalent to labeling its parameters. Then, suppose that:

• param_label (a, p) states that p is a valid label for a set of resources of component a (and therefore K a,ωp is a valid parameter);

• param(a, p, i) states that: K a,ωp = i;

• infered _param(a, p) states that the parameter inference of K a,ωp was conclusive (Proposition 3).

It is thereby possible to enumerate the possible values of all parameters for which Proposition 3 was not conclusive, with the following cardinality rule:

1 { param(A, P, I) | component_levels(A, I) } 1 ←
param_label (A, P ), ¬infered _param(A, P ).

Indeed, this rule applies to any possible parameter P of any component A (param_label ) whose value is still unknown (¬infered _param), and states that any expression level I of this component (component_levels) is a candidate value for the parameter (param). Furthermore, the lower and upper bounds are both 1, which forces each enumerated parameter to have exactly one value. In other worlds, this cardinality rule creates as many answer sets as there are candidate parametrizations so that if K a,ωp could not be inferred by Proposition 3, then K a,ωp ∈ 0; l a (thus completely disregarding the notion of admissible parametrizations given in Subsect. 5.2).

Example 14. In the scope of Example 9, K a,{b} and K a,{c} could not be inferred by Proposition 3. The previous cardinality rule allows to produce 9 parametrizations, in which these two parameters can take all possible values:

(K a,{b} ; K a,{c} ) ∈ {0, 1, 2} × {0, 1, 2}
and all the other parameters keep their inferred values.

Constraints

Finally, a constraint is a rule with no head part:

← A 1 , A 2 , . . . , A n , ¬B 1 , ¬B 2 , . . . , ¬B m .
A constraint is satisfied only if its body is not satisfied, which thus allows to invalidate answer sets containing some unwanted combinations of atoms. In the scope of parameters enumeration, for example, constraints are especially useful to filter out parametrizations that do not respect the assumptions of Subsect. 5.2. Indeed, suppose that:

• less_active(a, p, q) states that ω p is a set of resources of a with (loosely) less activators and more inhibitors than ω q ;

• param_inf (a, p, q) states that: K a,ωp ≤ K a,ωq .

Then, the monotonicity assumption (Property 8) is formulated as the following constraint:

← less_active(A, P, Q), ¬param_inf (A, P, Q).

Indeed, this constraint removes all parametrization results where parameters K A,ω P and K A,ω Q exist such that A is less activated by the set of resources ω P than it is by ω Q , but K A,ω Q < K A,ω P , thus violating the monotonicity assumption. Of course, other assumptions can be formulated in the same way.

Example 15. All the candidate values enumerated in Example 14 already respect Properties 5 to 8. Therefore, the constraints encoding these properties do not filter any solution.

This subsection succinctly described how ASP programs come in handy to represent a model and solve complex problems on it. It finds a particularly interesting application in the enumeration of parameters: all possible parametrizations are generated in separate answer sets, and integrity constraints are formulated to remove those that do not fit the assumptions of admissible parametrizations, thus reducing the number of candidate parametrizations to be considered in the end. However, all steps of the inference presented in this paper (Sect. 4 & 5) were implemented in and benefited from this programming paradigm.

Examples

This section aims at giving several applications of our work in order to understand its results and range of applications. First, Subsect. 7.1 gives a detailed application of our IG inference method by studying parts of a model of the phage lambda immunity control. Then, Subsect. 7.2 gives a practical application of our results on a biological model of epithelial growth factor receptor taken from the literature. Finally, Subsect. 7.3 gives data related to the results and computation times of our software when applied to several large biological models.

The inference and enumeration methods described in this paper have been implemented as part of Pint1 , which gathers PH related tools. Our implementation mainly consists in ASP programs that are solved using Clingo2 . The IG and parameter inferences can be performed using the following command:

ph2thomas -i model.ph --dot ig.dot where model.ph is the PH model file in Pint format, and ig.dot is an output file to write the inferred IG in DOT format. The (possibly partial) inferred parametrization will be returned on the standard output. Instead of the --dot ig.dot option, it is possible to specify an alternative IG to perform the parameters inference on (instead of using the inferred IG) with the --ig model.ig option, where model.ig is an IG model file in Pint format. The admissible parametrizations enumeration is performed by appending the --enumerate option to the command.

Detailed example: the bacteriophage lambda immunity control

In order to illustrate the IG inference developed as Sect. 4, this first subsection focuses on a model of the lambda phage immunity control, which has been widely studied mainly for its particular switch mechanism, allowing this virus to "chose" between lysis and lysogenization. These two possible responses are characterized by different dynamics that lead to two separate attractors.

We consider here the model with four genes proposed in [START_REF] Thieffry | Dynamical behaviour of biological regulatory networks -II. immunity control in bacteriophage lambda[END_REF], which will serve as example to detail the inference method proposed in the previous sections. This model can be straightforwardly represented under the form of a Process Hitting with 4 components and 4 cooperative sorts:

Σ = {cI, cro, cII, N Γ , cI-cro-cII, cI-cro-N, cI-cro, cro-cII ∆ }
The four components have respectively 3, 4, 2 and 2 expression levels, which we denote:

L cI = {cI 0 , cI 1 , cI 2 } L cro = {cro 0 , cro 1 , cro 2 , cro 3 } L cII = {cII 0 , cII 1 } L N = {N 0 , N 1 }
The full Process Hitting model is not represented here due to a lack of space, but excerpts are given in Fig. 5 and Fig. 6. The whole set of actions is given in appendix.

Regulation of cI on N

We first study the influence of cI on N. The relevant part of the Process Hitting model for this problem is depicted in Fig. 5. It involves a third component cro and the cooperative sort cI-cro; indeed, by studying the full model, we find:

pred(N) = {cI, cro, cI-cro} reg(N) = {cI, cro} X(N) = {{cI, cro}}
Thus, we only have to focus on the group of regulators {cI, cro} in order to analyze the regulation cI → N. The idea behind the IG inference developed in Sect. 4 is to make the level of the considered regulator vary (here, cI) while keeping constant all the other regulators in the same group (here, only cro) and to observe the results on the evolution of the observed component (here, N). In other words, if increasing the sole level of cI makes N tend to decrease, then a negative local influence of cI on N is inferred. Let us now focus instead on the group of regulators {cI}. One again, we trivially have:

cI 0 1 2 cro 0 1 2 3 N 0 1 cI-cro F F cI < 1 cro < 2 F T cI < 1 cro ≥ 2 T F cI ≥ 1 cro < 2 T T cI ≥ 1 cro ≥ 2
∀cI i ∈ L cI , ∀cro j ∈ L cro , ς {cI} cro ( cI i , cro j ) = cI i , cro j
Furthermore, from Eq. ( 15), we especially find:

B {cI} cro ( cI 0 , cro 2 ) = {cro 3 } B {cI} cro ( cI 0 , cro 3 ) = {cro 2 }
Therefore, Eq. ( 12) gives: cro 3 -→ cro ∈ Ê-. Finally, as no more information is added by observing the other levels of cro and cI, and given Proposition 2, it comes: cro -3

--→ cro ∈ E.

Conclusion on the phage lambda immunity control

The IG inference of Proposition 2 applied to the whole Process Hitting modeling the phage lambda immunity response produces the IG given in Fig. 7. This IG contains all edges included in the original model except the self-influence on cI (namely, cI +2 --→ cI / ∈ E); this is due to the fact that an equivalent dynamics exists without this edge.

Finally, we note that the parameters inference of Proposition 3 applied to the same model, and using the output of the previous IG inference is conclusive for all parameters. The complete parametrization produced is available in Table 1.

We note two main differences between the parametrization obtained by our inference and the one proposed in [START_REF] Thieffry | Dynamical behaviour of biological regulatory networks -II. immunity control in bacteriophage lambda[END_REF]: • The absence of the self-regulation of cI automatically reduces the size of the parametrization, but this does not modify the dynamics, as mentioned previously;

cI cro cII N -2 -2 -1 -1 -3 -3 -2 +1 +1
• Our inference returns K cro,{cI;cro} = 2, which is a value equivalent to the original value 0 according to Remark 1: the resulting dynamics is the same, and 1 would also have been an acceptable value.

Biological application: the epithelial growth factor receptor

This subsection focuses on the study of the epidermal growth factor (EGF) receptor model detailed in [START_REF] Sahin | Modeling ERBB receptor-regulated G1/S transition to find novel targets for de novo trastuzumab resistance[END_REF]. This model is represented by an IG containing 20 components and 52 edges. A protein named EGF, having no regulator, can be considered as the only input, and a chain of reactions leads to the activation of protein pRB which is responsible for regulating the cell division, therefore making it essential to prevent cancer development.

Three models are created from the original IG, with different levels of precision regarding the cooperation between components. The translation from an IG into a Process Hitting is not detailed here as it was previously covered in [START_REF] Paulevé | Refining dynamics of gene regulatory networks in a stochastic π-calculus framework[END_REF]. Model (1) represents a translation of the raw IG into Process Hitting, that is, without any knowledge of the Boolean rules (and therefore the cooperations) of the components. Model (2) implements some of the rules based on the results of several knockdown experiments. Model (3) is the totally refined model with all cooperations implemented given the Boolean functions of all components. Therefore, those three models can be considered as successive refinements of the original and most general one. The results of the IG and parameters inferences are detailed in Table 2 and discussed in the following alongside with details about their construction.

Model [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF] encompasses only sole interactions between components, that is, independent activations or inhibitions of a component on another given the regulations specified in the original IG. Therefore, the IG inferred from model [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF] Table 2: Results of the IG and parameters inference on three models derived from the EGF receptor model of [START_REF] Sahin | Modeling ERBB receptor-regulated G1/S transition to find novel targets for de novo trastuzumab resistance[END_REF] with different precisions in the definition of the cooperations. Model [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF] contains no cooperations between the components. Several cooperations were included in model ( 2) under the form of 14 cooperative sorts and all of them were included in Model (3) under the form of 22 cooperative sorts. The second column gives the number of edges in the IG inferred with Proposition 2 (the number of nodes is always the number of components in the model, that is, 20). The third column gives the number of parameters in the model (given the IG), the fourth column gives the number of parameters that could be inferred using Proposition 3, and the fifth column consequently gives the number of compatible models with the studied PH model, which exponentially depends on the number of parameters that could not be inferred. Finally, the last column gives the number of fixed points in the PH model, computed with another existing PH tool provided with Pint.

functions could in fact be simplified in a way that a component did not appear anymore in it. No edge have therefore been inferred by our method in this case. Finally, model (3) was build using all the Boolean functions provided in [START_REF] Sahin | Modeling ERBB receptor-regulated G1/S transition to find novel targets for de novo trastuzumab resistance[END_REF]. These functions take the form of 22 cooperative sorts into the model in order to match the desired behavior of the system. As all cooperations are fully defined in this model, all the parameters are inferred and only one Thomas model can be derived. We note also that this PH model is the only one containing at least one fixed point. In fact, the three found steady states include the two states that correspond to a complete propagation of the input signal, that is, in the case where EGF is active and in the case where it is not. The two other models contain no fixed point because some cooperations are not fully defined, leading to oscillations that are a consequence of the nondeterministic behavior.

Computation times on several large models

The current implementation can successfully handle large PH models of BRNs found in the literature such as:

• the EGF receptor model from [START_REF] Sahin | Modeling ERBB receptor-regulated G1/S transition to find novel targets for de novo trastuzumab resistance[END_REF] with 20 components presented in the last subsection3 ,

• a T cell receptor model described as an IG in [START_REF] Klamt | A methodology for the structural and functional analysis of signaling and regulatory networks[END_REF], which contains 40 components and 14 cooperative sorts.

For each model, IG and parameters inferences are performed together in less than a second on a standard desktop computer. Bigger models related to the aforementioned systems were also tested with our implementation:

• a model of the T cell receptor with 94 components, described in [START_REF] Saez-Rodriguez | A logical model provides insights into t cell receptor signaling[END_REF],

• a model of the EGF receptor with 104 components, described in [START_REF] Samaga | The logic of egfr/erbb signaling: Theoretical properties and analysis of high-throughput data[END_REF].

These two models were obtained in a previous work by an automatic translation from the CellNetAnalyzer [START_REF] Klamt | Structural and functional analysis of cellular networks with cellnetanalyzer[END_REF] formalism.

The composition of all models and the results of the inferences are summarized in Table 3. In every case, all parameters could be inferred and it was not necessary to enumerate compatible parametrizations.

Model

5s 27496

Table 3: Computation times and several pieces of information related to the IG and parametrization inferences of four biological models. The second column gives the number of components of each model and the third column gives the number of cooperative sorts used to model joint actions. The fourth (resp. fifth) column gives the computation times of the IG inference (resp. the parametrization inference). The last column gives the number of parameters in each model.

Conclusion and Discussion

This work establishes the abstraction relationship between PH and Thomas's approaches for qualitative BRN modeling. The PH allows an abstract representation of BRNs dynamics (allowing incomplete knowledge on the cooperation between components) that cannot be exactly represented in René Thomas's formalism by a single instance of BRN parametrization. This motivates the concretization of PH models into a set of compatible Thomas models in order to benefit of the complementary advantages of these two formal frameworks.

We first propose an original inference of the Interaction Graph (IG) from a BRN having its dynamics specified in the PH framework. An IG gives a compact abstract representation of the influence of the components between each others. Then, based on a prior inference of René Thomas's parametrization for BRNs from a PH model, we delimit the set of admissible Thomas's parametrizations that are compatible with the PH dynamics, and give arguments for their correctness. A parametrization is compatible with the PH if its dynamics (in terms of possible transitions) is included in the PH dynamics. The enumeration of such parametrizations is efficiently tackled using Answer Set Programming. We illustrate the overall method with several results on both small and large biological models.

Several extensions of the presented work are now to be considered. First, the link between successively refined models of a system could be formally studied. Indeed, it is convenient to refine a PH model by removing actions or adding cooperations; the study and formalization of such process would allow to predict behavioral changes and lead to more accurate models. Second, the inference of BRN multiplexes [START_REF] Bernot | Gene regulatory networks with multiplexes[END_REF] may be of practical interest as they allow to implicitly reduce the possible parametrizations by making cooperations appear in the IG. Because of its atomicity, the PH allows to specify a range of cooperations that cannot be completely captured by a single instance of BRN multiplexes, then encouraging the inference of a set of compatible ones. Finally, in order to improve the performances in the IG inference, we will consider projection operations on the PH structure to undo cooperations between components and reduce the cardinality of configurations to explore by making the interactions independent.

Figure 1 :

 1 Figure1: A Process Hitting (PH) example. Sorts are represented by labeled boxes, and processes by circles (ticks are the identifiers of the processes within the sort, for instance, a 0 is the process ticked 0 in the box a). An action (for instance b 1 → a 1 a 2 ) is represented by a pair of directed arcs, having the hit part (b 1 to a 1 ) in plain line and the bounce part (a 1 to a 2 ) in dotted line. Here, actions involving b 1 or c 1 are in thick lines.

Example 1 .

 1 Fig. 1 represents a PH (Σ, L, H) with Σ = {a, b, c}, L a = {a 0 , a 1 , a 2 }, L b = {b 0 , b 1 }, L c = {c 0 , c 1 }, and

Figure 2 :Figure 3 :

 23 Figure 2: A PH modeling a cooperativity between b 1 and c 1 to make a 1 bounce to a 2 . Actions involving b 1 or c 1 are in thick lines.

Figure 4 :

 4 Figure 4: (left) IG example. Components are represented by nodes labeled with a name and regulations by edges labeled with their sign and threshold. For instance, the edge from b to a is labeled +1, which stands for: b +1 --→ a. This means that if the expression level of b is equal to (i.e., above) 1, then b activates a; otherwise, b inhibits a. (right) Example parametrization on the left IG.

Definition 3 (

 3 Effective levels (levels)). If a s,t -→ b ∈ E, we define: levels(a → b) ∆ = t; l a and levels(a → b) ∆ = 0; t -1 . Example 3. Fig. 4(left) represents an Interaction Graph (Γ, E) where Γ = {a, b, c}, with l a = 2 and l b = l c = 1, and:

Property 1 .

 1 focals(a, S a , σ) = ∅ if and only if there exists a state s ∈ L where s[a] ∈ S a and ∀b ∈ hitters(a), b = a, s[b] = σ[b], such that ∀n ∈ N there exists a sequence of actions h 1 , . . . , h n+1 in H sequentially playable in s with target(h 1 ) ∈ S a and ∀m ∈ 1; n -1 , bounce(h m ) = target(h m+1 ).

Property 2 .

 2 If focals(a, S a , σ) = ∅, for all state s ∈ L where s[a] ∈ S a and ∀b ∈ hitters(a), b = a, s[b] = σ[b]

Definition 10 (

 10 Partition X(a) of reg(a)). Given a ∈ Σ, X(a) is the finest partition of reg(a) such that for any b, c ∈ reg(a), b = c, if there exists a cooperative sort υ ∈ hitters(a) ∩ ∆ such that {b, c} ⊆ reg(υ) then there exists g ∈ X(a) such that {b, c} ⊆ g. If reg(a) = ∅, X(a) ∆ = {∅}.

Figure 5 :Figure 6 :

 56 Figure5: A portion of the PH model used in this section modeling the phage lambda immunity response. This part of the whole model only depicts the interactions towards N from components cI and cro. All the actions updating the cooperative cI-cro are simply depicted by two zigzag arrows. Furthermore, as this cooperative sort cI-cro is deterministic, its processes are labeled by the configuration they represent, given by two Boolean variables indicating whether the values of cI and cro lie above their respective thresholds (T ) or not (F ). For example, the configuration cI 1 , cro 0 is represented by the process T F .

Figure 7 :

 7 Figure 7: Result of the IG inference performed on the Process Hitting model of the phage lambda immunity response taken from [34].

  is the same as the IG used to create the model, with one additional positive

	Model |E| |K|	Inferred parameters	Possible models	Fixed points
	(1) 52	196 20	2 176 9.6 • 10 52 0
	(2) 51	192 98	2 94 2.0 • 10 28	0
	(3) 51	192 192	1	3

  |Γ| |∆| IG inference K inference

					|K|
	EGF receptor [35] 20	22	<1s	<1s	192
	T cell receptor [36] 40	14	<1s	<1s	143
	T cell receptor [37] 94	39	100s	<1s	578
	EGF receptor [38] 104 89	200s	2.

Available at http://loicpauleve.name/pint

Available at http://potassco.sourceforge.net

All models mentioned in this section are available as examples distributed with Pint.

Acknowledgement. This work was partially supported by the Fondation Centrale Initiatives and the French National Agency for Research (ANR-10-BLANC-0218 BioTempo project).

We first note that in the case where cro 3 is present, we have from Eq. [START_REF] Twardziok | Stochasticity in reactions: a probabilistic boolean modeling approach[END_REF]:

We can then deduce the following directions of evolution from Eq. ( 15):

Therefore, Eq. ( 12) gives: cI 1 -→ N ∈ Ê-. The configurations using the other processes of cro (cro 0 , cro 1 and cro 2 ) do not add more information on the influences towards N. Thus, given Proposition 2, it comes: cI

We note that an equivalent work on cro also gives the second influence on

Self-regulation of cro due to cI

We now wish to study the self-regulation of cro. For this, we will show that studying only the actions self-hitting cro is not enough; indeed, it is required to also study the interactions with cI. We have the following results from the whole model:

The interactions of cro and cI on cro have been represented in Fig. 6. Let us first focus on the group of regulators {cro}. Of course:

and from Eq. ( 15) and ( 16):

We thus cannot infer any local influence from Eq. ( 13) by observing the group of regulators {cro}, because Φ(cro 1 ) = Φ(cro 2 ) = ∅. This is due to the fact that all processes of cI can hit cro 1 and cro 2 ; these processes can therefore not be considered stable in any configuration. auto-edge on the only input EGF (which is due to its absence of regulators). The only parameters that could be inferred are the parameters for the extreme cases of regulation (all activators present and all inhibitors absent, and the opposite). This first model therefore abstracts a large number of Thomas models as a lot of parameters are left undecided.

In order to build model (2), 14 cooperative sorts were added in order to model the Boolean functions of several components (consisting of AND and OR operands). To do so, the following components were noticed due to their importance in the chain of reactions: CDK4, CDK6, CycD1, ERα and c-MYC. Indeed, knockdown experiments have been conducted in [START_REF] Sahin | Modeling ERBB receptor-regulated G1/S transition to find novel targets for de novo trastuzumab resistance[END_REF] and the results showed that knocking down these components lead to an important decrease in the production of pRB. We therefore concluded that these components were involved in other components' Boolean functions in a way that the knockdown of the former was sufficient to prevent the production of the latter (which is typical of AND operands). In order to reproduce such requirements, the Boolean functions of their successors, that are CDK4, CDK6, prB, p21, p27, IGF1R, MYC, CycD1 and CycE1, were modeled as cooperative sorts, if needed. In theory, 9 cooperative sorts would have sufficed, but the chaining of cooperative sorts described in Subsect. 2.1 was used to reduce the number of processes in each cooperative sort. As a result, the added cooperations allowed to infer about half the parameters; however, the number of possible Thomas models that can be inferred from this PH is still significant because of the numerous remaining unknown parameters. Furthermore, we note that the inferred IG contains one edge less than the original IG. This is due to the fact that one of the Boolean Appendix. List of actions of the phage lambda immunity control model used in Subsect. 7.1

N0 → cI-cro-N1 cI-cro-N0 N0 → cI-cro-N3 cI-cro-N2 N1 → cI-cro-N0 cI-cro-N1 N1 → cI-cro-N2 cI-cro-N3 }