Adrien Richard

Jean-Paul Comet

Gilles Bernot

Formal Methods for Modeling Biological Regulatory Networks

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Biological systems are one of the most fascinating aspects in biology. They control such diverse dynamics phenomena as temperature control in warm-blooded animals; differentiation of a zygote into the various specialized organs, tissues and cells of the mature organism; the fate of certain viruses, called temperate bacteriophages, which upon infection of a bacterial population can behave in two extremely different ways. Most of these infected cells display a response called lytic: virus multiplies and kills cells. But, a fraction of the cells become lysogenic bacteria and carry the viral genome in a dormant form making the host immune towards infection of other virus.

Computational systems biology tries to establish methods and techniques that enable us to understand such systems as systems, including their robustness, design and manipulation. It means to understand : the structures and the dynamics of systems, methods to control, design and modify systems to cope with desired properties. The modelling contributes in a major way to reach these aims by introducing methods for understanding, simulating and predicting the behaviour of the systems. However, the modelling of biological systems is currently subject to two major difficulties: the biochemical reaction mechanisms underlying the interactions of systems are usually not or incompletely known and quantitative information on kinetic parameters or molecular concentration is rarely available. Thus the modelling activity needs an interaction with the experimental biology to confront models to biological objects. Consequently as in the design of large computing systems, two activities can be distinguished in the modelling step:

1. build a rigorous model of the system satisfying the assumed behaviour corresponding to biological knowledge, 2. design experiments to verify a posteriori the model predictions.

Here we would like to show that some methods from computer science can be reused in the context of system biology, as, for example, formal methods for validation and verification used for the design of large computing systems. For designing experiments, we just mention that the test methods via test generation from model theories may be an efficient way to propose experiments permitting biologists to validate or refute models. For building a rigorous model, the model checking verification tool is particularly suited. In this chapter we present an application of this formal method to build qualitative models of biological regulatory networks.

A biological regulatory network describes interactions between the biological entities, often macromolecules or genes, of a given system. It is statically represented by an interaction graph whose vertices abstract biological entities and arcs their interactions. For describing the evolution of the system, the concentration level of each entity is represented by a value associated to the corresponding vertex. The temporal evolution of these levels constitutes the dynamics of the system.

Ordinary differential equation systems have been first used for describing the dynamics of networks. They are powerful tools particularly to model metabolic processes [START_REF] Voit | Computational Analysis of biochemical systems: a practical guide for biochemists and molecular biologists[END_REF]. However, due to the non-linearity of biological regulations, these differential equation systems cannot often be solved analytically. They can be solved numerically to any desired precision, but this precision itself may be misleading because the values of the parameters and the shape of interactions often have to be guessed for lack of experimental data. This remark led Thomas to simplify the models: he introduced in the 70's a Boolean approach to capture the qualitative nature of the dynamics and he proved its usefulness in the context of immunity in bacteriophages [START_REF] Thomas | A complex control circuit. Regulation of immunity in temperate bacteriophages[END_REF][START_REF] Thomas | Logical analysis of systems comprising feedback loops[END_REF]. Later on, he generalized his formalism to multi-valued levels of concentration (the so called multi-valued logic or ``generalized logical approach'' [START_REF] Thomas | Dynamical behaviour of biological regulatory networks -I. Biological role of feedback loops an practical use of the concept of the loop-characteristic state[END_REF]) since the Boolean idealization may be too caricatured to correctly model biological systems. It has been proved that this qualitative description allows the representation of the essential qualitative features of an ordinary differential equation system provided that the differential equations are piece-wise linear [START_REF] Snoussi | Qualitative dynamics of a piecewise-linear differential equations : a discrete mapping approach[END_REF]. The underlying parameters of the qualitative description can be deduced from the kinetic parameters of the continuous system but can take only a finite number of values. Consequently, all possible qualitative features of the system can be reduced to a finite number of models i.e. parameterisations.

Certainly the most important concepts of the generalized logical analysis are those of positive and negative feedback circuits. If an entity tends to favour (resp. decrease) its own production via the feedback circuit, the circuit is said positive (resp. negative). It has been conjectured by Thomas [START_REF] Thomas | On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations[END_REF] and then proved in different contexts [START_REF] Plathe | Feedback loops, stability and multistationarity in dynamical systems[END_REF][START_REF] Snoussi | Necessary conditions for multistationarity and stable periodicity[END_REF][START_REF] Cinquin | Roles of positive and negative feedback in biological systems[END_REF][START_REF] Cinquin | Positive and negative feedback: striking a balance between necessary antagonists[END_REF][START_REF] Soulé | Graphical requirements for multistationarity[END_REF] that at least one positive circuit is necessary to generate multi-stationarity whereas at least one negative circuit is necessary to obtain a stable oscillatory behaviour. These concepts are especially important since when modelling biological systems, differentiation and homeostasis have often to be taken into consideration. In such cases, these biological constraints can reduce drastically the set of models to consider. These properties can be reinforced by introducing some more complex properties on the dynamics of the system extracted from the biological knowledge or hypotheses. It becomes necessary to construct models which are coherent not only with the previous conditions of multistationarity/homeostatis but also with the additional ones. Formal methods from computer science should be able to help modeller to automatically perform this verification [START_REF] Pérès | Contribution of computational tree logic to biological regulatory networks: example from pseudomonas aeruginosa[END_REF][START_REF] Bernot | A fruitful application of formal methods to biological regulatory networks: Extending Thomas' asynchronous logical approach with temporal logic[END_REF] and to select exhaustively all suitable models. The chapter is organized as follows. Section 2 introduces the formalism due to Thomas for modelling the dynamics of a biological regulatory network. The resulting dynamics corresponds to a Kripke structure which can be deduced easily from the interaction graph. Section 3 describes the link between this transition system and the dynamics obtained with the classical modelling using piece-wise linear ordinary differential equation systems. Section 4 explains how formal methods can improve the modelling process of regulatory networks. The temporal properties have first to be translated into a temporal specification language. Then one has to answer automatically the question: does a given model satisfy the given temporal specifications ? Model checking makes this stage automatic and its principle is also presented in this section. Section 5 illustrates the use of model checking to model the well studied genetic regulatory network of temperate bacteriophage lambda rapidly described before. A model of this system has already been proposed by Thieffry and Thomas in [START_REF] Thieffry | Dynamical behaviour of biological regulatory networks -II. Immunity control in bacteriophage lambda[END_REF]. We show that our approach, using model checking, automatically selects this model as well as other models satisfying the same criteria of validation.

Qualitative dynamics of biological regulatory networks

The multi-valued modelling of Thomas is able to represent the qualitative dynamics of biological regulatory networks whose entities can be molecules, macromolecules, cells, organs, or organisms, if no societies. In fact, all systems whose regulations have a sigmoid shape can be modelled in this formalism. The regulations of genetic regulatory networks have almost always a sigmoidal nature, that explains why this formalism has been introduced in this context and why its main application domain remains the genetic regulatory networks. In such systems, the concentration of a protein encoded by a gene u may activate or inhibit the synthesis of proteins encoded by other genes or itself (figure 1). Figure 1: A genetic regulatory network. The gene u synthesizes a protein which activates the expression of gene v and itself by binding the promoters P v and P u respectively. In turn, the protein of gene v inhibits the expression of gene u when it binds P u . Then, the arrow from a gene to its protein represents the transcription and translation processes and the arrow from a protein to a promoter abstracts the diffusion and the fixation of the protein on the promoter.

If the protein of u activates (resp. inhibits) the expression of a gene v, we said that u is a positive (resp. negative) regulator of v. In such situation an increasing of the concentration of the protein encoded by u induces an increasing of the rate of synthesis of the protein encoded by v. Generally, the relation between the concentration of a regulator and the rate of synthesis of its target is, as we have seen before, sigmoidal. When the sigmoid is steep, as in figure 2-(a), u has a little effect on v if it is below the concentration threshold θ uv and at higher concentration a plateau is reached representing the maximal rate of synthesis of v under the effect of u. Naturally, for an negative regulator, the sigmoid is decreasing. This section presents successively how the regulations can be summarized into a regulatory graph corresponding to the static part of the modelling, then introduces the parameters which encode the effects of regulators on their targets in all possible situations, and finally presents how the dynamics can be deduced from these parameters.

Biological regulatory graphs

To formally define the static part of biological regulatory networks, we use labelled directed graphs. Vertices represent the biological entities of the network and arcs their regulations.

Definition 1 [Biological regulatory graph] A biological regulatory graph is a labelled directed graph G=(V,E) where

• each vertex v of V, called variable, is provided with a boundary b v ∈IN less or equal to the out-degree of v in G.

• each arc u → v of E is labelled with a couple (t uv ,α uv ) where t uv is an integer between 1 and b v , called qualitative threshold and where α uv ∈{+,-} is the sign of the regulation. Moreover it is required that for any variable u with b u > 0, ∀ i∈{1,2,...,b u }, there exists a successor v of u such that t uv = i.

