Ludovic Mailleret
email: ludovic.mailleret@antibes.inra.fr

Frédéric Grognard
email: frederic.grognard@sophia.inria.fr

Optimal Release Policy for Prophylactic Biological Control

In this paper, we study the prophylactic biological control strategy for greenhouse crops protection. The method consists in the preventive installation of natural enemies to fight against an invading pest, using discrete augmentative (i.e. inondative) releases of the natural enemies. We consider a simple non negative prey (pest) -predator (natural enemy) model in ordinary differential equations together with discrete augmentation of the predator population at constant frequency. Assuming we have a fixed budget to spend in natural enemies releases per time unit, we show the stability and efficiency of the prophylactic biological control strategy (i.e. the pest is eradicated) if this budget is larger than some value. Then we show that the optimal strategy to minimize worst case damage is to use the most frequent (and thus smallest) releases.

Introduction

Due to the high costs induced by greenhouse crops production, a particular attention must be given to reduce to the lowest the damage caused by pests. Within this context, prophylatic methods (i.e. prevention), able to fight pests at the time of their introduction (not after some delay) are more effective than curative ones. Till now, these methods are most widely used for spraying of chemicals; however, due to social (growers and consumers health), environmental (water pollution) and ecological (development of resistance by pests) considerations, chemical pesticides should be used with great care [START_REF] Jarvis | Managing Diseases in Greenhouse Crops[END_REF]. An alternative strategy is to use biological control, i.e. to introduce natural enemies to fight against pests, instead of chemicals.

Assuming that crop damage are low and then that crop is not limiting for pest growth, we consider a prey (pest) -predator (natural enemy) system modeled with two non negative ordinary differential equations. We consider that the crop never ends and study on this horizon the augmentative biological control method: the increase, at a fixed time period and with a fixed number, of the predator population. We consider that we have a constant budget (i.e. number of natural enemies) to spend for biological control per time unit.

We show the existence and stability of a periodic pest-eradicated solution (i.e. prophylactic solution) iff the budget is larger than a value determined by the model parameters. Moreover, as the prophylactic solution runs, the natural enemies are able to eradicate any level of unforeseen pest infestation. It is to be noted that this first part uses a method quite similar to the one proposed by [START_REF] Shulgin | Pulse vaccination stategy in the SIR endemic model[END_REF] in the context of mathematical epidemiology and rephrased by [START_REF] Liu | The dynamical behaviors of a Lotka-Volterra predator-prey model concerning integrated pest management[END_REF] in the context of mathematical modelling of crop protection.

We then suppose that the prophylactic solution runs while an infestation occurs and look for the releasing period that minimizes the worst case (according to the moment of infestation) crop damage due to the pest. We show that most frequent releases is the more efficient strategy with respect to uncertainties on the invading pest population level and on the model parameters.

Model Description and Analysis

Model Description

Since we consider the crop as non-limiting for pest growth, we simplify the tri-trophic ecosystem "crop -pest -natural enemies" as a bi-trophic one with only a prey (pest) and a predator population (natural enemies). According to classical population dynamics, we have the following non negative model:

ẋ = f (x) -g(x)y ẏ = kg(x)y -my
x denoting the preys and y the predators. f (x) denotes the growth speed of the preys, g(x) the predation rate, k the yield associated with the predation for the predators and m the mortality rate of the predators. Since biological processes are always difficult to model, we only assume weak qualitative hypotheses on the (smooth enough) functions f (.) and g(.):

Hypotheses 1 • f (0) = 0, f ′ (0) > 0 and ∀x ≥ 0, f ′′ (x) ≤ 0 • g(0) = 0, g ′ (0) > 0 and ∀x ≥ 0, g ′′ (x) ≥ 0
We now model the periodic release of predators; let us suppose that we have a fixed budget (i.e. fixed number of natural enemies) per time unit to spend for crop protection with prophylactic biological control. Let µ be the number of predators that can be released per time unit, then at each time period T , µT predators are added to the predator population y; we have the following non negative system of impulsive differential equations [START_REF] Bainov | Impulsive Differential Equations[END_REF]:

   ẋ = f (x) -g(x)y ẏ = kg(x)y -my ∀n ∈ N, y(nT + ) = y(nT ) + µT (1)

Model Analysis

Theorem 1. Under Hypotheses 1, model (1) possesses a prophylactic periodic solution (corresponding to the pest eradication):

(x p (t), y p (t)) = 0, µT e -m(t mod T ) 1 -e -mT (2)
which is (locally) stable iff:

µ > mf ′ (0) g ′ (0) (3)
Moreover, as the prophylactic solution runs and (3) holds, the natural enemies are able to eradicate any level of unforeseen pest infestation and there exists a release period T min under which the pest population is always decreasing.

Proof. We first focus on the existence and stability of the prophylactic solution i.e. on the dynamics of y(t) as no prey is present; we get:

   ẋ = 0 ẏ = -my ∀n ∈ N, y(nT + ) = y(nT ) + µT (4)
that yields: y((n + 1)T + ) = y(nT + )e -mT + µT . It is clear that the sequence (y(nT + )) n∈N has a single and globally stable equilibrium y ⋆ = µT /(1-e -mT ).

Then, system (4) possesses the single globally stable periodic solution [START_REF] Chicone | Ordinary Differential Equations with Applications[END_REF]. We now prove the (local) stability of the prophylactic solution for system [START_REF] Bainov | Impulsive Differential Equations[END_REF], that is to say with the presence of a prey population. We must then consider small perturbations around the periodic solution (x p (t), y p (t)), i.e. we consider (x(t), y(t)) = (x p (t), y p (t)) + (x(t), ỹ(t)). We get:

ẋ = f (x) -g(x)(ỹ + y p (t)) ẏ = kg(x)(y p (t) + ỹ) -mỹ (5)
At first order in x and ỹ, we have:

ẋ = (f ′ (0) -g ′ (0)y p (t))x ẏ = kg ′ (0)y p (t)x -mỹ (6)
which is a linear system (in x and ỹ) with periodic coefficients. We then use Floquet's theory which ensures the local stability of (0, 0) for system (6) (and thus of (x p (t), y p (t)) for ( 1)) iff the Floquet multipliers of the system have their absolute values lower than one (see e.g. [START_REF] Chicone | Ordinary Differential Equations with Applications[END_REF] for the theory). In our case, the two Floquet multipliers are e -mT ∈ (0, 1) and e

R T 0 (f ′ (0)-g ′ (0)yp(τ ))dτ > 0 that is lower than one iff µ > mf ′ (0)/g ′ (0). Then (x p (t), y p (t)) is locally stable iff (3) holds.
We now come back to system [START_REF] Shulgin | Pulse vaccination stategy in the SIR endemic model[END_REF]. From Hypotheses 1 and since y p (t) ≥ 0, system (5) is non-negative. We now consider a forward trajectory of system (5) initiated at time t 0 ≥ 0 at (x 0 , 0) ≥ 0. From ỹ ≥ 0 and Hypotheses 1, we have:

ẋ

≤ (f ′ (0) -g ′ (0)y p (t))x (7)
Thus, since y p (t) is T -periodic:

x(t) ≤ x 0 e R lT t 0 (f ′ (0)-g ′ (0)yp(τ ))dτ e (k-l) R T 0 (f ′ (0)-g ′ (0)yp(τ ))dτ e R t kT (f ′ (0)-g ′ (0)yp(τ ))dτ
with k the integer part of t/T and l the integer part of t 0 /T plus one. Suppose (3) holds, then T 0 (f ′ (0)-g ′ (0)y p (τ ))dτ < 0. Since k goes to infinity as t does, x(t) converges asymptotically to 0.

We come back to equation ( 7), we clearly have:

ẋ ≤ (f ′ (0) -g ′ (0) min t y p (t))x (8)
y p reaches its minimum at t = T -mod T and is equal to µT /(e mT -1) which is a decreasing function of T on (0, +∞) from +∞ to 0. Let us define T min > 0 such that min y p (T min ) = f ′ (0)/g ′ (0); then ∀T < T min , the right hand side of (8) is negative and x(t) is a decreasing function of t. ⊓ ⊔

3 Optimization of the Release Policy

Statement of the Problem

From now on, we consider that the prophylactic solution (2) runs as a pest population x 0 invades the crop at a time t 0 ∈ [0, T ) (without loss of generality). We consider that the invading population x 0 is small so that the linear system ( 6) is a good approximation of system (1). Moreover we assume that we have chosen T so that the pest population is always decreasing i.e. T ∈ (0, T min ).

To evaluate crop damage due to the pest, we use the concept of "Economic Injury Level" (EIL) that has been introduced from the early bases of theoretical biological control [START_REF] Stern | The integrated control concept[END_REF]. EIL (denoted x in the sequel) is defined as the lowest (positive) pest population level that will cause economic losses on the crop. For a fixed x 0 > x, we claim that the lower the time the pest population is above x the lower are the crop damage. Then we look for the release period T < T min that minimizes the time spent by x above x (denoted Π in the following) for its worst t 0 .

Main Result

We first look for the t 0 that maximizes the damage time Π, we have: Lemma 1. Suppose x 0 > x and T ∈ (0, T min ) are fixed, then one of the following holds:

i-∃k ∈ N, x 0 exp ((f ′ (0) -µg ′ (0) m )kT ) = x then Π = kT ii-∃k ∈ N, max t0 Π(t 0 ) = Π(t * 0 ) = (k + 1)T -t * 0 (i.e. x((k + 1)T ) = x)
We will first show case (i-) and prove that, otherwise, either case (ii-) holds or Π is maximum at t * 0 = 0, this latter case being impossible.

Proof. Suppose that Π(t 0 ) is maximum for t * 0 ∈ (0, T ) and such that Π(t * 0 ) ∈ (kT, (k + 1)T ) for some integer k. Pick t m 0 < t * 0 < t M 0 in (0, T ) such that Π(t 0 ) ∈ [kT, (k + 1)T ) for all t 0 ∈ [t m 0 , t M 0 ]. Within this set, define τ (t 0 ) as:

Π(t 0 ) = kT -t 0 + τ (t 0 ) (9)
Integrating ( 6) between t 0 and t 0 + Π(t 0 ), we get (with τ standing for τ (t 0 )):

x(t 0 +Π(t 0 )) = x 0 e f ′ (0)(kT -t0+τ ) e " g ′ (0)µT m(1-e -mT )
(-e -mt 0 +k(e -mT -1)+e -mτ ) " (10) which is, from the definition of Π, equal to x.

To have a maximum of Π at t * 0 , we need dΠ dt0 (t * 0 ) = 0. Thus from (9) we must have dτ dt0 (t * 0 ) = 1. We differentiate (10) with respect to t 0 and get:

dτ dt 0 (t 0 ) = g ′ (0)µT 1-e -mT e -mt0 -f ′ (0) g ′ (0)µT 1-e -mT e -mτ (t0) -f ′ (0) (11)
Clearly we must have τ (t * 0 ) = t * 0 . Thus from (10) there must exist an integer k such that:

x 0 e (f ′ (0)-µg ′ (0) m )kT = x and Π = kT does not depend on t 0 (case (i-)). Otherwise Π has no extremum at t * 0 such that t * 0 ∈ (0, T ) and Π(t * 0 ) ∈ (kT, (k + 1)T ). Two cases remain to be studied, either Π is maximum for t * 0 = 0 or t * 0 is such that t * 0 + Π(t * 0 ) = (k + 1)T ; t * 0 = T or t * 0 + Π(t * 0 ) = kT might be studied by k reparametrization.

Assume t * 0 = 0. Then from (11) and (9) the right derivative of Π at t * 0 is:

dΠ dt 0 (t * 0 = 0 + ) = g ′ (0)µT 1-e -mT -f ′ (0) g ′ (0)µT 1-e -mT e -mτ (t * 0 ) -f ′ (0) - 1
Since T < T min , we have from (8) f ′ (0) < g ′ (0)µT e -mT /(1 -e -mT ). Then both the numerator and denominator are positive, the former being larger than the latter. dΠ dt0 (t * 0 = 0 + ) is then positive and Π is minimum at t 0 = 0. Assume now that t * 0 is such that t * 0 + Π(t * 0 ) = (k + 1)T . Then from (9) τ (t * 0 ) = T and the left derivative of Π at t * 0 is:

dΠ dt 0 (t * - 0 ) = g ′ (0)µT 1-e -mT e -mt * 0 -f ′ (0) g ′ (0)µT 1-e -mT e -mT -f ′ (0) - 1
In a very same way as in the previous case, since T < T M we show that this derivative is positive. Similarly one can show that the right derivative of Π at t * 0 is negative. This can be performed through the reparametrization of k as k + 1 while noticing that this corresponds to τ (t * 0 ) = 0. Then for the considered x 0 and T , Π is maximum for t * 0 s.t.:

t * 0 + Π(t * 0 ) = (k + 1)T . ⊓ ⊔
Lemma 1 is quite natural. Indeed, the decrease that takes place between times T and kT is independent of t 0 . The worst case should then contain the end of the first-time interval where the predators are scarce rather than the beginning of the last interval where they are abundant.

We now focus on the evolution of the maximum of Π according to t 0 as the release period T varies in (0, T min ).

Theorem 2. Suppose x 0 is fixed. Let:

T 1 = m(ln (x 0 ) -ln (x)) µg ′ (0) -mf ′ (0) (12)
Then the following hold: i-there exists a n 0 ∈ N * such that min T <Tmin max t0 Π = T 1 is reached at

T = T n = T1 n for all integer n > n 0 . ii-∀n ≥ n 0 , ∆Π n = max T ∈[Tn+1,Tn] max t0 Π -min T max t0 Π is a decreas- ing function of n.
Moreover it tends to 0 as n tends to infinity.

Proof. Consider T 1 ; it is clear that there exists an integer n 0 such that for all integer n > n 0 , T n < T min . Now consider a T n with n > n 0 ; using (12), we have:

x

0 e (f ′ (0)-µg ′ (0) m )nTn = x 0 e (f ′ (0)-µg ′ (0) m )T1 = x
Such a T n corresponds to case (i-) in Lemma 1, thus Π = T 1 . Now we show that Π = T 1 is the minimum (according to T ) of max t0 Π. Suppose that we have T = T n (∀n ≥ n 0 ). According to case (ii-) in Lemma 1, ∃k ∈ N such that max t0 Π = (k + 1)Tt * 0 and x((k + 1)T ) = x. Using (10) we have:

x 0 e f ′ (0)((k+1)T -t * 0 ) e g ′ (0)µT m " 1-e -mt * 0
1-e -mT +(k+1)

« = x
Which, with a little effort, yields:

(k + 1)T -t * 0 = T 1 + µg ′ (0) µg ′ (0) -mf ′ (0) 1 -e -mt * 0 1 -e -mT T -t * 0 (13)
From (3), µg ′ (0) -mf ′ (0) > 0 and since t * 0 < T it can easily be shown that:

1 -e -mt * 0 1 -e -mT T -t * 0 ≥ 0
Then, (k + 1)Tt * 0 = max t0 Π is larger than T 1 which gives part (i-) of Theorem 2. Now we prove part (ii-). From (i-) and (13) we have for the corresponding t * 0 ∈ [0, T n ]:

∆Π n = max T ∈[Tn+1,Tn] max t0 Π -min T max t0 Π = µg ′ (0) µg ′ (0)-mf ′ (0) max T ∈[Tn+1,Tn] 1-e -mt * 0 1-e -mT T -t * 0 (14)
We now look for disjoint sets, depending on n, that contain ∆Π n . We first notice that 1-e -mt * 0

1-e -mT T belongs to 1-e -mt * 0

1-e -mT n+1 T n+1 , 1-e -mt * 0

1-e -mTn T n for all T ∈ [T n+1 , T n ] and with t * 0 , defined by case (ii-) of Lemma 1 for k = n, depending on T . Since (n + 1)Tt * 0 (T ) = T 1 (case (i-)) at T = T n+1 and T = T n , we have t * 0 (T n+1 ) = 0 and t * 0 (T n ) = T n . Then, for all s ∈ [0, T n ], ∃T ∈ [T n+1 , T n ] such that t * 0 (T ) = s. Thus:

max T ∈[Tn+1,Tn] 1-e -mt * 0 1-e -mT T -t * 0 ≤ max T ∈[Tn+1,Tn] 1-e -mt * 0 1-e -mTn T n -t * 0 = max t * 0 ∈[0,Tn] 1-e -mt * 0 1-e -mTn T n -t * 0 Differentiating 1-e -mt * 0
1-e -mTn T n -t * 0 with respect to t * 0 , we show that it reaches its maximum for:

t * 0 = 1 m ln mT n 1 -e -mTn Then: max T ∈[Tn+1,Tn] 1-e -mt * 0 1-e -mT T -t * 0 ≤ 1 m mTn 1-e -mTn -1 -ln mTn 1-e -mTn
1-e -mT T -t * 0 ≥ 1-e -mt * 0 ( T ) 1-e -m T T -t * 0 ( T ) ≥ 1-e -mt * 0 ( T ) 1-e -mT n+1 T n+1 -t * 0 ( T ) = H(T n+1 )
We have shown that:

H(T n+1 ) ≤ max T ∈[Tn+1,Tn] 1 -e -mt * 0 1 -e -mT T -t * 0 ≤ H(T n )
This, together with (14), implies that ∆Π n is a decreasing function of n. Moreover H(0) = 0 and T n tends to 0 as n tends to infinity. Then ∆Π n tends to 0 as n tends to infinity. ⊓ ⊔

Discussion

Our main result directly comes from the interpretation of Theorem 2. We have found the values T n of the release period T that solve our min max problem through point (i-). However we use point (ii-) and a robustness argument to chose a particular T n .

Indeed, we would like the optimality property to be robust to the level of invading pest population x 0 (that we do not know in advance), as well as to the model parameters. Our choice of T = T n will be made on the basis of the assumed nominal values of these parameters, but the uncertainties will make that our choice will not be a true T n . Most probably our choice of T will lie strictly between two true values of T k and then Π would not be minimum. Therefore we should pick T so that, in its neighborhood, Π(T ) does not attain too large values. We conclude from point (ii-) of Theorem 2 that we should then choose T small (through n large).

To summarize our results, provided the budget spent by time unit is large enough, prophylactic biological control is able to protect crops from invading pests. Moreover the most frequent releases of natural enemies minimize the worst case damage caused by pests on the crop, robustly to the initial pest population level and to model parameters.

  -mT n+1 T n+1 -t * 0 = max t * 0 ∈[0,Tn] 1-e -mt * 0 1-e -mT n+1 T n+1 -t * 0 = H(T n+1 ) This maximum is reached for T = T ∈ [T n+1 , T n ]. Then: max T 1-e -mt * 0

Fig. 1 .

 1 Fig.1. Graphical illustration of Theorem 2: Π as a function of the release period T ∈ (0, Tmin) for some n > n0 and with γ = µg ′ (0)/(µg ′ (0) -mf ′ (0)).


Frédéric Grognard
email: frederic.grognard@inria.fr

Hidde De Jong
email: hidde.de-jong@inrialpes.fr

Jean-Luc Gouzé
email: gouze@sophia.inria.fr

Piecewise-linear models of genetic regulatory networks: theory and example

Summary. The experimental study of genetic regulatory networks has made tremendous progress in recent years resulting in a huge amount of data on the molecular interactions in model organisms. It is therefore not possible anymore to intuitively understand how the genes and interactions together influence the behavior of the system. In order to answer such questions, a rigorous modeling and analysis approach is necessary. In this chapter, we present a family of such models and analysis methods enabling us to better understand the dynamics of genetic regulatory networks. We apply such methods to the network that underlies the nutritional stress response of the bacterium E. coli.

The functioning and development of living organisms is controlled by large and complex networks of genes, proteins, small molecules, and their interactions, socalled genetic regulatory networks. The study of these networks has recently taken a qualitative leap through the use of modern genomic techniques that allow for the simultaneous measurement of the expression levels of all genes of an organism. This has resulted in an ever growing description of the interactions in the studied genetic regulatory networks. However, it is necessary to go beyond the simple description of the interactions in order to understand the behavior of these networks and their relation with the actual functioning of the organism. Since the networks under study are usually very large, an intuitive approach for their understanding is out of question. In order to support this work, mathematical and computer tools are necessary: the unambiguous description of the phenomena that mathematical models provide allows for a detailed analysis of the behaviors at play, though they might not exactly represent the exact behavior of the networks.

In this chapter, we will be mostly interested in the modeling of the genetic regulatory networks by means of differential equations. This classical approach allows precise numerical predictions of deterministic dynamic properties of genetic regulatory networks to be made. However, for most networks of biological interest the application of differential equations is far from straightforward. First, the biochemical reaction mechanisms underlying the interactions are usually not or incompletely known, which complicates the formulation of the models. Second, quantitative data on kinetic parameters and molecular concentrations is generally absent, even for extensively-studied systems, which makes standard numerical methods difficult to apply. In practice, the modeler disposes of much weaker information on the network components and their interactions. Instead of details on the mechanisms through which a protein regulates a gene, we typically only know whether the protein is an activator or an inhibitor. And even if it had been shown, for example, that the protein binds to one or several sites upstream of the coding region of the gene, numerical values of dissociation constants and other parameters are rarely available. At best, it is possible to infer that the regulatory protein strongly or weakly binds to the DNA, with a greater affinity for one site than for another.

Due to those uncertainties, we cannot hope to build a model that is guaranteed to reproduce the exact behavior of the considered genetic regulatory network. No model will be quantitatively accurate. It is therefore necessary to concentrate on the construction of models that reproduce the qualitative dynamical properties of the network, that is, dynamical properties that are invariant for a range of parameter values and reaction mechanisms. The qualitative properties express the intimate connection between the behavior of the system and the structure of the network of molecular interactions, independently from the quantitative details of the latter.

Consequently, qualitative approaches have been developed for the modeling, analysis, and simulation of genetic regulatory networks and other networks of biological interactions: Boolean networks [START_REF] Kauffman | The Origins of Order: Self-Organization and Selection in Evolution[END_REF][START_REF] Thomas | Biological Feedback[END_REF], Petri nets [START_REF] Koch | Application of Petri net theory for modelling and validation of the sucrose breakdown pathway in the potato tuber[END_REF][START_REF] Reddy | Qualitative analysis of biochemical reaction systems[END_REF], process algebras [START_REF] Regev | Representation and simulation of biochemical processes using the π-calculus process algebra[END_REF], qualitative differential equations [START_REF] Heidtke | Design and implementation of a qualitative simulation model of λ phage infection[END_REF], hybrid automata [START_REF] Ghosh | Symbolic reachable set computation of piecewise affine hybrid automata and its application to biological modelling: Delta-Notch protein signalling[END_REF],... In this chapter, we concentrate on one particular class of qualitative models of genetic regulatory networks, originally proposed by Glass and Kauffman [START_REF] Glass | The logical analysis of continuous non-linear biochemical control networks[END_REF]: piecewise-linear (PL) differential equations. In Section 1, we describe this family of models and give a small example. In Section 2, we show qualitative results that have been obtained for the analysis of such systems. We then illustrate these models on the nutritional stress response of E. coli in Section 3, before discussing remaining challenges for the analysis and control of such models in Section 4.

Models of genetic regulatory networks

Among the many emerging families of models (see [START_REF] Jong | Modeling and simulation of genetic regulatory systems: a literature review[END_REF]), a class of piecewise-linear (PL) models, originally proposed by Glass and Kauffman [START_REF] Glass | The logical analysis of continuous non-linear biochemical control networks[END_REF], has been widely used in modeling genetic regulatory networks. The variables in the piecewise-linear differential equation (PLDE) models are the concentrations of proteins encoded by the genes, while the differential equations describe the regulatory interactions in the network by means of step functions. The use of step functions is motivated by the switch-like behavior of many of the interactions in genetic regulatory networks [START_REF] Ptashne | A Genetic Switch: Phage λ and Higher Organisms[END_REF], but it leads to some mathematical difficulties. The vector field for the PLDE model is undefined when one of the variables assumes a value where the step function is discontinuous, referred to as a threshold value. Recent work by Gouzé and Sari [START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF] uses an approach due to Filippov to define the solutions on the threshold hyperplanes.

The approach involves extending the PLDE to a piecewise-linear differential inclusion (PLDI). As is well known, such discontinuities can lead to sliding modes. The definitions and results of this section are mainly taken from [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF].

The family of PL-models is best illustrated with an example: the schematic diagram in Figure 1 describes a simple genetic regulatory network. In this example, the genes a and b code for the proteins A and B, which in turn control the expression of the two genes a and b. Protein A inhibits gene a and activates gene b above certain threshold concentrations, which are assumed to be different. Similarly protein B inhibits gene b and activates gene a above different threshold concentrations. This two-gene regulatory network is simple but represents many features of regulation found in real networks: auto-regulation, cross-regulation and inhibition/activation. Such a two-gene network could be found as a module of a more complex genetic regulatory network from a real biological system. The equations modeling the example network in Figure 1 can be written down as