In a biological regulatory graph G, the set of the regulators of a variable v corresponds to the set of its predecessors in G, denoted by G -(v), and the set of its targets corresponds to the set of its successors in G, denoted by G + (v). Each regulation u → v is labelled by a sign α uv which indicates if u is an activator or an inhibitor of v, and by a qualitative threshold t uv . Thresholds t are integers and do not correspond to biological thresholds θ∈IR, most often difficult to measure, but they give the order of the continuous thresholds: if t uv = i then the corresponding continuous threshold θ uv is the i th lowest threshold among {θ uv | v∈G + (u)}. That explains the requirement on qualitative thresholds of the previous definition which implies that b v is the number of different thresholds "outgoing" from v. gives an example of biological regulatory graph which can be deduced from the genetic regulatory networks described in figure 1. Figure 2 assumes that θ uv < θ uu , and consequently t uv = 1 and t uu = 2. Obviously concentration levels are associated to variables. For describing the evolution of the concentration level of each variable, it is necessary to know which regulators have an effect on the variable. Only the position of the regulator concentration with regard to their thresholds is sufficient. The concentrations are then discretized according to thresholds and each variable can take a finite number of values called abstract qualitative levels. For example, in figure 2-(b), the variable u has three different behaviours with regard to its targets:

• In the first region (the concentration of u is less than θ uv ), u acts neither on v nor on itself. • In the second region (the concentration of u is between θ uv and θ uu ), u acts on v but it still not act on itself. • In the last region (the concentration of u is greater than θ uu ), u acts both on v and on itself.

Three qualitative levels emerge, 0, 1 and 2, corresponding to the three previous regions and constitute the only relevant information from a qualitative point of view. More generally, a variable v can take b v +1 qualitative levels, from 0 to b v , and the qualitative level q means that v acts on all targets v' such that t vv' ≤ q. A state of the system is then defined as a vector constituted by qualitative levels of variables.

Definition 2 [Qualitative state] Let G=(V,E) be a biological regulatory graph. A qualitative state of G is a vector q=(q v ) v∈V such that for all v∈V, q v ∈{0,1,…,b v }. The set Q of states of G is then defined by Q=∏ v∈V {0,1,...,b v }.

In the sequel, we write v = l for denoting q v = l if it does not cause confusion.

Models of biological regulatory graphs

Remember that the sigmoid nature of a regulation u → v leads to distinguish two different situations: if u ≥ t uv , then the regulation is active and if u < t uv , it is not. We said that u is a resource of v when u induces an increasing of v:

• if the regulation u→v is positive, u is a resource of v when the regulation is active,

• if the regulation u→v is negative, u is a resource of v when the regulation is not active.

From the point of view of resources, the absence of an inhibitor acts as the presence of an activator.

Definition 3 [Resources]

Let G=(V,E) be a biological regulatory graph, v∈V and q∈Q. The set ω v (q) of resources of v at the state q is the subset of G -(v) defined by: ω v (q)={u∈G -(v) | (q u ≥ t uv and α uv = +) or (q u < t uv and α uv = -)}.

At the state q, the evolution of variable v depends on its resources ω v (q). It remains to define in which direction evolves v at the state q. The parameter K v,ω v (q) , called the attractor of v when the resources are ω v (q), denotes the level towards which v is attracted:

• if q v < K v,ω v (q)
then v tends to increase,

• if q v = K v,ω v (q) then v does not evolve and

• if q v > K v,ω v (q) then
v∈V, • ω ⊆ G -(v), • K v,ω ≤ b v .
It is often additionally required that:

K v,ω ≤ K v,ω'
for all v∈V, and for all ω,ω' ⊆ G -(v) such that ω ⊆ ω'.

These constraints, called Snoussi's constraints in the remainder, mean that the more a variable has resources the greater is the level towards which it is attracted. In other words, neither the presence of an activator nor the absence of an inhibitor can induce a decrease of the considered target (see the following section for the mathematical grounds of these constraints). This property, as well as signs of regulations, can often be deduced from biological knowledge and when it can be used, the number of models to consider for a given biological regulatory graph decreases drastically.

Model

K u,{} = 0 K u,{u} = 2 K u,{v} = 0 K u,{u,v} = 2 K v,{} = 0 K v,{u} = 1 u v Attractors Tendencies 0 0 K u,{v} = 2 K v,{} = 0 0 1 K u,{} = 0 K v,{} = 0 1 0 K u,{v} = 2 K v,{u} = 1 1 1 K u,{} = 0 K v,{u} = 1 2 0 K u,{u,v} = 2 K v,{u} = 1 2 1 K u,{u} = 2 K v,{u} = 1

Dynamics of models

The classical approach to describe the dynamics of models is to define the state of the system at time t +1 from its state at time t. One possibility is to consider that the next state is directly the attractor of the current state: if q is the current state then q' = (K v,ω v (q) ) v∈V is the next one and we said that there is a transition from q to q' (figure 5-(a)). This description raises serious problems for its application to biological systems:

1. From any initial state, the system will follow a well-defined path, without any branching or possibility of choice whereas biological systems typically include choices among several pathways (as illustrated for example by the numerous different pathways leading to various cell lines from a zygote during embryonic development).

2. Suppose that v is a gene which can take two values (b v = 1) and that the current state is q. If q v = 0, then K v,ω v (q) = 1 means that resources of v induce the production of the corresponding protein. This protein will appear after a time delay corresponding, for example, to the time of diffusion of its regulators (figure 4). Similarly the same phenomenon is observed when q v = 1 and K v,ω v (q) = 0 with an a priori different delay.

However, when q differs from q' = (K v,ω v (q) ) v∈V by at least two components, the corresponding variables change simultaneously (dashed arrow in figure 5-(a)). This synchronous description thus assumes that time delays are equal which is unlikely.

v (K v,ω v (q u ) = K v,{u} = 1
). The resulting protein appears after the delay v t ∆ . Finally, the protein u disappears, the gene v is no more stimulated (K v,ω v (q u ) = K v,{} = 0), and the protein v disappears after the different delay v t ∆ .

If the attractor of a variable is sufficiently away from its current value, one can have |q

v -K v,ω v (q) | > 1.
In such cases the qualitative level increases abruptly and jumps several thresholds (dotted arrow in figure 5-(a)). Since the dynamics of the model abstracts a continuous phenomenon, during a transition, each variable can pass through at most one threshold.

These points lead us to introduce the following asynchronous description.

Definition 5 [Asynchronous state graph] Let G=(V,E) be a biological regulatory graph and M(G) be a model of G. The asynchronous state graph of M(G) is a directed graph whose set of vertices is the set Q of states of G, and such that there is an edge from q to q' if:

• for all variables v∈V, q v = q' v = K v,ω v (q) or • there exist a variable v∈V such that:

o for any variable u ≠ v, q u =q' u and o q v < K v,ω v (q) and q' v = q v +1 or q v > K v,ω v (q) and q' v = q v -1.
In this definition, a state q which has itself as successor, is a stable steady state of the asynchronous state graph:

q v = K v,ω v (q)
for all v∈V. Otherwise, if q is a state for which n variables tend to evolve (n variables v such that q v ≠ K v,ω v (q) ), q has n successors and each of them differs from q by only one component corresponding to one of these n variables. Thus, when time delays are unknown, the asynchronous state graph contains all the a priori possible transitions. Some of them can be removed when time delays are taken into consideration.

Figure 5 shows the synchronous and asynchronous dynamics of the model of table 1. The attractors are the same in both descriptions but paths differ: the asynchronous state graph contains a circuit (0,0) → (1,0) → (1,1) → (0,1) → (0,0) which is absent in the synchronous description.

Differential modelling

We have seen that the asynchronous dynamics is more suited than the synchronous one for describing biological regulatory networks. This section proves how this asynchronous description can be deduced from a discretization of a particular class of ordinary differential equation systems classically used for describing biological regulatory networks.

Ordinary differential equation systems

Classically, the dynamic of a biological regulatory graph G=(V,E) is modelled by ordinary differential equation systems [START_REF] Glass | The logical analysis of continuous non linear biochemical control networks[END_REF][START_REF] Snoussi | Qualitative dynamics of a piecewise-linear differential equations : a discrete mapping approach[END_REF][START_REF] Thomas | Biological Feedback[END_REF], called here underlying differential systems (UDSs for short) of G, whose form is

v x = f v (x) -? v x v ,
for all v∈V, where x = (x v ) v∈V is a vector, whose components x v ∈IR + give the concentrations of variables v. The vector x is called the quantitative state of G. The previous equations define the rate of change of each concentration x v as the difference of the synthesis rate f v (x) and the degradation rate λ v x v of v. The function f v expresses how the synthesis rate of v depends on the concentrations x u of its regulators u∈G -(v). It can be defined as,

f v (x) = k v + α ∑ u∈G -(v) k uv s α uv (x u ,θ uv ),
where:

• k v ∈IR + and k uv ∈IR +* are kinetic parameters,

• the function s α uv gives the effect of a regulator u and its target v. This function is usually a sigmoid depending on the sign α uv and on the quantitative threshold θ uv ∈IR +* of the interaction.

Since the qualitative thresholds t of G give the order of the continuous thresholds θ (see the section 2.1), it is required that for all targets v'∈G

+ (u) of u different than v, θ uv < θ uv' if t uv < t uv' and θ uv > θ uv' if t uv > t uv'
. By denoting θ u l the threshold(s) θ uv such that v∈ G + (u) and t uv =l we then have:

θ u 1 < θ u 2 < . . . < θ u 1 < θ u b u
The sigmoidal function s α uv is often approximated by a step function in order to make possible the analytical analysis of the system. s α uv is then defined as a Boolean function which indicates if u is or not a resource of v:

) , ( 1 ) , ( and if , 0 if , 1 ) , ( uv u uv u uv u uv u uv u x s x s x x x s θ - = θ θ < θ >    = θ + - +
Notice that these functions are not defined for x u = θ uv and that the system becomes a piecewise linear equation system. Figure 6 gives an example of an UDS of the biological regulatory graph of figure 3.

     × × = × × + × + = + + - 2 ) 16 ( 25
5 ) 20 ( 40 ) 10 ( 35 20 v u v u u v u x - , x s x x - , x s , x s x

Discretization map and domains

Since the step functions s α uv are not defined for x u = θ uv , the differential equation system is not defined for the states x for which at least one component x v equals a threshold θ vv' , v'∈G + (v). Such states are called singular states. Consequently, the properties of the system can be analysed in the |V|-dimensional phase space Ω defined by

Ω = ∏ v∈V Ω v with Ω v = IR + \{θ vv' | v'∈G + (v)} for all v∈V.
Ω corresponds to the set of regular states. We are now in position to define the discretization

map d : Ω → Q by d(x) = (d v (x v )) v∈V with, for every v∈V, d v : Ω v → Q v defined by d v (x v ) = |{θ vv' | v'∈G + (v) and θ vv' < x v }|.
This discretization map gives directly the cardinal of the set of thresholds less than the concentration of v. If d v (x v ) = l, then x v is greater than the l smallest thresholds and less than others. For all v∈G + (u) we have

x u > θ uv ⇒ d u (x u ) ≥ t uv and x u < θ uv ⇒ d u (x u ) < t uv .
Consequently, for all state x∈ Ω:

s α uv (x u ,θ uv ) = 1 ⇒ u∈ω v (d(x)).
Then, for all x∈Ω, f v can be rewritten as

f v (x) = k v + ∑ u∈ω v (d(x)) k uv .
(

The infinite set of continuous states whose discretization gives q∈Q is an hyper-rectangular region D(q) of Ω, called domain, defined by:

D(q) = ∏ v∈V D v (q v ) with D v (q v )= {x v ∈ Ω v | d v (x v ) = q v } for all v∈V.
A domain D(q) is bounded by hyperplanes corresponding to thresholds: for all variable v, if q v >0 then θ v q v is the lower bound of D v (q v ) and if q v < b v then θ v q v +1 is the upper bound of D v (q v ) (see the figure 7).

Dynamics of differential equation systems

In a domain D(q), each function f v reduces to the constant k v +∑ u∈ω v (q) k uv (see equation 2).

The system thus simplifies to a linear and uncoupled differential equation system whose solution in D(q), starting at x 0 ∈D(q), is given by

x v (t) = κ v (x 0 ) -(κ v (x 0 ) -x v 0 ) e -λ v t
, for all v∈V,

with κ v (x) = f v (x)/λ v for all x∈Ω. Function κ v is also reduced to a constant (k v +∑ u∈ω v (q) k uv )/λ v in D(q), denoted κ v (q)
by abuse of notation. The state κ(q) = (κ v (q)) v∈V acts as an attractor in D(q). Indeed, it is easy to verify that in D(q), x(t) has the following properties:

1. if κ v (q)∈D v (q v ) then, x v (t) monotonically converges from x v 0 towards κ v (q) and reaches κ v (q) in infinite time. Thus, if κ(q)∈D(q), x(t) does not leave D(q) and the state κ(q) is the unique stable steady state in D(q).

2. if κ v (q)∉D v (q v ), then, if x v 0 < κ v (q) (resp. x v 0 > κ v (q)
), x v (t) monotonically increases (resp. decreases) from x v 0 until to reach the threshold value θ v q v +1 (resp. θ v q v ). The threshold is reached in a finite time iff κ v (q) is different from it. Throughout this section, we suppose that the parameters k and λ are taken such that κ(q)∈Ω for all q∈Q. Consequently, if κ(q)∉D(q), x(t) leaves in a finite time D(q) by reaching a threshold hyperplane.

If κ v (q)∉D v (q v ), x v (t) reaches its corresponding threshold at time t v given by

        - κ θ - κ λ - = α + 0 ) ( ) ( ln 1 v v q v v v v x q q t v with α = 1 if x v 0 < κ v (q) and α = 0 if x v 0 > κ v (q).
If at least two components x u (t) and x v (t) reach their thresholds simultaneously, one can deduce that x 0 belongs to an at most (|V|-1)dimensional surface of zero Lebesgue measure in D(q). Therefore, we do not consider this case and reason now as for almost every x 0 ∈D(q). Suppose that t v is the smallest value in {t u | κ u (q)∉D u (q u )}, in other words, suppose that v is the variable whose concentration first leaves the domain. The component x v (t) reaches

α + θ v q v
at the singular state x 1 given by

x v 1 = α + θ v q v and x u 1 = x u (t v ) for all u ≠ v.
At this time, the trajectory exits from the domain D(q) and enters into D(q') defined by:

• D u (q' u ) = D u (q u ) for all u ≠ v, since only v reached its threshold, • D v (q' v ) = D v (q v +1) if α = 1 and D v (q' v ) = D v (q v -1) if α = 0.
But at the singular state x 1 , the differential equation system is not defined as well as κ v (x 1 ). The linear differential equation system of D(q') is then extended by continuity to the hyperplane

x v = α + θ v q v
. Thus κ v (x 1 ) is defined, and the trajectory is extended with the solution of the differential system of D(q') from the new starting point x 1 .

Discrete and differential modellings are coherent

Let M(G) be a model of a biological regulatory graph G=(V,E). The UDSs of G such that :

d v (κ v (q)) = K v,ω v (q)
for all v∈V and q∈Q are called underlying differential systems of M(G). A model M(G) has UDSs if and only if it satisfies the Snoussi's constraints (equation 1) since we have:

) ( and all for , / , v G V v K k k d v v u uv v v - ω ∈ ⊆ ω ∈ =         λ         + ω ∑ and thus ω⊆ω' implies ∑ u∈ω k uv ≤ ∑ u∈ω' k uv which implies K v,ω ≤ K v,ω' .
For example the UDS of figure 6 is an UDS of the model described in Table 1. The following propositions show the coherence between the asynchronous dynamics of M(G) and the dynamics of its UDSs.

Proposition 1

• If there is an UDS of M(G) such that x∈D(q) is a stable steady state, then q is a stable state of the asynchronous state graph S of M(G). • Conversely, if q is a stable state of S then, for all UDSs of M(G), there is a stable steady state in the domain D(q).

Proof. A state x∈D(q) is a stable steady state iff x v = κ v (q) for all v∈V. That implies d v (x v ) = d v (κ v (q)) ⇒ q v = K v,ω v (q)
for all v∈V and thus, q is a stable state of S. Conversely, if q∈Q is a stable state, then

q v = K v,ω v (q) = d v (κ v (q)
) for all v∈V. Thus, κ v (q)∈D v (q v ) for all v∈V and consequently, κ(q)∈D(q) is a stable steady state.

We define now the boundary of a domain as the set of singular states whose distance to the domain is null.

Proposition 2

• If there is an UDS of M(G) for which there is a trajectory starting in D(q) which reaches directly from D(q) the hyperplane separating D(q) and an adjacent domain D(q'), then q → q' is a transition of the asynchronous state graph S of M(G). • Conversely, there exist UDSs of M(G) such that, for each successor q' of q in S, there is a trajectory starting in D(q) which reaches directly from D(q) the hyperplane separating D(q) and D(q').

Proof. We have seen in section 3.3 that if a trajectory starting at x 0 ∈D(q) reaches the hyperplane separating D(q) and an adjacent domain D(q'), then there is a unique variable v∈V such that q' v ≠ q v and we have q

' v = q v +1 if x v 0 < κ v (q) or q' v =q v -1 if x v 0 > κ v (q). Moreover, κ v (q)∉D v (q v ) thus x v 0 < κ v (q) iff d v (x v 0 ) < d v (κ v (q)) which is equivalent to q v < K v,ω v (q) .
Similarly

x v 0 > κ v (q) iff d v (x v 0 ) > d v (κ v (q)) which is equivalent to q v > K v,ω v (q)
. According to definition 5, q → q' is a transition of S. Now, we prove the second part of the proposition. Consider the UDSs of M(G) such that λ u =λ for all u∈V and an initial state x 0 ∈D(q). The trajectory starting at x 0 describes the part of the segment connecting x 0 to κ(q) which belongs to D(q). Let q' be a successor of q in S. We have κ(q)∉D(q). Let us choose a point x 1 of the boundary of D(q) belonging to the hyperplane separating D(q) from the domain D(q') and whose only one component equals a threshold. The trajectories starting at a point of the line connecting x 1 and κ(q) which belongs to D(q), reach x 1 .

We deduce from the previous propositions that all the regular stable steady states of an UDS of M(G) are represented in its asynchronous state graph S. Moreover if a trajectory of an UDS of M(G) passes successively through the domains D(q 0 ), D(q 1 ), ..., D(q n ) then q 0 → q 1 → ... → q n is a path of S. But if q 0 → q 1 → ... → q n is a path of S, it does not mean that there is a trajectory passing successively through the domains D(q 0 ), D(q 1 ), ..., D(q n ). Using the terminology of [START_REF] Kuipers | Qualitative reasoning: modeling and simulation with incomplete knowledge[END_REF], the qualitative modelling is said sound. A graphical comparison between the asynchronous dynamics of a model and a trajectory of one of its UDS is given in figure 8. Any UDS of a biological regulatory graph G is an UDS of a model of G satisfying the Snoussi's constraints. Thus the trajectories of the infinite set of UDSs of G are summarized by a finite set of asynchronous state graphs (for the biological regulatory graph of figure 3, we have 42 different state graphs deduced from the 60 different models satisfying the Snoussi's constraints).

Feedback circuit functionality

Feedback circuits play a major role for the dynamics of systems since they can generate multistationarity or homeostasis. A positive (resp. negative) circuit is said functional if it generates multi-stationarity (resp. homeostasis). The functionality of circuits is strongly related to the the stationarity of particular singular states and to discontinuities of the UDS. To deal with them, we first introduce the differential inclusion systems.

Differential inclusion systems

To deal with ordinary differential equation systems with discontinuous right-hand sides, Filippov [START_REF] Filippov | Differential equations with discontinuous right-hand sides[END_REF] proposed to extend them to systems of differential inclusions. For the regulatory networks, the UDSs can be extended to the following differential inclusions systems:

x v ∈H v (x), for all v∈V, (3)
where H v is a set-valued function, defined as follow:

• for all regular state x,

H v (x) = { f v (x) -λ v x v }. For all x∈D(q), since f v (x)/λ v = κ v (q), H v (x) can be rewritten as H v (x) = {λ v (κ v (q) -x v )}.
• for all singular state x,

H v (x) = ___ co ({λ v (κ v (q) -x v ) | q∈N(x)}).
• where ___ co (E) designs the smallest closed convex set of a set E which is the intersection of all closed convex sets containing E, and where N(x) is the set of qualitative states which correspond to domains whose boundary contains x:

. ) ( with if or 1 if ), ( , ) (                ∈ θ = - Ω ∈ = ∈ ∀ ∈ = + u G v x t t x x d q V u Q q x N uv u uv uv u u u u u • Obviously we have ( ) ( )       - κ λ - κ λ = ∈ ∈ v v v x N q v v v x N q v x q x q x H ) ( max , ) ( min ) ( ) ( ) (
Consider the example of figure 7. For x such that x u = θ u 2 = θ uu and x v > θ u 2 , we have N(x) = {(1,1), (2,1)} and for these states,

ω u (1,1) = {} is included in ω u (2,1) = {u} and ω v (1,1) = ω v (2,1) = {u}. We deduce that H u (x) = [λ u (κ u (1,1) -x u ), λ u (κ u (2,1) -x u )] and H v (x) = {λ v (κ v (1,1) -x v )}.
Intuitively, at the singular state x, the regulation u → v is clearly defined: s + (x u ,θ uv ) = 1. This is why the set H v (x) of the possible derivatives of x v is single-valued. However, as x u = θ uu the self regulation of u remains undefined and H u (x) is a priori not single-valued: the derivative of x u is comprised between the derivatives obtained with s + (θ uu ,θ uu ) = 0 and s + (θ uu ,θ uu ) = 1. An absolutely continuous function x(t) is solution of the system (3) in the sense of Filippov if ) (t x v ∈H v (x(t)) for all v∈V and for almost all t ≥ 0. The qualification ``for almost all t ≥ 0'' means that the set time-points for which the condition does not holds if of measure 0. In particular, the condition is not satisfied at time-points when the solution arrives or leaves a threshold hyperplane.

We do not analyse the solutions in the sense of Filippov in this section (see [START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF][START_REF] Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF] for a detailed analysis), but the previous formalism will be useful for analysis of the steadiness of singular states.

Steadiness of singular states

It is not surprising that a state x, regular or singular, is an equilibrium point (in the sense that there is a solution x(t) such that x(t) = x for all t ≥ 0) when 0∈H v (x) for all v∈V. For a regular state x∈D(q), we have, as for differential equation systems:

0∈H v (x) ⇒ 0∈{ λ v (κ v (q) -x v ) } ⇒ x v = κ v (q).
In this case, x is a regular stable steady state. For a singular state, the inclusion can be written as an inequality:

( ) (
)

) ( max ) ( min ) ( max 0 ) ( min 0 ) ( ) ( ) ( ) ( q x q x q x q x H v x N q v v x N q v v v x N q v v v x N q v κ ≤ ≤ κ ⇔ - κ λ ≤ ≤ - κ λ ⇔ ) ( ∈ ∈ ∈ ∈ ∈ and if x v ∉Ω v the inequality becomes strict: ) ( max ) ( min 0 ) ( ) ( q x q x H v x N q v v x N q v κ < < κ ⇔ ) ( ∈ ∈ ∈
because κ v (q)∈Ω v for all q∈N(x). Among all singular equilibrium points, those for which we have min q∈N(x) κ v (q) = x v = max q∈N(x) κ v (q) for x v ∈ Ω v , are singular steady states [24, 7]. Figure 9 shows a graphical representation of the conditions for the steadiness of singular states.

Proposition 3 Let x be a singular state and v a variable. If for all u∈G -(v) x u ≠ θ uv , then κ v (q) is constant for all q∈N(x).

Proof. For all u∈G -(v) we have x u ∈Ω u or x u = θ uv' ≠ θ uv with v'∈G + (u). In the first case, it is evident that q u = q' u for all q and q' in N(x). In the second case, for all q and q' in N(x), q u and q' u belong to {t uv' -1,t uv' } and t uv' ≠ t uv . Then q u and q' u are on the same side of t uv . Consequently, for all q and q' in N(x) we have ω v (q)=ω v (q') which implies κ v (q)=κ v (q').

Figure 9: Equilibrium points and their steadiness. Grey regions, a rectangle in (a) and a segment in (b), correspond to the Cartesian product Ψ(x)= [min q∈N(x) κ u (q), max q∈N(x) κ u (q)] × [min q∈N(x) κ v (q), max q∈N(x) κ v (q)] for a singular state x.

In (a) x = (θ uv ,θ vu ) is an equilibrium point (x∈Ψ(x)) and since all variables are singular, it is steady. In (b) the singular state is such that x u > θ uu and x v = θ vu and it is not an equilibrium point (x∉Ψ(x)).

Circuit characteristic states Definition 6 [Circuit] Let G=(V,E) be a biological regulatory graph. A circuit of G is a finite sequence of distinct elements of V, denoted C

= v 1 ,v 2 ,...,v n , such that v n → v 1 ∈E and v i → v i+1 ∈E for all i∈{1,...,n-1}.
In the sequel, {C} denotes the set of variables of a circuit C and i+1 (resp. i-1) is always computed modulo n: v i+1 (resp. v i-1 ) denotes the successor (resp. predecessor) of v i in C. Two circuits C and C' are disjointed if they have no variable in common. In a pedagogical objective, we focus here on the properties of a single circuit, but all results can be extended to a union of disjointed circuits [START_REF] Snoussi | Logical identification of all steady states : the concept of feedback loop characteristic states[END_REF]. Moreover we take into consideration only regulatory graphs where for any variable the out-thresholds are distinct (b

v = |G + (v)|, ∀ v∈V).
A singular state x is said characteristic of a circuit C = v 1 ,v 2 ,...,v n if the concentration x v i of each variable v i of the circuit is equal to the threshold θ v i v i+1 and if the concentrations of other variables are regular: x u ∈Ω u , ∀ u∉{C}. Proposition 4 A singular steady state is a characteristic state of a circuit of the biological regulatory graph G.

Proof. Let x be a singular state and S = {v | x v ∉Ω v } be the set of variables equal to a threshold at the state x. If x is steady, we have for all v∈S :

) ( max ) ( min ) ( ) ( q x q v x N q v v x N q κ < < κ ∈ ∈ .
According to the proposition 3, if for all u∈G -(v) we have x u ≠ θ uv then min q∈N(x) κ v (q) = max q∈N(x) κ v (q) and x is not steady. Thus v has at least one predecessor u such that x u = θ uv , which implies that u∈S. Moreover, because θ uv' ≠ θ uv for all v'∈G + (u), the successor v of u is the only one such that x u = θ uv . Each variable v of S has then a unique predecessor u in S such that x u = θ uv . Now, we prove that all the steady singular states can be identified in the qualitative modelling.

Proposition 5 Let G be a biological regulatory graph containing a circuit C

= v 1 ,...,v n .

Consider a UDS of a model M(G) and a characteristic state x of C. Let q∈N(x). If x is steady, then M(G) is such that

K v,ω v (q) = q v for all v∉{C}

K v i ,ω v i (q) \ {v i-1 } < t v i v i+1 ≤ K v i ,ω v i (q) ∪ {v i-1 } for all i∈{1,…,n}
Proof. Let v∉{C}. Since x is characteristic of C, we have x v ∈Ω v . If x is steady, then min q∈N(x) κ v (q) = x v = max q∈N(x) κ v (q). That means that κ v (q) is constant for all q∈N(x) and we have

d v (x v ) = d v (κ v (q)) which is equivalent to q v = K v,ω v (q) . Let v i ∈{C}. As x is characteristic of C, v i-1 is the unique predecessor of v i such that x v i-1 = θ v i-1 v i .
Thus, ω v i (q) \ ω v i (q') equals {} or {v i-1 } for all q and q' in N(x). Moreover there is at least one state q∈N(x) such that q v i-1 = t v i-1 v i and another one such that q v i-1 = t v i-1 v i -1. Thus, there is a state q + ∈N(x) such that v i-1 ∈ω v i (q + ) and a state q -with v i-1 ∉ω v i (q -). We deduce that for all q∈N(x), ω v (q)∪{v i-1 } = ω v i (q + ) and

ω v (q)\{v i-1 } = ω v i (q -). So max q∈N(x) κ v i (q) = κ v i (q + ) and min q∈N(x) κ v i (q) = κ v i (q -). Since x v i = θ v i v i+1 , if
x is steady, we have for all q∈N(x):

κ v i (q -) < θ v i v i+1 < κ v i (q + ) ⇔ d v i (κ v i (q -)) < θ v i v i+1 < d v i (κ v i (q + )) ⇔ K v i ,ω v i (q -) < t v i v i+1 ≤ K v i ,ω v i (q + ) ⇔ K v i ,ω v i (q) \ {v i-1 } < t v i v i+1 ≤ K v i ,ω v i (q) ∪ {v i-1 } Definition 7 [Quasi-characteristic qualitative states] Let G=(V,E) be a biological regulatory graph containing a circuit C = v 1 ,...,v n . A state q∈Q is quasi-characteristic of C if q v i = t v i v i+1 for all v i ∈{C}.
The quasi-characteristic states are useful to locate the singular characteristic states of the UDS.

Proposition 6 Let G be a biological regulatory graph containing a circuit C = v 1 ,...,v n and a quasi-qualitative characteristic state q of C. If a model M(G) satisfies the Snoussi's constraints and if K v,ω v (q) = q v for all v∉{C}

K v i ,ω v i (q) \ {v i-1 } < t v i v i+1 ≤ K v i ,ω v i (q) ∪ {v i-1 } for all i∈{1,…,n}
then, for all the UDSs of M(G), there exists a unique steady characteristic state x of C such that d u (x u ) = q u for all u∉{C}.

As the proof is quite similar to the previous one, it is omitted.

The previous proposition makes easy the determination of all steady singular states underlying of a qualitative model. Let us consider for instance the model M(G) of Table 1. The corresponding biological regulatory graph G (Figure 3-(a)) contains two circuits, C 1 = u,v and C 2 = u. The unique quasi-characteristic state of C 1 is (t uv ,t vu ). It satisfies

K u,ω u (t uv ,t vu ) \ {v} < t uv = 1 ≤ K u,ω u (t uv ,t vu ) ∪ {v} and K v,{} < t vu = 1 ≤ K v,{u} .
Indeed the first inequality is verified because ω u (t uv ,t vu ) = {}, K u,{} = 0 and K u,{v} = 2, the second is also verified since K v,{} = 0 and K v,{u} = 1. Consequently, the characteristic state uv ,θ uv ) is steady in all the UDSs of M(G). For circuit C 2 , there are two quasi-characteristic states: (t uu ,0) and (t uu ,1).

• The first one, (t uu ,0), does not satisfy K u,ω u (t uu ,0) \ {u} < t uu = 2 ≤ K u,ω u (t uv ,0) ∪ {u} and K v,ω v (t uu ,0) = 0. since ω v (2,1) = {u} and K v,{u} = 1. Thus there is not any steady characteristic state of C 2 such that x v < θ vu .

• The second quasi-characteristic state, (t uu ,1), satisfies

K u,ω u (t uu ,1) \ {u} < t uu = 2 ≤ K u,ω u (t uv ,1) ∪ {u} and K v,ω v (t uu ,1) = 1. since ω v (2,1) = {u}, K v,{u} = 1, ω u (2,1) = {u}, K u,{} = 0 and K u,{u} = 2. For all UDSs of M(G) there is a unique steady characteristic state x of C 2 such that x v > θ vu .
The detected singular states are represented in the asynchronous state graph of M(G) in figure 10.

(4) Figure 10: Representation of the steady singular states of model of Table 1.

Circuit functionality

Each variable of a feedback circuit has an influence on its target but also an indirect effect on all following variables including itself. A circuit is said positive (resp. negative) if each variable has a positive (resp. negative) influence on itself. The sign of a circuit is determined by the number of inhibitions: if it is odd, the circuit is negative and otherwise, the circuit is positive. Negative and positive circuits have different typical behaviours.

• In a negative circuit, a high level of a variable tends to make decrease itself and conversely. Thus the circuit makes the level of each variable to tend to (or oscillate around) an equilibrium concentration. It generates stable oscillation behaviour corresponding to homeostasis in biology. • In a positive circuit, a high (resp. low) level of a variable tends to make it increase (resp. decrease). Thus each variable stays either at a low or high concentration and the positive circuit generates multi-stationarity corresponding to differentiation in biology.

A circuit which presents a typical behaviour is said functional. Several authors have proved that at least one positive circuit is necessary to generate multi-stationarity whereas at least one negative circuit is necessary to obtain a stable oscillatory behaviour [START_REF] Plathe | Feedback loops, stability and multistationarity in dynamical systems[END_REF][START_REF] Snoussi | Necessary conditions for multistationarity and stable periodicity[END_REF][START_REF] Cinquin | Roles of positive and negative feedback in biological systems[END_REF][START_REF] Cinquin | Positive and negative feedback: striking a balance between necessary antagonists[END_REF][START_REF] Soulé | Graphical requirements for multistationarity[END_REF]. Snoussi and Thomas realized that when a characteristic state is steady, the corresponding circuit is functional [START_REF] Snoussi | Logical identification of all steady states : the concept of feedback loop characteristic states[END_REF]. In the qualitative formalism, the circuit functionality is then defined as follow.

Definition 8 [Functional circuit] Let M(G) be a model of a biological regulatory graph G containing a circuit C. If there is a quasi-characteristic state q of C satisfying the constraint (4) then C is functional.

We deduce from the proposition 6 that if a circuit C is functional, there is, for all underlying differential systems, a steady characteristic state x of C such that x u = d u (q u ) for all u∉{C}. In the model of Table 1, both circuits u → u and u → v → u are functional. As a result, multistationarity and homeostasis are present in the corresponding asynchronous state graph (figure 10). Summing up, homeostasis and/or multi-stationarity are dynamical properties almost always present in biological systems. Circuit functionality is then useful for modelling such systems. For example, it has been used to model immunity control in lambda phage [START_REF] Thieffry | Dynamical behaviour of biological regulatory networks -II. Immunity control in bacteriophage lambda[END_REF], pattern formation during the embryonic development of Drosophilae [START_REF] Sánchez | Establishment of the dorso-ventral pattern during embryonic development of drosophila melanogaster: a logical analysis[END_REF][START_REF] Sánchez | A logical analysis of the drosophila gap-gene system[END_REF] and flower morphogenesis in Arabidopsis thaliana [START_REF] Mendoza | Genetic control of flower morphogenesis in arabidopsis thaliana: a logical analysis[END_REF].

Formal methods

To study the behaviour of the genetic regulatory network, the ordinary differential equation systems are well adapted if all the parameters are well known. Unfortunately they are most often unknown and are difficultly measurable in vivo. The discrete approach of Thomas and co-workers simplifies the problem of determining the suitable parameters since the number of possible models is finite. Indeed finding suitable classes of those parameters constitutes a major issue of the modelling activity. Even if the Snoussi's constraints on parameters are used, the number of remaining models is too large to analyse them by hand. Then biological knowledge or hypotheses on the behaviour of the system can be used as an indirect criterion to constrain the parameters. For example homeostasis (resp. multi-stationarity) is experimentally observable and it indicates that a negative (resp. positive) feedback circuit is functional, this functionality leading to some constraints on the parameters (see section 3).

To go further, conditions of multi-stationarity and homeostasis can be reinforced by introducing other conditions on the dynamics of the system. The available knowledge on the evolution of the system, as temporal properties, can be taken into consideration for constraining the values of parameters. Among all suitable models only a part of them are coherent with these temporal properties. Since numerous models have to be checked against those properties, a formal language is needed to perform automatically these checkings.

Temporal logic

The properties as the deadlock can be easily checked by exploring the transition system, called asynchronous state graph in section 2. For more complex properties on the dynamics of the system it is necessary to use a well adapted formal language: a temporal language which allows the specification of properties along the execution paths of the transition system. The step of the specification of the properties can then be distinguished from the specification of the system since it is not necessary to know the dynamic structure of the system to be checked for specifying the properties.

Expressing temporal properties on a transition system needs to define the atomic propositions which depends of the considered regulatory graph G=(V,E). Generally the set of atomic propositions is denoted by AP. The subset of AP containing all the atomic propositions which are true in a state q, is given by the labelling function L:

L(q) = { (v = q v ) | v ∈ V }
where (v = q v ) signifies that the variable v has the concentration level q v . The pair composed of a transition system and a labelling function is called a Kripke structure.

Execution traces of the transition system model implicitly a discrete time: if an execution passes from the state s 0 to s 1 , the instant associated to the state s 1 follows the one corresponding to the state s 0 . The temporal logics allow one to specify dynamical properties referring to this discrete time [START_REF] Emerson | Handbook of theoretical computer science, Volume B : formal models and semantics, chapter Temporal and modal logic[END_REF]. The Linear Temporal Logic, LTL, is used to specify properties on an execution of the system. If the system is determinist, from any initial state there is a unique execution, LTL is appropriated to specify properties of the system. Nevertheless the qualitative behaviour of a biological regulatory network is represented by an asynchronous state graph which is non determinist: the current state can have several possible futures. Since time has a tree structure, we prefer the Computation Tree Logic, CTL, in which it is possible to express properties of the form ''it is possible in the future that...''.

Definition 9 [Syntax of CTL] A CTL formula on the set of atomic propositions AP is inductively defined by:

• and any atomic proposition of AP are formulae • if φ and ψ are formulae, then (¬φ), (φ∧ψ), (φ∨ψ), (φ⇒ψ), (φ⇔ψ), AXφ, EXφ, A[φUψ], E[φUψ], AGφ, EGφ, AFφ, EFφ are formulae.

The semantics of CTL is defined on the execution trees of the transition system which are completely defined by their initial state and the transition relation. The semantics is given by the definition of the satisfaction relation s φ meaning that the formula φ is satisfied on the execution tree starting at s.

Definition 10 [Semantics of CTL] Let s 0 be a state. The semantics of CTL is defined inductively by:

• s 0 and s 0 , • ∀ p∈AP, s 0 p iff p∈L(s 0 ),

• s 0 ¬ϕ iff s 0 ϕ, • s 0 ϕ 1 ∧ϕ 2 iff s 0 ϕ 1 and s 0 ϕ 2 , • s 0 ϕ 1 ∨ϕ 2 iff s 0 ϕ 1 or s 0 ϕ 2 , • s 0 ϕ 1 ⇒ϕ 2 iff s 0 ϕ 1 or s 0 ϕ 2 , • s 0 ϕ 1 ⇔ϕ 2 iff s 0 (ϕ 1 ⇒ϕ 2 )∧(ϕ 2 ⇒ϕ 1 ),
• s 0 AXϕ iff for all successors s 1 of s 0 , we have s 1 ϕ, • s 0 EXϕ iff for any successor s 1 of s 0 , we have s 1 ϕ, • s 0 AGϕ iff for all paths s 0 ,s 1 ...s i ..., and for all s i along the path we have s i ϕ, • s 0 EGϕ iff for a particular path s 0 ,s 1 ...s i ... we have for all s i along the path s i ϕ, • s 0 AFϕ iff for all paths s 0 ,s 1 ...s i ..., there exists s i along the path such that s i ϕ, • s 0 EFϕ iff for a particular path s 0 ,s 1 ...s i ..., there exists s i along the path such that s i ϕ, • s 0 A[ϕ 1 Uϕ 2 ] iff for all paths s 0 ,s 1 ...s i ..., there exists s i along the path such that s i ϕ 2 and for each j<i we have s j ϕ 1 ,

• s 0 E[ϕ 1 Uϕ 2 ] iff for a particular path s 0 ,s 1 ...s i ..., there exists s i along the path such that s i ϕ 2 and for each j<i we have s j ϕ 1 .

is the always true formula; is the always false formula; a state s satisfies all the atomic formulae of L(s); ¬, ∧, ∨, ⇒, ⇔ are the usual connectives (respectively not, and, or, implication, equivalence). All the temporal connectives are pairs of symbols: the first element is A or E followed by X, F, G or U whose meanings are given in the next table and illustrated in Figure 11.

A for All paths choices X neXt state E for at least one path choice (Exist) F some Future state G all future states (Globally) U Until

Consider the example of Figure 5-(b) where variables are u and v. The atomic proposition are AP={(u = 0),(u = 1),(u = 2),(v = 0),(v = 1)}. AX(v = 1) means that in all next states accessible from the current state in the asynchronous state graph, the concentration level of v is 1. This formula is true iff the current state is (1,1), (2,0) or (2,1). EG(¬(u = 2)) means that there exists at least one path starting from the current state where the concentration level of u is constantly strictly less than 2. In Figure 5-(b), all states for which u is strictly less than 2 satisfy the formula. Then ¬(u = 2) ⇒ EG(¬(u = 2)) is satisfied for all states. A[(v = 1)U(v = 0)] means that for any possible path from the current state there exists a future state where v = 0 and in between v remains equal to 1. Note that (2,1) is the only state which does not satisfy the formula. And so on for other temporal connectives.

It is now possible to translate a biological temporal property into a CTL formula. Classically a biological system can have several steady states corresponding to distinct phenotypes. Let us suppose that two distinct stable states, ss 1 and ss 2 , are possible and that formulae ψ 1 and ψ 2 characterize the states ss 1 and ss 2 respectively. If the system is able to go from state s 0 , characterized by the formula ϕ 0 , either to state ss 1 or to state ss 2 , these temporal properties can be translate into formulae:

ψ 1 ⇒ AG ψ 1 stability of state ss 1 ψ 2 ⇒ AG ψ 2 stability of state ss 2 (ϕ 0 ⇒ EFψ 1 ) (ϕ 0 ⇒ EFψ 2 )
reachability of ss 1 and ss 2 from s 0 Such formulae are used in the concrete example of section 5 for expressing biological knowledge on the immunity control in bacteriophage lambda.

Model checking

The model checking is a verification method that proves automatically if a Kripke structure satisfies a temporal formula [START_REF] Huth | Logic in Computer Science: Modelling and reasoning about systems[END_REF]. We briefly present the basic algorithm of model checking for a CTL formula. Since the connectives ∨, ⇒ and ⇔ can be rewritten in term of ¬ and ∧, and since we have the following equivalence: AXϕ ≡ ¬EX(¬ϕ)

EGϕ ≡ ¬AF(¬ϕ) EFϕ ≡ E( U ϕ) AGϕ ≡ ¬EF(¬ϕ) A[ϕ 1 Uϕ 2 ] ≡ ¬ ( E[¬ϕ 2 U(¬ϕ 1 ∧¬ϕ 2 )] ∨ EG(¬ϕ 2
) ) we consider in the sequel formulae containing only the connectives: ¬, ∧, EX, AF and EU. Obviously any CTL formula can be transformed into a semantically equivalent CTL formula which uses only those connectives.

The model checking for a CTL formula ϕ consists in labelling each state s of the transition system with sub-formulae of ϕ which are satisfied at the state s. These sub-formulae are added to L(s) containing initially the atomic propositions true in s. Suppose that ψ is a subformula of ϕ and that states satisfying all the immediate sub-formulae of ψ have already been labelled. The labelling algorithm for ψ uses a case analysis to label states with ψ:

• if ψ∈AP, then the labelling is given directly by L(s)

• if ψ = p∧q, then L(s) = L(s) ∪ {p∧q} for all s such that p, q∈L(s) • if ψ = ¬p, then L(s) = L(s) ∪ {¬p} for all s such that p∉L(s) • if ψ = EXq, then L(s) = L(s) ∪ {EXq} for all predecessors s of a state t such that q∈L(t) • if ψ = AFq, then 1. L(s) = L(s)
∪ {AFq} for all s such that q∈L(s) 2. Repeat: L(s) = L(s) ∪ {AFq} for all states s such that all successors are labelled with AF q, until there is no change. It can be proved that this labelling algorithm ends and that states are labelled with all subformulae of ϕ that they satisfy. Thus s ϕ if the state s is labelled with ϕ. By extension if all states are labelled with ϕ, we say that the considered Kripke structure satisfies ϕ.

• if ψ = E[
The model checking algorithm is linear with the size of the system and the size of the formula. Unfortunately, practical applications lead to transition systems with an enormous number of states, and the previous algorithm is often inefficient. To push back these limits, symbolic model checking [START_REF] Mcmillan | Symbolic Model Checking[END_REF] has been developed. It consists in computations on symbolic representation of subspaces of states.

To sketch the symbolic model checking, let us introduce the operator Pre. Let S be the set of states and x be a subset of S. Pre(x) gives the set of states which have a successor in x. The set sat(ϕ) of states satisfying ϕ can then be defined inductively:

• if ϕ∈AP, sat(ϕ) = {s∈S | ϕ∈L(s)} • sat(¬ϕ) = S \ sat(ϕ) • sat(ϕ ∨ ψ) = sat(ϕ) ∪ sat(ψ) • sat(ϕ∧ψ) = sat(ϕ) ∩ sat(ψ) • sat(EXϕ) = Pre(sat(ϕ))
• sat(AXϕ) = S \ Pre(S \ sat(ϕ))

• The connectives AFϕ and E[ϕ 1 Uϕ 2 ] are more difficult to define. Let us remark that we have the following equivalence:

AF ϕ ≡ ϕ ∨ (AX (AFϕ)) E[ϕ 1 Uϕ 2 ] ≡ ϕ 2 ∨ ( ϕ 1 ∧ EX(E[ϕ 1 Uϕ 2 ]) ).
Then sat(AFϕ) and sat(E[ϕ 1 Uϕ 2 ]) can be defined as the smallest fixed points of equations:

f 1 (x) = sat(ϕ) ∪ sat(AX x) f 2 (x) = sat(ϕ 2 ) ∪ (sat(ϕ 1 ) ∩ sat(EX x )).
Since functions f 1 and f 2 are monotone and that the set of states is finite, the iterative computation of the smallest fixed point ends.

The Binary Decision Diagrams, or BDD for short, are data structures allowing the representation of Boolean expressions in a very compact way. Then subsets of states can be coded with such Boolean expressions and necessary operations for computing sat can be defined on these structures. Numerous works detail utilization of BDDs for the verification of systems, see for example [START_REF] Mcmillan | Symbolic Model Checking[END_REF][START_REF] Huth | Logic in Computer Science: Modelling and reasoning about systems[END_REF].

A tool for the selection of models: SMBioNet

We have designed a software for a computer aided modelling based on the previous described formal methods [START_REF] Bernot | A fruitful application of formal methods to biological regulatory networks: Extending Thomas' asynchronous logical approach with temporal logic[END_REF]. This software, SMBioNet1 , helps the biologist and/or the modeller to verify systematically the coherence of models of a given biological system, and to select suitable models which satisfy the temporal properties extracted from knowledge or hypothesis. More precisely inputs of SMBioNet consist in • a biological regulatory graph representing the interactions of the biological system and • a CTL formula expressing its known or hypothetical dynamical properties.

Then it generates all the models of the biological regulatory graph and gives as output those satisfying the CTL formula. For each generated model, SMBioNet calls the model checker NuSMV [START_REF] Cimatti | NuSMV: a reimplementation of SMV[END_REF] and selects it if the formula is satisfied. For each selected model, the asynchronous state graph and the steady states (regular and singular) are given. Depending on the available biological knowledge, the user can • reduce the domain of variation of some parameters,

• apply general constraints on parameters as, for example, the Snoussi's and observability 2 constraints, • specify a set of steady states (regular and singular) and a set of functional circuits.

These direct constraints on parameters decrease significantly the number of models to generate and consequently increase the efficiency of the selection. However, one can test directly the coherence of the regulatory graph (i.e. is there at least one suitable model ?), without enumeration of models by using a symbolic description of the set of all models.

In the next section, we shows how SMBioNet can be used for modelling the immunity control in bacteriophage lambda.

Immunity control in bacteriophage lambda

One of the most studied genetic regulatory networks is probably the one controlling immunity in temperate bacteriophage lambda which is a temperate virus. As described in figure 12, after infection of a bacterial population, many bacteria soon lyse and produce new phages but some survive and carry lambda genome in a dormant form. The first response is called lytic and the second lysogenic. In the lysogenic bacteria, viral DNA has integrated into the bacterial chromosome and will be faithfully transmitted to the bacterial progeny. In this condition, the viral gene cI, produces a repressor which blocks the expression of all the other genes of the phage, thus making the viral genome harmless for the bacterium. Moreover, cI makes lysogenic bacteria immune towards other infections. Lysogenization necessitates two events, integration of the viral DNA into the bacterial chromosome and development of immunity due to the expression of the repressor. The choice between the lytic and lysogenic pathways is very similar to cell differentiation, in the sense that a given virus, infecting apparently identical cells, can behave in two extremely different ways. It is actually in the context of this biological system that Thomas started to develop his formalism. Although he proposed various models of the immunity control [START_REF] Thomas | A complex control circuit. Regulation of immunity in temperate bacteriophages[END_REF][START_REF] Thomas | Biological Feedback[END_REF][START_REF] Thieffry | Dynamical behaviour of biological regulatory networks -II. Immunity control in bacteriophage lambda[END_REF], we focus in this section on the model developed by Thieffry and Thomas in [START_REF] Thieffry | Dynamical behaviour of biological regulatory networks -II. Immunity control in bacteriophage lambda[END_REF], which is denoted M(G) in the sequel. We will show that SMBioNet allows one to select, automatically and with very few biological knowledge, a set of models containing M(G) and satisfying the validation criteria given by Thieffry. All models of this set have to be considered since they have a priori the same prediction capacity than M(G).

Biological regulatory graph

The biological regulatory graph G summarizes the main regulations of the immunity control (Figure 13). Obviously it contains gene cI, but also three others (cro, cII, and N) which play a predominant role. Gene cI is activated by cII. Once on, gene cI remains on because its product activates its own synthesis, but at the same time, gene cI switches off the other lambda genes, including cII which had just switched it on. In addition gene cro exerts a negative control on cI, directly and indirectly, by repressing gene cII. Finally, gene N exerts a positive control on cII and is itself under negative control of cI and cII. According to the thresholds fixed by Thieffry, variables cI, cro, cII and N are 3-,4-,2-and 2-valued respectively, leading to 48 possible states. In the remainder, the state of the system is represented by the vector (cI,cro,cII,N).

Temporal properties

When the viral genome integrates a cell, all the viral proteins are initially absent. Thus (0,0,0,0) corresponds to the initial state of the system. The existence of both responses, lytic and lysogenic, implies that there exist two paths starting from the initial state leading respectively to the lytic state and to the immune one. The lytic state is known to be characterized by high concentration of cro and a low concentration of cI, cII and N whereas immune state is characterized by high concentration of cI and low concentration of cro, cII and N. In [START_REF] Thieffry | Dynamical behaviour of biological regulatory networks -II. Immunity control in bacteriophage lambda[END_REF], both states (0,2,0,0) and (0,3,0,0) correspond to the lytic state and (2,0,0,0) is the only state corresponding to the immunity. Without change of the environment, the choice between the lytic and the lysogenic pathways is irreversible, thus the lytic and immune states are steady. Then if the system reaches one state of the sets A={(0,2,0,0),(0,3,0,0)} or B={(2,0,0,0)}, then it will never leave it. These sets of states are said steady sets.

Summing up, dynamics of models to consider have to contain paths from (0,0,0,0) to the steady sets of states A and B. These properties are translated into the CTL formula Φ as follow:

init = ((cI = 0) ∧ (cro = 0) ∧ (cII = 0) ∧ (N = 0)) lytic = ((cI = 0) ∧ (cro ≥ 2) ∧ (cII = 0) ∧ (N = 0)) immune = ((cI = 2) ∧ (cro = 0) ∧ (cII = 0) ∧ (N = 0)) Φ A = lytic ⇒ AG(lytic) Φ B = immune ⇒ AG(immune) Φ r = init ⇒ ( EF(lytic) ∧ EF(immune) ) Φ = Φ A ∧ Φ B ∧ Φ r
The sub-formulae init, lytic and immune characterize the initial state, and the sets A and B. The steadiness of A and B is translated by Φ A and Φ B . The formula Φ r expresses reachability of A and B from the initial state and Φ represents the temporal properties to use for the selection of models.

Selected models

There is near 7 thousands of millions of models associated to G leading to about 3 millions of different asynchronous state graphs. If we consider the Snoussi's constraints (equation 1) as Thieffry and Thomas did, it remains 151200 models. Moreover, we use the activity constraints [START_REF] Bernot | Semantics of biological regulatory networks[END_REF]:

for each regulation u → v there is a set ω⊂G -(v) such that K v,ω ≠ K v,ω∪{u}
which stands for the observability of any regulation. If u → v does not satisfy the constraints, the attractor of v does not depend on the level of u. It seems then quite obvious that any model should satisfy this property in order that all regulations play a role in the dynamics. Taking into account these constraints, SMBioNet selects among the 882 remaining models, 88 models satisfying the formula Φ. The model M(G) proposed by Thieffry and Thomas is one of them. Table 2 shows the possible values of parameters for the selected models. 17 parameters among 24 are fixed by formula Φ (in particular, all the parameters associated to N).

K cI,{} = 0 K cI,{cI} = 1 or 2 K cI,{cro} = 0,1 or 2 K cI,{cII} = 0,1 or 2 K cI,{cI,cro} = 2 K cI,{cI,cII} = 1 or 2 K cI,{cro} = 2 K cI,{cI,cro,cII} = 2 K cro,{} = 0 K cro,{cI} = 2 K cro,{cro} = 0 K cro,{cI,cro } = 2 or 3 K cII,{} = 0 K cII,{cI} = 0 K cII,{cro} = 0 K cII,{N} = 0 K cII,{cI,cro} = 0 or 1 K cII,{cI,N} = 0 or 1 K cII,{cro,N} = 0 or 1 K cII,{cI,cro,N} = 1 K N,{} = 0 K N,{cI} = 0 K N,{cro} = 0 K N,{cI,cro } = 1
Table 2: Possible values of parameters for the selected models. Bold numbers correspond to the model M(G).

Validation of models

Thieffry and Thomas exhibited one model whose coherence is analysed through the likelihood of some paths of the asynchronous state graph of M(G) and through the pertinence of predictions on the dynamics of some mutants. Our approach leads to select 88 models which have to be evaluated with the same biological criteria of validation.

• Although 4 positive feedback circuits are present in the regulatory graph, the 88 selected models present only two steady states (regular or singular): (2,0,0,0) is always steady and the other one is either (0,2,0,0) or a singular state adjacent to (0,2,0,0) and (0,3,0,0). These steady states correspond to the lytic and immune states, and no other stable behaviour (phenotype) can be observed.

• Even if several pathways are possible from the initial state to immune state, all selected models present the most likely pathway in M(G) from initial state to A (see Figure 14). Similarly the pattern of dynamics present in M(G) allowing the system to evolve from initial state to lytic state, is also present in all selected models.

• Biological knowledge on mutants is available and can be used for validating models.

The considered mutations correspond to the inactivation of different combinations of genes. Then simulations of the behaviour of these mutants can be performed and confronted to the biological knowledge. For example, the dynamics of the mutant λ cIcro -, where genes cI and cro are inactivated, is obtained from M(G) by setting to 0 all parameters associated to cro or cI. Consequently, from an initial state where cI and cro are absent, they will never appear. The dynamics of this mutant is given in Figure 15. The dynamics of mutants obtained from M(G) are coherent from a biological point of view, since the remaining basins of attraction allow the prediction of the behaviour of mutants. For any selected model, results are the same and are given in table 3.

Among the 88 selected models, some differences can be highlighted. For example, 2 states are unreachable from the initial state in M(G) whereas for some selected models 15 states are unreachable. In such models, all states with cI = 2 and cro = 3 are not reachable, which is reasonable because high concentration of cI and cro is rarely observed. Moreover, such models do not contain the path (0,0,0,0) → (1,0,0,0) → (2,0,0,0) present in the dynamics of M(G) and which is unlikely in view of the low expression of cI when cII is absent.

Mutants Basins of attraction λ

cI - A λ cro - B λ cII - A and B λ N - A and B λ cI -cro - {(0,0,1,1)} λ cII -N - A and B λ cro -N - A λ cro -cII - A λ cro -cII -N - A
Table 3: Basins of attraction for a collection of mutants.

In conclusion, these 88 selected models satisfy the same criteria of validation that M(G) and have also to be considered. These models have been selected using a formula Φ expressing the well known properties of the system. Thieffry and Thomas have exhibited their model with the circuit functionality and some hypothesis. We can notice that the used constraints for functionality are not necessary to reproduce the biological properties (expressed by Φ) because some of the models selected by Φ do not satisfy these functionality constraints. Moreover some parameters are valuated according to hypotheses (K cII,{cI,N} = 2 for example) which have to be slacken since some models selected by Φ propose different values for these parameters.

Conclusion

We have defined a formal description of biological regulatory networks which allows a computer aided manipulation of the semantics of the discrete modelling of Thomas, this manipulation being proved correct by construction. Our approach allows biology to take advantage of the whole corpus of formal methods from computer science. Model checking is a first powerful tool offered by the formalization of biological regulatory networks. In particular, temporal properties can be added into the specifications of the system, and the modelling task consists in exhibiting one or more generally all models that are coherent with the previous specifications expressing a part of the biological knowledge concerning the dynamics of the system. All potential models have to be checked against temporal formulae, and this task can be done automatically using model checking. This brute force approach permits one to exhibit exhaustively all suitable models, i.e. all models satisfying the temporal formulae. Information provided by a new experiment or a new theoretical point of view will refine the set of selected models. The available temporal properties concern generally the homeostasis, the multistationarity, stable steady states and the accessibility of some stable steady states from a partially specified initial state. Unfortunately the stable steady states are some time singular and not formally represented in the asynchronous state graph of Thomas. Then the specifications cannot easily contain temporal properties concerning such singular states. This would necessitate to rewrite these temporal properties with only atomic propositions of regular states, and this task is generally difficult.

De Jong et al. [START_REF] Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF] introduced the singular states into their qualitative dynamics. Their qualitative modelling of genetic regulatory networks is also based on piecewise-linear differential equations. Authors propose a mathematically well founded method to deal with singular states using differential inclusions [START_REF] Filippov | Differential equations with discontinuous right-hand sides[END_REF][START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF]. Our approach consisting in adding temporal properties into the specifications for determining the suitable parameter values, would allows in this context to treat regular states as well as singular states.

More generally the formal methods can be applied in the field of biological regulatory networks and systems biology in order to explicit some behaviours or to take into account biological knowledge which have been ignored for the moment. The cooperation between biology and formal methods from computer science opens a large horizon of research perspectives.

• The introduction of transitions in the regulatory graph could help to specify how the different regulators cooperate for inducing or repressing their common target [START_REF] Bassano | Marked regulatory graphs: a formal framework to simulate biological regulatory networks with simple automata[END_REF]. One can also separate inhibitors from activators [START_REF] Bernot | Semantics of biological regulatory networks[END_REF] to increase the expressivity of the approach, or take into account time delays [START_REF] Thomas | Multistationarity, the basis of cell differentiation and memory. II. Logical analysis of regulatory networks in terms of feedback circuits[END_REF] between the beginning of the activation order and the synthesis of the product and conversely for the turn-off delays.

• Automatic generation of experiment schema from suitable models. In order to reduce again the set of suitable models, we would like to propose the biologist to perform an determining experiment. The result is then confronted to each model and only those which are coherent with the experiment, have to be kept. An experiment often consists to put the system in a particular state (partially specified) and to observe after a while if one or several gene products are present or not. This implies to extract the specificities of the biological application domain in order to define patterns of formulae expressing feasible experiments. • The modelling of a regulatory network concerns generally only a small part of the global regulatory network of the cell. It becomes crucial to prove that the dynamical properties of this sub-network are preserved when it is embedded into the global network. This is correlated to the treatment of knock-out mutants, identification of functional patterns [START_REF] Shen-Orr | Network motifs in the transcriptional regulation network of Escherichia coli[END_REF] as well as the structure of huge regulatory networks.

To achieve such development several directions have to be considered. High-level Petri nets are graphical oriented languages for design, specification, simulation and verification of systems. They are in particular well-suited for systems in which communication, synchronization and resource sharing are important. Clearly, biological systems present these characteristics, and modelling by such nets would allows us to take advantage of all results and tools in the field of high-level Petri nets. Hybrid automata can take into account the continuous aspects of a regulatory network: it is a mathematical model for hybrid systems, which combines, in a single formalism, automaton transitions for capturing discrete changes with differential equations for capturing continuous changes. Symbolic model checkers, as HyTech [START_REF] Henzinger | HyTech: A model checker for hybrid systems[END_REF], have been developed for the subclass of linear hybrid automata. It becomes possible to perform parametric analysis, i.e. to determine the values of parameters for which a linear hybrid automaton satisfies a temporallogic requirement.

These research perspectives aim to link modelling and experiments together, by furnishing to biologists model structuring methods and model validation tools from current researches in theoretical computer science. The resulting formal models are not only a posteriori explanations of biological results, they are guides for biological experiments whose success will be in fine the discriminating criterion.

Figure 2 :

 2 Figure 2: (a) Sigmoid relations between the concentration of u and the rate of synthesis of v and itself. As u is an activator of v and itself, see figure 1, both sigmoids are increasing. (b) Resulting qualitative levels of u.

Figure 3 :

 3 Figure 3: (a) Biological regulatory graph deduced from the genetic regulatory network of figure 1. (b) States of the biological regulatory graph.

Figure 3 -

 3 (b) shows the possible states of the biological regulatory graph of figure 3-(a).

Figure 4 :

 4 Figure 4: Time delays. Gene v has a unique regulator u which is an activator. Initially, both u and v are absent. Then, protein u appears and stimulates the expression of v (K v,ω v (q u ) = K v,{u} = 1). The resulting

Figure 5 :

 5 Figure 5: Synchronous and asynchronous dynamics for the model, given in table 1, of the biological regulatory graph of figure 3.

Figure 6 :

 6 Figure 6: Example of UDS of the biological regulatory graph of figure 3. Parameters are: k u = 20, k vu = 35, k uu = 40, θ vu = 10, θ uu = 20, λ u = 5 for the first equation and k v = 0, k uv = 25, θ uv = 16, and λ v = 2 for the second. Notice that θ uv < θ uu since 1 = t uv < t uu = 2.

Figure 7 : 2 =

 72 Figure 7: Domains of the phase space Ω of the UDS of figure 6 (θ u 1 = θ uv and θ u 2

Figure 8 :

 8 Figure 8: (a) The asynchronous state graph of the model M(G) of Table 1. (b) A trajectory of an UDS of M(G). The dotted arrows represent the extensions of solutions towards the attractors.

Figure 11 :

 11 Figure 11: Semantics of temporal connectives of CTL.

  qUr], then 1. L(s) = L(s) ∪ {E[qUr] } for all s such that r∈L(s), 2. Repeat: L(s) = L(s) ∪ {E[qUr]} for all states s such that q∈L(s) and which have a successor labelled with E[qUr], until there is no change.

Figure 12 :

 12 Figure 12: The life cycle of bacteriophage lambda.

Figure 13 :

 13 Figure 13: Biological regulatory graph G for immunity control.

Figure 14 :

 14 Figure 14: Likely paths from the initial state to the lytic and immune states (in bold). The dotted arrow is absent for the 44 models such that K cro,{cro,cI} = 3, M(G) included, whereas the dashed ones are absent for others.

Figure 15 :

 15 Figure 15: Dynamics of the mutant λ cIcro -obtained from M(G).

Table of Contents

 of • 1 Introduction • 2 Qualitative dynamics of biological regulatory networks o 2.1 Biological regulatory graphs o 2.2 Models of biological regulatory graphs o 2.3 Dynamics of models • 3 Differential modelling o 3.1 Ordinary differential equation systems o 3.2 Discretization map and domains o 3.3 Dynamics of differential equation systems o 3.4 Discrete and differential modelling are coherent

	o 3.5 Feedback circuit functionality
	3.5.1 Differential inclusion systems
	3.5.2 Steadiness of singular states
	3.5.3 Circuit characteristic states
	3.5.4 Circuit functionality

• 4 Formal methods o 4.1 Temporal logic o 4.2 Model checking o 4.3 A tool for the selection of models: SMBioNet • 5 Immunity control in bacteriophage lambda o 5.1 Biological regulatory graph o 5.2 Temporal properties o 5.3 Selected models o 5.4 Validation of models • 6 Conclusion

  of v tends to decrease. Different values for these parameters are possible and we call model of a biological regulatory graph, a possible instantiation of these parameters.

	Definition 4 [Model of a biological regulatory graph] Let G=(V,E) be a regulatory graph. A
	model of G, denoted M(G) by abuse of notation, is a family of natural numbers K v,ω indexed
	by the set of couples (v,ω) such that
	•

Table 1 :

 1 One possible model for the biological regulatory graph of figure3. This number increases exponentially with the number of predecessors of each variable and even if static constraints on parameters are used, as the Snoussi's constraints, it remains huge. Moreover, since parameters K are most often not measurable in vivo, additional properties deduced from biological experiments are needed to eliminate the models whose dynamics do not satisfy them.

	The table gives for each state the corresponding attractors and tendencies
	deduced from the model.
	Models M(G) of the biological regulatory graph G of figure 3 are all possible instantiations of
	six parameters: K u,{} , K u,{v} , K v,{} , K v,{u} , K v,{v} , K v,{u,v} . Because b u = 1 (resp. b v = 2) each
	K u,… (resp. K v,… ) can take 2 (resp. 3) different values. So 2 2 × 3 4 =324 different models can be
	a priori associated to G, but only 60 of them respect the Snoussi's constraints. Table 1 gives
	the tendencies of variables resulting from such a model.
	More generally there are ∏ v∈V (b v +1) 2 |G-(v)|	models associated to a biological regulatory
	graph G.

Selection of Models for Biological Networks, see http://smbionet.lami.univ-evry.fr

Presented in the next section.

Acknowledgements

The authors thank genopole ® -research in Evry (H. Pollard and P. Tambourin) for constant supports. We gratefully acknowledge the members of the genopole ® working groups observability and G 3 for stimulating interactions.