ẋa = κ a s + (x b , θ 1 b )s -(x a , θ 2 a ) -γ a x a ẋb = κ b s + (x a , θ 1 a )s -(x b , θ 2 b ) -γ b x b (1)
where s + (x s , θ s ) is equal to 0 when x s < θ s and equal to 1 when x s > θ s and s -(x s , θ s ) = 1s + (x s , θ s ). In this model, gene a is expressed at a rate κ a if the concentration x b of protein b is above the threshold θ 1 b and the concentration x a of protein A is below the threshold θ 2 a . Similarly, gene b is expressed at a rate κ b if the concentration x a of protein A is above the threshold θ 1 a and the concentration x b of the protein B is below the threshold θ 2 b . Degradation of both proteins is assumed to be proportional to their own concentrations, so that the expression of the genes a and b is modulated by the degradation rates γ a x a and γ b x b respectively. Such a model is readily generalized to models containing both expression and degradation terms for each gene:

ẋi = f i (x) -γ i x i
where f i (x) represents the expression rate of gene i, depending on the whole state x = (x 1 , • • • , x n ) T and γ i x i is the degradation rate. However, the expression rates of (1) have the additional property of being constant for values of x a and x b belonging to intervals that do not contain thresholds values θ j i . This can be rewritten by detailing f i (x) as follows:

f i (x) = Li l=1 κ il b il (x)
where b il (x) is a combination of step-functions s ± (x r , θ j r ) and κ il > 0 is a rate parameter. The generalized form of (1) is a piecewise linear model

ẋ = f (x) -γx (2)
where the model is linear within hyper-rectangles of the state-space.

The dynamics of the piecewise-linear system (2) can be studied in the ndimensional state-space

Ω = Ω 1 × Ω 2 × • • • × Ω n , where each Ω i is defined by Ω i = {x ∈ IR + | 0 ≤ x i ≤ max i } for some positive parameter max i > max x ∈ Ω fi(x) γi
. A protein encoded by a gene will be involved in different interactions at different concentration thresholds, so for each variable x i , we assume there are p i ordered thresholds θ 1 i , • • • , θ pi i (we also define θ 0 i = 0 and θ pi+1 i = max i ). The (n -1)-dimensional hyperplanes defined by these thresholds partition Ω into hyper-rectangular regions we call domains. Specifically, a domain D ⊂ Ω is defined to be a set

D = D 1 × • • • × D n ,
where D i is one of the following:

D i = {x i ∈ Ω i |0 ≤ x i < θ 1 i } D i = {x i ∈ Ω i |θ j i < x i < θ j+1 i } for j ∈ {1, • • • , p i -1} D i = {x i ∈ Ω i |θ pi i < x i ≤ max i } D i = {x i ∈ Ω i |x i = θ j i } for j ∈ {1, • • • , p i } A domain D ∈ D is

Solutions in regulatory domains

For any regulatory domain D ∈ D r , the function f (x) is constant for all x ∈ D, and it follows that the piecewise-linear system (2) can be written as a linear vector field

ẋ = f D -γx (3)
where In Figure 2, example (1) is used to illustrate this concept: the considered regula-

f D is constant in D. Restricted to D,
tory domain D is {x a ∈ Ω a |θ 1 a < x a < θ 2 a } × {x b ∈ Ω b |θ 1 b < x b < θ 2 b }, so that system (1) becomes ẋa = κ a -γ a x a ẋb = κ b -γ b x b
and the corresponding focal point is κa γa , κ b γ b . In the figure, this focal point is supposed to be outside of D: every solution starting in D therefore exits this domain in finite time.

Solutions in switching domains

In switching domains, the PL system (2) is not defined, since in a switching domain of order k ≥ 1, k variables assume a threshold value. If solutions do not simply go through a switching domain, it is necessary to give a definition of what a solution can be on that domain. Classically, this is done by using a construction originally proposed by Filippov [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF] and recently applied to PL systems of this form [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].

The method consists of extending the system (3) to a differential inclusion,

ẋ ∈ H(x), (4)
where H is a set-valued function (i.e. H(x) ⊆ IR n ). If D is a regulatory domain, then we define H simply as

H(x) = {f D -γx}, (5)
for x ∈ D. If D is a switching domain, for x ∈ D, we define H(x) as

H(x) = co({f D ′ -γx | D ′ ∈ R(D)}), (6)
where

R(D) = {D ′ ∈ D r |D ⊆ ∂D ′ }
is the set of all regulatory domains with D in their boundary, and co(X) is the closed convex hull of X. For switching domains, H(x) is generally multi-valued so we define solutions of the differential inclusion as follows.

Definition 3. A solution of ( 4) on [0, T ] in the sense of Filippov is an absolutely continuous function (w.r.

t. t) ξ t (x 0 ) such that ξ 0 (x 0 ) = x 0 and ξt ∈ H(ξ t ), for almost all t ∈ [0, T ].
In order to more easily define these Filippov solutions, it is useful to define a concept analogous to the focal points defined for regulatory domains, extended to deal with switching domains.

Definition 4. Let D ∈ D s be a switching domain of order k. Then its focal set Φ(D) is Φ(D) = supp(D) ∩ co({φ(D ′ ) | D ′ ∈ R(D)}). (7)
Hence Φ(D) for D ∈ D s is the convex hull of the focal points φ(D ′ ) of all the regulatory domains D ′ having D in their boundary, as defined above, intersected with the threshold hyperplane supp(D) containing the switching domain D (Figure 3). We have shown that

H(x) = γ(Φ(D) -x) (8)
which is a compact way of writing that H(x) = {y ∈ IR n | ∃φ ∈ Φ(D) such that y = γ(φx)}. The Filippov vector field is defined by means of the focal set.

If Φ(D) = { }, with D a switching domain, solutions will simply cross D; otherwise, sliding mode is possible and convergence takes place "in the direction" of Φ(D). If Φ(D)∩D = { }, solutions eventually leave D. In the case where Φ(D)∩D is not empty, it can be assimilated to an equilibrium set within D towards which all solutions will converge in the following sense ) is simply the segment that links them, so that [START_REF] Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF] implies that φ(D) is the intersection of this segment with supp(D).

D 1 D 2 φ(D 1 ) φ(D 2 ) D φ(D)
Lemma 1. [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF] For every regulatory domain D ∈ D r , all solutions ξ t in D monotonically converge towards the focal set Φ(D). For every switching domain D ∈ D s , the non-switching component (ξ t ) i of the solution ξ t in D monotonically converges towards the closed interval

π i (Φ(D)) = {φ i ∈ Ω i | φ ∈ Φ(D)}, the projection of Φ(D) onto Ω i , if (ξ 0 ) i ∈ π i (Φ(D)). Every switching component (ξ t ) i of the solution ξ t in D is a constant (ξ t ) i = π i (Φ(D)) = θ qi i .
Basically, this means that convergence does not take place towards Φ(D), but towards the smallest hyper-rectangle that contains Φ(D). Indeed, if Φ(D) is neither empty, nor a singleton, and ξ t0 belongs to Φ(D), the Filippov vector field at this point is defined as H(ξ t0 ) = γ(Φ(D)ξ t0 ) and there is no guarantee that no element of H(ξ t0 ) points outside of Φ(D) (we know however that a solution stays at ξ t0 ). Due to the structure of the differential equations, it is on the other hand certain that the transient solution does not leave the smallest hyper-rectangle containing Φ(D). This phenomenon is illustrated in Figure 4 We then have the following corollary

Corollary 1. [3] All solutions ξ t in D converge towards Π(D), if ξ 0 ∈ Π(D). For all solutions ξ t in D, Π(D) is invariant. φ1 φ2 Π(D) Φ(D)
Φ(D) ∩ supp(D ′ ) = {}, ∀D ′ ⊆ ∂D. (9)
it has been possible to develop stability results for this family of systems.

Stability and qualitative properties of PL models

The stability analysis of the various equilibria is a direct consequence of the analysis in the previous section. It is easily seen that equilibria xr in some D ∈ D r are asymptotically stable. Indeed, they are the focal points of the domains in which they are contained, so that the convergence that was described in the previous section, leads to asymptotic stability. The more difficult part consists in defining and handling the stability of Filippov equilibria that lie in switching surfaces. In a switching domain D ∈ D s , recall that solutions are defined by considering the differential inclusion H(x). We say that a point y ∈ Ω is an equilibrium point for the differential inclusion if 0 ∈ H(y), [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF] where H is computed using the Filippov construction in [START_REF] Jong | Genetic Network Analyzer: Qualitative simulation of genetic regulatory networks[END_REF]. In other words, there is a solution in the sense of Filippov, ξ t , such that ξ t (y) = y, ∀t > 0. We call such a point a singular equilibrium point. It is easily seen that, for y to be an equilibrium point inside D, it must belong to Φ(D). Also, since Assumption 1 prevents Φ(D) from intersecting the border of D, we then have that Φ(D) ⊂ D. Every element φ of Φ(D) is then an equilibrium when Φ(D) ⊂ D so that, for every φ ∈ Φ(D), there exists a solution ξ t (φ) = φ for all t.

One of the interesting results of [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF] concerns the link between the configuration of the state transition graph and the stability of an equilibrium. This discrete, qualitative description of the dynamics of the PL system that underlies the qualitative simulation of genetic regulatory networks was originally due to Glass. It indicates the passages between the different domains making up the phase space. A state transition graph is a directed graph whose vertices are the domains of the system and whose edges are the possible transitions between these domains (easily determined by examining the PL model [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF]). The transition graph of system (1) is illustrated in Figure 5. For a two-dimensional system, we show how this graph indicates the stability of singular equilibria: Theorem 1. [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF] Let the dimension of the PL model be 2, and let D be a switching domain containing a singular equilibrium point φ(D). If for all regulatory domains D ′ ∈ R(D) (that is, adjacent to D), there exists a transition from D ′ to D in the state transition graph, then φ(D) is asymptotically stable. This result is purely qualitative: the actual value of the parameters is not needed. It can be directly applied to show that the singular equilibrium (x a , x b ) = (θ 2 a , θ 2 b ), corresponding to D 19 on Figure 5, is asymptotically stable because there are transitions to D 19 from D 13 , D 15 , D 23 and D 25 , the regulatory domains adjacent to D 19 .

x b D 2 D 7 D 12 D 17 D 22 0 x a max b θ 1 a θ 2 a max a D 1 D 3 D 5 D 6 D 8 D 13
A generalization, but in a weaker form, of this theorem to dimension n is also available. These results are very helpful for the qualitative analysis of the genetic regulatory networks. However, some stable equilibria cannot be identified through those criteria. Some less restrictive criteria are therefore under development.

Besides this method, we can discover stable equilibria that would not have been directly identified by our criteria, through a rigorous simplification of the model. This can be done through model reduction or identification of regions of the statespace that cannot be reached by the solutions (maybe after some finite time). We will illustrate the kind of things that can be done on an example in the following section. In that section, since the resulting models are very simple, we do not need to go back to transition graph analysis at the end of the reduction procedure, but we could have done so, and it will be necessary to do so if the model reduction procedure does not yield very small models.

Carbon starvation response of E. coli

We will present a specific model reduction and stability analysis for the model in dimension 6 of the carbon starvation response of E. coli of Ropers et al. [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]. In their natural environment, bacteria like Escherichia coli rarely encounter conditions allowing continuous, balanced growth. While nutrients are available, E. coli cells grow quickly, leading to an exponential increase of their biomass, a state called exponential phase. However, upon depletion of an essential nutrient, the bacteria are no longer able to maintain fast growth rates, and the population consequently enters a non-growth state, called stationary phase (Figure 6). During the transition from exponential to stationary phase, each individual E. coli bacterium undergoes numerous physiological changes, concerning among other things the morphology and the metabolism of the cell, as well as gene expression [START_REF] Huisman | Morphological and physiological changes during stationary phase[END_REF]. These changes enable the cell to survive prolonged periods of starvation and be resistant to multiple stresses. This carbon starvation response can be reversed and growth resumed, as soon as carbon sources become available again.

? Exponential phase

Signal of nutrient deprivation

Stationary phase

Fig. 6. Nutrient-stress response of bacteria during the transition from exponential to stationary phase.

On the molecular level, the transition from exponential phase to stationary phase is controlled by a complex genetic regulatory network integrating various environmental signals [START_REF] Hengge-Aronis | The general stress response in Escherichia coli[END_REF][START_REF] Martinez-Antonio | Identifying global regulators in transcriptional regulatory networks in bacteria[END_REF][START_REF] Wick | Molecular components of physiological stress responses in Escherichia coli[END_REF]. The molecular basis of the adaptation of the growth of E. coli to carbon starvation conditions has been the focus of extensive studies for decades [START_REF] Hengge-Aronis | The general stress response in Escherichia coli[END_REF]. However, notwithstanding the enormous amount of information accumulated on the genes, proteins, and other molecules known to be involved in the stress adaptation process, there is currently no global understanding of how the response of the cell emerges from the network of molecular interactions. Moreover, with some exceptions [START_REF] Botsford | Cyclic AMP in prokaryotes[END_REF][START_REF] Harman | Allosteric regulation of the cAMP receptor protein[END_REF][START_REF] Wang | Modeling of inducer exclusion and catabolite repression based on a PTS-dependent sucrose and non-PTS-dependent glycerol transport systems in Escherichia coli K-12 and its experimental verification[END_REF], numerical values for the parameters characterizing the interactions and the molecular concentrations are absent, which makes it difficult to apply traditional methods for the dynamical modeling of genetic regulatory networks. The notation follows, in a somewhat simplified form, the graphical conventions proposed by Kohn [START_REF] Kohn | Molecular interaction maps as information organizers and simulation guides[END_REF]. The contents of the boxes labeled 'Activation' and 'Supercoiling' are detailed in [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF].

The above circumstances have motivated the qualitative analysis of the carbon starvation response network in E. coli [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]. The objective of the study was to simulate the response of an E. coli bacterium to the absence or presence of carbon sources in the growth medium. To this end, an initial, simple model of the carbon starvation response network has been built on the basis of literature data. It includes six genes that are believed to play a key role in the carbon starvation response (Figure 7). More specifically, the network includes genes encoding proteins whose activity depends on the transduction of the carbon starvation signal (the global regulator crp and the adenylate cyclase cya), genes involved in the metabolism (the global regulator fis), cellular growth (the rrn genes coding for stable RNAs), and DNA supercoiling, an important modulator of gene expression (the topoisomerase topA and the gyrase gy-rAB).

Model of carbon starvation response

The graphical representation of the network has been translated into a PL model supplemented with parameter inequality constraints. The resulting model consists of seven variables, one concentration variable for the product of each of the six genes ((x c , x y , x f , x g , x t , x r ) for (crp, cya, fis, gyrAB, topA, rrn)) and one input variable u s representing the presence or absence of a carbon starvation signal [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]. The 38 parameters are constrained by 54 parameter inequalities, the choice of which is largely determined by experimental data.

The model of Ropers et al. is:

                         ẋc = κ 1 c + κ 2 c s -(x f , θ 2 f )s + (x c , θ 1 c )s + (x y , θ 1 y )s + (u s , θ s ) +κ 3 c s -(x f , θ 1 f ) -γ c x c ẋy = κ 1 y + κ 2 y 1 -s + (x c , θ 3 c )s + (x y , θ 3 y )s + (u s , θ s ) -γ y x y ẋf = κ 1 f + κ 2 f s + (x g , θ 1 g )s -(x t , θ 2 t ) 1 -s + (x c , θ 1 c )s + (x y , θ 1 y )s + (u s , θ s ) s -(x f , θ 5 f ) -γ f x f ẋg = κ g 1 -s + (x g , θ 2 g )s -(x t , θ 1 t ) s -(x f , θ 4 f ) -γ g x g ẋt = κ t s + (x g , θ 2 g )s -(x t , θ 1 t )s + (x f , θ 4 f ) -γ t x t ẋr = κ 1 r s + (x f , θ 3 f ) + κ 2 r -γ r x r
with u s = 0 in the presence of carbon sources and u s = 1 in a depleted environment (and θ s = 0.5). In order to uniquely determine the situation of the various focal points in the state-space, the following constraints on the parameters are needed:

                                 0 < θ 1 c < θ 2 c < θ 3 c < max c , θ 1 c < κ 1 c γc < θ 2 c , θ 1 c < (κ 1 c +κ 2 c ) γc < θ 2 c , θ 3 c < (κ 1 c +κ 3 c ) γc < max c 0 < θ 1 y < θ 2 y < θ 3 y < max y , θ 1 y < κ 1 y γy < θ 2 y , θ 3 y < (κ 1 y +κ 2 y ) γy < max y 0 < θ 1 f < θ 2 f < θ 3 f < θ 4 f < θ 5 f < max f , θ 1 f < κ 1 f γ f < θ 2 f , θ 5 f < (κ 1 f +κ 2 f ) γ f < max f 0 < θ 1 g < θ 2 g < max g , θ 2 g < κg γg < max g 0 < θ 1 t < θ 2 t < max t , θ 2 t < κt γt < max t 0 < θ r < max r , 0 < κ 2 r γr < θ r , θ r < (κ 1 r +κ 2 r ) γr < max r
A qualitative analysis of this model has been carried out in [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF] by using GNA (Genetic Network Analyzer [START_REF] Jong | Genetic Network Analyzer: Qualitative simulation of genetic regulatory networks[END_REF]), a computer tool that automatically generates the statetransition graph and possible trajectories in that graph, that is, qualitative solutions that are possible for this system. The following simulations are produced for the transition to the stationary phase (Figure 8) and to the exponential phase (Figure 9). In the first case, we see that the solution converges towards a single region of the state space, where we can guess that convergence towards an equilibrium takes place. In the second case, the behavior of the solution is not as clear: oscillations can be detected between various regions but it is impossible to say, based on the transition graph alone, if those oscillations are damped or not. Therefore, it is useful to try and analyze the model further to check what kind of oscillations take place (and in the same time if convergence actually takes place towards an equilibrium in the case of the entry in stationary phase).

Asymptotic dynamics

Since the 6-dimensional model, with all its constraints, is too complex to handle directly, we first check if some kind of simplifications can be made. Independently of the case that we will study (stationary phase or exponential phase conditions), we notice that • x r is a variable whose evolution depends on, but does not influence the rest of the system. As a consequence, it can be removed from the analysis. Once the analysis of the remaining 5-dimensional system is completed, we will be able to easily identify the consequence of its behavior on the concentration of stable RNAs (x r ). • There exists a finite time after which x t (t) ≤ θ 1 t since, as long as x t > θ 1 t , the x t dynamics reduces to ẋt = -γ t x t .

Once x t reaches θ 1 t , we cannot a priori rule out a sliding mode along x t = θ 1 t . Since θ 1 t < θ 2 t , this indicates that we can replace s -(x t , θ 2 t ) with 1 for the purpose of our analysis. We simply consider that the aforementioned finite time has already occurred.

• Similar studies show that x c (t) ≥ θ 1 c and x y (t) ≥ θ 1 y after some finite time. We can then replace s + (x c , θ 1 c ) and s + (x y , θ 1 y ) with 1 in our analysis. The system that we need to analyze has now become

             ẋc = κ 1 c + κ 2 c s -(x f , θ 2 f )s + (u s , θ s ) + κ 3 c s -(x f , θ 1 f ) -γ c x c ẋy = κ 1 y + κ 2 y 1 -s + (x c , θ 3 c )s + (x y , θ 3 y )s + (u s , θ s ) -γ y x y ẋf = κ 1 f + κ 2 f s + (x g , θ 1 g ) s -(u s , θ s )s -(x f , θ 5 f ) -γ f x f ẋg = κ g 1 -s + (x g , θ 2 g )s -(x t , θ 1 t ) s -(x f , θ 4 f ) -γ g x g ẋt = κ t s + (x g , θ 2 g )s -(x t , θ 1 t )s + (x f , θ 4 f ) -γ t x t
The next simplification step consists in seeing that x y does not influence the rest of the model, so that it can be removed, and that x c does not influence the rest of the model either (except x y ) so that it can also be removed. These actions are in the same line of thought as the removal of x r . As a consequence of these simplifications, we are able to see that the core of the long term dynamics is not really influenced by x r , x y and x c . We now have the three-dimensional system: x cya

x crp u s x f is x topA x rrn 0 0 θ 5 f is θ 3 f is 0 max crp max cya max f is CRP Cya Fis 0 Signal max s θ s θ 3 crp θ 1 crp θ 2 crp θ 3 cya θ 2 cya θ 1 cya θ 4 f is θ 2 f is θ 1 f is 0 GyrAB θ 2
c , xy = θ 3 y , x f < θ 1 f , xg = θ 2 g , xt < θ
     ẋf = κ 1 f + κ 2 f s + (x g , θ 1 g ) s -(u s , θ s )s -(x f , θ 5 f ) -γ f x f ẋg = κ g 1 -s + (x g , θ 2 g )s -(x t , θ 1 t ) s -(x f , θ 4 f ) -γ g x g ẋt = κ t s + (x g , θ 2 g )s -(x t , θ 1 t )s + (x f , θ 4 f ) -γ t x t (11)
Once we have analyzed the behavior of the solutions of this model, we will be able to reconstruct what happens with x c , x y and x r . For this analysis, we still suppose that x t ≤ θ 1 t .

Asymptotic dynamics in the absence of carbon sources

The analysis of the case u s = 1, the stationary phase solution in a depleted environment, is very straightforward. System [START_REF] Ghosh | Symbolic reachable set computation of piecewise affine hybrid automata and its application to biological modelling: Delta-Notch protein signalling[END_REF] becomes

   ẋf = -γ f x f ẋg = κ g 1 -s + (x g , θ 2 g )s -(x t , θ 1 t ) s -(x f , θ 4 f ) -γ g x g ẋt = κ t s + (x g , θ 2 g )s -(x t , θ 1 t )s + (x f , θ 4 f ) -γ t x t
so that x f goes to 0. It is then directly seen that, after a finite time (the time taken for x f to fall below θ 4 f ), we have ẋt = -γ t x t so that x t also goes to zero. The x g dynamics then reduce to

ẋg = κ g s -(x g , θ 2
g )γ g x g so that x g reaches θ 2 g in finite time. The three dimensional system thus has a very simple behavior: the state goes to (x f , x g , x t ) = (0, θ 2 g , 0). Since the solutions of the 6-dimensional system are bounded, the behavior of the other three states can be deduced from the analysis of the corresponding equations with (x f , x g , x t ) approaching their equilibrium (so that x f < θ 1 f , x t < θ 1 t and x g > θ 1

g ). We then have:

   ẋc = κ 1 c + κ 2 c + κ 3 c -γ c x c ẋy = κ 1 y + κ 2 y 1 -s + (x c , θ 3 c )s + (x y , θ 3 y ) -γ y x y ẋr = κ 2 r -γ r x r
It is then directly seen that, once (x f , x g , x t ) is close to its equilibrium, the variables (x c , x r ) exponentially converge towards (

κ 1 c +κ 2 c +κ 3 c γc , κ 2
r γr ) while x y reaches θ 3 y in finite time.

Asymptotic dynamics in the presence of carbon sources

The case u s = 0, the behavior of the model in an environment rich in carbon sources, is more intricate to analyze. System [START_REF] Ghosh | Symbolic reachable set computation of piecewise affine hybrid automata and its application to biological modelling: Delta-Notch protein signalling[END_REF] becomes

     ẋf = κ 1 f + κ 2 f s + (x g , θ 1 g ) s -(x f , θ 5 f ) -γ f x f ẋg = κ g 1 -s + (x g , θ 2 g )s -(x t , θ 1 t ) s -(x f , θ 4 f ) -γ g x g ẋt = κ t s + (x g , θ 2 g )s -(x t , θ 1 t )s + (x f , θ 4 f ) -γ t x t
As stated earlier, we know that x t ≤ θ 1 t after some finite time; this does not help us for further simplifications of this model. In the following, we will show that, after some finite-time, we have x t < θ 1 t , which will help us eliminate the x t equation. In order to do that, we first show that, after some finite time, x g ≤ θ 2

g . Indeed, if we suppose that x g > θ 2

g for all times, system (11) would become

     ẋf = κ 1 f + κ 2 f s -(x f , θ 5 f ) -γ f x f ẋg = κ g s + (x t , θ 1 t )s -(x f , θ 4 f ) -γ g x g ẋt = κ t s -(x t , θ 1 t )s + (x f , θ 4 f ) -γ t x t
which shows that x f reaches θ 5 f in finite time so that ẋg becomes equal to

ẋg = -γ g x g
This leads to the convergence of x g to 0 and thus to below θ 2 g , which is a contradiction. This shows that x g should reach θ 2 g in finite time when x g (0) > θ 2

g . An ensuing case-by-case analysis shows that the region where x g ≤ θ 2

g is invariant [START_REF] Madjebi | Etude de modèles de réseau de régulation génique[END_REF].

We will now show that x t is decreasing almost all of the time when x g ≤ θ 2 g and x t ≤ θ 1 t , that is in a region which we have shown to be reached in finite time and invariant. Detailing three cases, we have:

x g < θ 2 g or x f < θ 4 f : ẋt = -γ t x t . x g = θ 2
g and x f > θ 4 f : We have ẋg = -γ g x g < 0 at such a point and in a neighborhood surrounding each such point so that any solution directly enters the region where x g < θ 2

g (and consequently ẋt = -γ t x t , as we have seen).

x g = θ 2
g and x f = θ 4 f : We have ẋf = κ 1 f + κ 2 fγ f x f > 0 at this point and in a neighborhood surrounding it, so that any solution directly goes in one of the two previously described regions, where we have seen that x t is decreasing.

For any solution of (11), x t could only increase if x stayed in the second or third region, which we have shown not to be possible. We then have ẋt = -γ t x t for almost all times in the region of interest. After elimination of x t , we have to analyze the following system:

ẋf = κ 1 f + κ 2 f s + (x g , θ 1 g ) s -(x f , θ 5 f ) -γ f x f ẋg = κ g s -(x g , θ 2 g )s -(x f , θ 4 f ) -γ g x g (12)
At first sight, this analysis is not straightforward because this is a second order piecewise linear system with two thresholds in each direction, which theoretically gives rise to 9 regions. However, as is illustrated on Figure 10, some of the regions have the same dynamics and can be grouped together, giving rise to six regions. g and x f < θ 4 f , we have ẋg = κ g s -(x g , θ 2 g )γ g x g with κ g γ g > θ 2

g so that the line is attractive (black wall). Moreover,

ẋf = κ 1 f + κ 2 f -γ f x f > 0
so that x f is increasing and all solutions reach the end-point (x f , x g ) = (θ 4 f , θ 2 g ) in finite time. In some sense, each time the solution reaches this black wall, there is a reset taking place that sends the system to the end-point (θ 4 f , θ 2 g ) Along the line where x f = θ 5 f and x g > θ 1 g , we have

ẋf = κ 1 f + κ 2 f s -(x f , θ 5 f ) -γ f x f
so that this line also is a black wall (bearing in mind that

κ 1 f +κ 2 f γ f > θ 5 f ).
In addition, ẋg = -γ g x g so that x g is decreasing and all solutions reach the end-point (x f , x g ) = (θ 5 f , θ 1 g ) in finite time.

The observation of Figure 10 (as well as a detailed analysis of the linear systems in each of the regions) indicate that, eventually, the solutions oscillate around (x f , x g ) = (θ 4 f , θ 1 g ). Whether this oscillation is damped, neutrally stable or unstable is still unclear. It is clear, though, that the oscillation is bounded, as it cannot go beyond the black walls.

In order to analyze the oscillations, we will compute the first return map from and to the segment that links (θ 4 f , θ 1 g ) to (θ 4 f , θ 2 g ). We will therefore consider some (θ 4 f , x) as initial condition and compute the function f (x) such that (θ 4 f , f (x)) is the image of (θ 4 f , x) on the segment after one cycle around (θ 4 f , θ 1 g ). The computation of this first-return map can be handled in four steps, corresponding to the passages in the four regions surrounding (θ 4 f , θ 1 g ). The first step consists in computing the image of (θ 4 f , x), belonging to the initial segment, on the horizontal segment that links (θ 4 f , θ 1 g ) to (θ 5 f , θ 1 g ). The transition takes place in the region ∆ (4) so that [START_REF] Glass | The logical analysis of continuous non-linear biochemical control networks[END_REF] becomes

ẋf = κ 1 f + κ 2 f -γ f x f ẋg = -γ g x g ( 13
)
whose solution is

x f (t) = θ 4 f e -γ f t + κ 1 f +κ 2 f γ f (1 -e -γ f t ) x g (t) = x e -γgt (14)
In a first computation, we will suppose the absence of the vertical black wall, and use the dynamics (13) for both regions ∆ (4) and ∆ (5) : it is then straightforward to see that the solution impacts the target segment when x g (t) = θ 1 g , that is, at

t = t 1 (x) = ln(x)-ln(θ 1 g ) γg , so that x f (t 1 (x)) = θ 4 f e -γ f t1(x) + κ 1 f +κ 2 f γ f (1 -e -γ f t1(x) ) = θ 4 f θ 1 g x γ f γg + κ 1 f +κ 2 f γ f 1 - θ 1 g x γ f γg
However, we must account for the black-wall and it is possible that the actual solution hits this wall before reaching the target segment, so that the previously computed x f (t 1 (x)) > θ 5 f . In that case, the actual solution stays on the vertical black wall until it reaches the point (θ 5 f , θ 1 g ). Therefore the target of the point (θ 4 f , x) on the horizontal segment is

(f 1 (x), θ 1 g ) =   min   θ 4 f θ 1 g x γ f γg + κ 1 f + κ 2 f γ f   1 - θ 1 g x γ f γg   , θ 5 f   , θ 1 g   Similarly we can define (θ 4 f , f 2 (x)) as the image of (x, θ 1 g ) (with x ∈ [θ 4 f , θ 5 f ]
) on the vertical segment below the equilibrium, (f 3 (x), θ 1 g ) as the image of (θ 4 f , x) (with

x ∈ [0, θ 1 g ]) on the horizontal segment on the left of the equilibrium and (θ 4 f , f 4 (x)) as the image of (x, θ 1 g ) (with x ∈ [0, θ 4 f ]) on the initial segment. This yields

f 1 (x) = min θ 4 f θ 1 g x γ f γg + κ 1 f +κ 2 f γ f 1 - θ 1 g x γ f γg , θ 5 f f 2 (x) = θ 1 g θ 4 f - κ 1 f γ f x- κ 1 f γ f γg γ f f 3 (x) = θ 4 f θ 1 g - κg γg
x-

κg γg γ f γg + κ 1 f γ f   1 - θ 1 g - κg γg
x-

κg γg γ f γg   f 4 (x) = min   θ 1 g θ 4 f - κ 1 f +κ 2 f γ f x- κ 1 f +κ 2 f γ f γg γ f + κg γg   1 - θ 4 f - κ 1 f +κ 2 f γ f x- κ 1 f +κ 2 f γ f γg γ f    , θ 2 g    and f (x) = f 4 (f 3 (f 2 (f 1 (x)))) which has x = θ 1
g as a fixed point. It was then shown in [START_REF] Madjebi | Etude de modèles de réseau de régulation génique[END_REF] that f ′ (x) < 1 when x > θ 4 f , so that the sequence x n+1 = f (x n ), which represents the successive impacts on the initial segment converges to x = θ 4 f . We can then conclude that the cyclic solutions that surround (θ 4 f , θ 1 g ) are damped. This point is therefore a globally attractive equilibrium of (12) (cf. [START_REF] Farcot | Periodic solutions of piecewise affine gene network models: the case of a negative feedback loop[END_REF] that gives more general results in n dimensions for a negative feedback loop).

Having elucidated the dynamical behavior of the (x f , x g ) subsystem, we can now deduce the behavior of all other states. From the moment that we have x g < θ 2 g , it comes from ( 11) that ẋt = -γ t x t so that x t goes to 0. Once those three states are close to their equilibrium value, the remaining three equations become

   ẋc = κ 1 c -γ c x c ẋy = κ 1 y + κ 2 y -γ y x y ẋr = κ 1 r + κ 2 r -γ r x r
so that convergence of (x c , x y , x r ) towards

Comparison of the equilibria

It is interesting to compare both equilibria: we have x c x y x f x g x t x r

We see that most genes settle at different levels depending on the absence or presence of carbon sources. The most illustrative of the difference between the two states (carbon starved or not) is x r , which represents the concentration of stable RNAs and is a good indicator of the cellular growth. As expected, when carbon sources are depleted, the equilibrium level of x r is smaller than when carbon sources are abundant: when carbon shortage occurs, x r stays at a "house-keeping"-level whereas κ 1 r +κ 2 r γr , the equilibrium value in the presence of carbon sources, allows for fast cell growth. Also to be noted is the fact that x t = 0 in both cases; this does not mean that topA, the gene corresponding to x t , is useless. Indeed, when the carbon sources are either continuously present or absent, the effect of topA eventually dies down. However, in a time-varying environment, where nutrients are alternatively present and absent, an increase of the x t concentration can occur whenever x g > θ 2

g and x f > θ 4 f . TopA thus influences the transients.

Abstraction of the reduction method

We have seen that the preliminary model reduction has allowed for a simplification of the model analysis. Indeed, a global stability analysis of a 6-order model is no easy task, whereas there are various methods for the analysis of second order models. The reduction of the dimension of dynamical models is critical in the further development of the mathematical methods for genetic regulatory networks analysis because the networks typically are very large, so that it is rarely possible to study them directly. Classically, it has been attempted to apply time-scale separation methods, but these are mainly efficient for eliminating the fast metabolic components from mixed metabolic-genetic networks. Also balanced truncation methods have been introduced for genetic regulatory networks where inputs signals (action on the network) and output signals (measurements) are clearly identified ( [START_REF] Hardin | System reduction of nonlinear positive systems by linearization and truncation[END_REF][START_REF] Keil | Model reduction of modular systems using balancing methods[END_REF]). In this example, we have exploited the hierarchical triangular structure of the model arising after a finite time (this finite time allowed us to get rid of some of the interactions interfering with the triangular structure). We notice from graph theory that the identification of such a structure in the graph corresponding to the network is equivalent to the search for the strongly connected components of the graph. There are efficient algorithms to do so on large graphs, so that this model reduction method is tractable for the huge graphs that represent genetic regulatory networks (preliminary work on that subject has been done in [START_REF] Cristescu | Algorithmic study on genetic regulatory networks[END_REF] with links to GNA). Combining this approach with thresholds elimination allows for a progressive simplification of the graphs.

Challenges in PL models analysis

One of the major challenges in the analysis of models of genetic regulatory networks lies in the difficulty of obtaining accurate parameters. Therefore, one has to develop methods to identify the qualitative behavior of the system: when the parameters are linked together through inequalities (instead of being fixed at given values), we would like to be able to say something about the stability of the equilibria. Some interesting results have been obtained on that subject in [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF], as was shown in Section 2, and we would like to identify other cases where stability results can be deduced.

As we have seen in the analysis of E. coli, we are able to mathematically analyze PL models that are not trivial (dimension 6). However, actual genetic regulatory networks are much larger than that. It is therefore of paramount importance to develop methods that will help analyzing such large systems. Two major research directions are explored for that purpose: the model reduction approach (through balancing or through singular perturbations, in the linear case) and the separation of the original model into smaller, interconnected pieces that can be easily analyzed, as we have shown here.

Moreover, experimental techniques (e.g. gene deletion) are now available and allow to modify the production or degradation terms of some genes of the networks. This leads to problems of mathematical control of piecewise affine genetic networks, similar to more general problems for hybrid affine systems [START_REF] Habets | A control problem for affine dynamical systems on a full-dimensional polytope[END_REF]. The global problem is to control the trajectories through some prescribed sequence of rectangular regions. Some preliminary results have been obtained in [START_REF] Farcot | How to control a biological switch: a mathematical framework for the control of piecewise affine models of gene networks[END_REF]. For example, we have shown that a simple two-gene inhibitor system with a single equilibrium can be controlled to a bistable switch. We believe that interesting and original control problems are still to be solved in this domain.

Index

Fig. 1 .

 1 Fig. 1. Example of a genetic regulatory network of two genes (a and b), each coding for a regulatory protein (A and B).

Definition 1 .

 1 called a regulatory domain if none of the variables x i has a threshold value in D. In contrast, a domain D ∈ D is called a switching domain of order k ≤ n if exactly k variables have threshold values in D [25]. The corresponding variables x i are called switching variables in D. For convenience, we denote the sets of regulatory and switching domains by D r and D s respectively. It is also useful to define the concept of a supporting hyperplane for a domain. For every domain D ∈ D s of order k ≥ 1, define supp(D) to be the (nk)-dimensional hyperplane containing D. If D ∈ D r then we define supp(D) to be equal to Ω.

Definition 2 .Fig. 2 .

 22 Fig. 2. Illustration of the focal point φ(D 13) of a domain D13 in example[START_REF] Botsford | Cyclic AMP in prokaryotes[END_REF]

Fig. 3 .

 3 Fig. 3. Illustration of the definition of the focal set on a switching surface D according to the Filippov definition of solutions. The convex hull of the points φ(D 1 ) and φ(D 2) is simply the segment that links them, so that[START_REF] Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF] implies that φ(D) is the intersection of this segment with supp(D).

Fig. 4 .Assumption 1

 41 Fig. 4. Illustration of the non invariance of Φ(D): solutions with initial condition on Φ(D) stay inside the box Π(D) but do not necessarily stay in Φ(D)

25 D 4 D 9 D 14 1 D 2 D 3 D 4 D 5 Fig. 5 .

 254914123455 Fig. 5. Subdivision of the state-space in 25 domains and transition graph of system (1)

Theorem 2 .

 2 Assume Ω ⊂ IR n . Let D ∈ D s be a switching domain of order p ≥ 1 containing a singular equilibrium set Φ(D) that satisfies Assumption 1. If for all D ′ ∈ R(D), there is a transition from D ′ to D in the state transition graph, then Π(D) is asymptotically stable.

1 Fig. 7 .

 17 Fig.7. Network of key genes, proteins, and regulatory interactions involved in the carbon starvation network in E. coli. The notation follows, in a somewhat simplified form, the graphical conventions proposed by Kohn[START_REF] Kohn | Molecular interaction maps as information organizers and simulation guides[END_REF]. The contents of the boxes labeled 'Activation' and 'Supercoiling' are detailed in[START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF].

Fig. 8 .

 8 Fig. 8. Entry into stationary phase: qualitative temporal evolution of the proteins and stable RNA concentration in a depleted environment with the organisms being at the equilibrium of the exponential phase at the initial time. Convergence to one domain is detected (the domain where xc > θ 3 c , xy = θ 3 y , x f < θ 1 f , xg = θ 2 g , xt < θ 1 t and xr < θr)

Fig. 9 .

 9 Fig.9. Entry into exponential phase: qualitative temporal evolution of the proteins and stable RNA concentration in a rich environment with the organisms being at the equilibrium of the stationary phase at the initial time. Oscillations of the x f and xg states is detected.

Fig. 10 .

 10 Fig. 10. Illustration of the vector field and the various regions for system (12). The thick black lines indicate where sliding modes can occur.

  QS 43 QS 71 QS 69 QS 67 QS 65 QS 35 QS 33 QS 34 QS 5 QS 7 QS 43 QS 71 QS 69 QS 67 QS 65 QS 35 QS 33 QS 34 QS 5 QS 7 QS 43 QS 71 QS 69 QS 67 QS 65 QS 35 QS 33 QS 34 QS 5 QS 7 QS 43 QS 71 QS 69 QS 67 QS 65 QS 35 QS 33 QS 34 QS 5 QS 7 QS 43 QS 71 QS 69 QS 67 QS 65 QS 35 QS 33 QS 34 QS 5 QS 7 QS 43 QS 71 QS 69 QS 67 QS 65 QS 35 QS 33 QS 34 QS 5 QS 7 QS 43 QS 71 QS 69 QS 67 QS 65 QS 35 QS 33 QS 34 QS 5 QS 7

  QS 427 QS 255 QS 107 QS 109 QS 69 QS 71 QS 39 QS 41 QS 53 QS 95 QS 427 QS 255 QS 107 QS 109 QS 69 QS 71 QS 39 QS 41 QS 53 QS 95 QS 427 QS 255 QS 107 QS 109 QS 69 QS 71 QS 39 QS 41 QS 53 QS 95 QS 427 QS 255 QS 107 QS 109 QS 69 QS 71 QS 39 QS 41 QS 53 QS 95

	QS 2 QS 3	QS 427	QS 7	QS 255 QS 107 QS 109 QS 69 QS 71	QS 49	QS 51	QS 39 QS 41	QS 43	QS 53 QS 95	QS 105
	QS 2 QS 3		QS 7		QS 49	QS 51		QS 43		QS 105
	QS 2 QS 3		QS 7		QS 49	QS 51		QS 43		QS 105
	QS 2 QS 3		QS 7		QS 49	QS 51		QS 43		QS 105
	QS 2 QS 3		QS 7		QS 49	QS 51		QS 43		QS 105
	QS 2 QS 3	QS 427	QS 7	QS 255 QS 107 QS 109 QS 69 QS 71	QS 49	QS 51	QS 39 QS 41	QS 43	QS 53 QS 95	QS 105
	QS 2 QS 3	QS 427	QS 7	QS 255 QS 107 QS 109 QS 69 QS 71	QS 49	QS 51	QS 39 QS 41	QS 43	QS 53 QS 95	QS 105
					1 t and xr < θr)

Acknowledgments

This research is supported by the HYGEIA project of the NEST-Adventure program of the European Union (NEST-004995) on Hybrid Systems for Biochemical Network Modeling and Analysis.