domain: info.info-au math spi sdv

text:
Alfonso Carta
email: alfonso.carta@gmail.fr

Madalena Chaves
email: madalena.chaves@inria.fr

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

Continuous-switch piecewise quadratic models of biological networks: application to bacterial growth

An extension of the class of piecewise linear (PL) systems is proposed to model gene expression dynamics dependent on dilution due to cell growth rate. The growth rate is modeled as the weighted minimum of two or more limiting gene products responsible for bacterial growth. The production terms are still piecewise constant, but now the degradation terms are piecewise quadratic (PQ). This new mathematical formalism exhibits continuous switches between PQ modes. We first study the novel dynamical behavior generated by the nonlinear terms at the regions of discontinuity of the vector fields, showing that the sliding motion configurations occurring in PL systems can further lead to damped convergent oscillations or periodic behavior in PQ systems. As an application, a core model of the bacterial gene expression machinery is studied with the goal of externally tuning the growth rate of cells. This system may exhibit several behaviors including bi-mode bistability or damped oscillatory behavior.

Introduction

The concepts and methods from systems and control theory have a large potential to significantly help in systems and synthetic biology research, where the goal is to develop and apply engineering tools to control cellular behavior and achieve desired functions [START_REF] Mukherji | Synthetic biology: understanding biological design from synthetic circuits[END_REF]. Synthetic biology aims at constructing novel biological circuits in the cell and most recent designs focus on assembling components from the cell transcription machinery, which includes the genes to be expressed, their promoters, RNA polymerase, ribosomes and transcription factors, all serving as potential individual engineering components [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF][START_REF] Gardner | Construction of a genetic toggle switch in Escherichia coli[END_REF][START_REF] Wang | Engineering modular and orthogonal genetic logic gates for robust digital-like synthetic biology[END_REF].

In this context, control-based approaches are increasingly being used in synthetic biology [START_REF] Ang | Considerations for using integral feedback control to construct a perfectly adapting synthetic gene network[END_REF][START_REF] Leduc | How do control-based approaches enter into biology?[END_REF][START_REF] Menolascina | Analysis, design and implementation of a novel scheme for in-vivo control of synthetic gene regulatory networks[END_REF][START_REF] Yang | External control of the GAL network in S. cerevisiae: A view from control theory[END_REF] where some control theoretical results are applicable, although with various limitations due to biological constraints [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: The bistable switch example[END_REF][START_REF] Sontag | Some new directions in control theory inspired by systems biology[END_REF].

Here, our goal is to introduce a novel mathematical formalism to qualitatively model gene expression and dilution due to cell growth. In fact, one of the aims of systems and synthetic biology is to link molecular-level mechanisms (e.g. gene expression) to cell-level behavior (e.g. growth rate) [START_REF] Mukherji | Synthetic biology: understanding biological design from synthetic circuits[END_REF]. To this aim, we introduce a model where bacterial growth rate is limited by different factors, which ultimately leads to a continuous-switch piecewise quadratic (CSPQ) formalism-derived from piecewise linear (PL) systems [START_REF] Jong | Modeling and simulation of genetic regulatory systems : A literature review[END_REF][START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF].

An application of this modeling formalism is the control the cells' growth rate by designing appropriate laws. Growth control is essential in industrial biotechnology and fundamental research of this kind could pave the way to novel types of antimicrobial strategies [START_REF] Khalil | Synthetic biology: applications come of age[END_REF]. We analyze a core model of the gene expression machinery of the bacterium Escherichia coli, where the growth rate is controlled externally by tuning the synthesis of a component of the gene expression machinery (RNA Polymerase). This type of control can be easily implemented, for instance, by means of inducers that activate synthetic inducible promoters [START_REF] Kaern | The engineering of gene regulatory networks[END_REF].

A first formulation and preliminary results on CSPQ systems were introduced in [START_REF] Carta | A class of switched piecewise quadratic systems for coupling gene expression with growth rate in bacteria[END_REF]. Here, new properties of these systems are investigated. Section 2 recalls and extends the formulation of CSPQ systems. In Section 3 we state and prove some theorems on the stability of each piecewise quadratic (PQ) subsystem; in particular, new results are given on the properties of sliding modes and their differences with respect to the "standard" PL system. Section 4 studies local stability of the CSPQ system. Section 5 discusses an application to an open-loop control system of bacterial gene expression. Section 6 summarizes our results and conclusions.

Figure 1: Example of gene expression data points [START_REF] Brewster | The transcription factor titration effect dictates level of gene expression[END_REF] and possible activation functions. Step (dash-dotted), logoid (dashed), or a C 1 Hill-type function (solid).

Continuous-switch piecewise quadratic systems

As introduced in [START_REF] Carta | A class of switched piecewise quadratic systems for coupling gene expression with growth rate in bacteria[END_REF], the switched piecewise quadratic systems are an extension of piecewise linear (PL) systems, a class of well known qualitative models originally introduced by [START_REF] Glass | The logical analysis of continuous, non-linear biochemical control networks[END_REF] and further studied in [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF][START_REF] Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF][START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF]. The PL systems are very useful to study genetic networks in the absence of detailed data or parameter knowledge (see, for instance, [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]) due to their suitableness for theoretical analysis and simulations using algorithmic tools and software. However, PL systems consider only linear protein degradation without taking into account cell growth, required to represent the connections between molecular and cellular mechanisms.

To overcome this problem, we consider a general model structure to describe the dynamics of the concentrations of a system of n species, x = (x 1 , ...,

x n ) T ∈ R n ≥0 : ẋi = f i (x) -(µ(x) + γ i ) x i , 1 ≤ i ≤ n . (1)
where f i (x) is the synthesis rate and γ i represents the typical linear degradation constant for species i, and µ(x) represents dilution in the concentrations, proportional to the bacterial growth rate and equal for all species. These terms are next defined. The function f i : R n ≥0 → R ≥0 represents the expression rate of the gene i depending on state x. Several experiments [START_REF] Brewster | The transcription factor titration effect dictates level of gene expression[END_REF][START_REF] Wang | Engineering modular and orthogonal genetic logic gates for robust digital-like synthetic biology[END_REF] show that gene expression often follows an increasing (or decreasing) curve h(x) as a function of an activator (or inhibitor) x, between two essentially flat regions:

h(x) =    0, x < x low h(x), x low < x < x high 1, x > x high (2)
where h(x) may take any form. The logoid (linear h) is studied in [START_REF] Mestl | A mathematical framework for describing and analysing gene regulatory networks[END_REF]. A C 1 alternative is a Hill-type function b+ax 2 /(x 2 +θ 2 ) (see Fig. 1). To facilitate theoretical analysis, in the class of piecewise linear systems, activation functions are approximated by step functions, so that the function f i has the general form:

f i (x) = n ∑ l=1 k il b il (x)
where k il ≥ 0 is a rate parameter and b il (x) is a sum of products of step functions, s + , s -. These are defined in terms of a threshold parameter θ > 0:

s + (x i , θ ) = 1 if x i > θ 0 if x i < θ ; s -(x i , θ ) = 1 -s + (x i , θ ).
Note that, in the case µ(x) ≡ 0, system (1) recovers the classical piecewise linear (PL) systems [START_REF] Jong | Modeling and simulation of genetic regulatory systems : A literature review[END_REF]. Assuming that bacterial growth rate is limited by the amount of one or more generic cellular components [START_REF] Carta | A class of switched piecewise quadratic systems for coupling gene expression with growth rate in bacteria[END_REF] which are necessary to sustain the gene expression machinery of the cell, let the growth rate µ : R n ≥0 → R ≥0 be given as:

µ(x) = min 1≤i≤n { µ i x i } (3)
where µ i ≥ 0 are proportion factors depending, for instance, on the carbon source used. If µ k = 0, then growth doesn't depend on species k. Similar expressions have been used recently for ribosomal regulation in E. coli [START_REF] Shachrai | Cost of unneeded proteins in E. coli is reduced after several generations in exponential growth[END_REF]. This formulation appears also in other contexts, such as in ecology where the specific growth rate of species is often determined by the resource that is most limiting according to Liebig's "law of the minimum" [START_REF] Huisman | Biodiversity of plankton by species oscillations and chaos[END_REF]. System (1) belongs to the class of switched systems [START_REF] Liberzon | Basic problems in stability and design of switched systems[END_REF] in which the growth rate µ acts as the rule that orchestrates the switching between the different subsystems (4):

mode-r: ẋi = f i (x) -[µ r x r + γ i ] x i , 1 ≤ i ≤ n. ( 4
)
Each mode is piecewise quadratic (PQ), hence to study the dynamics of the (full) CSPQ system (1) we will first characterize the dynamics of its PQ modes (4), and then investigate the properties arising from the switching condition.

3 The PQ subsystem: dynamical study A dynamical study is provided for mode-r in [START_REF] Brewster | The transcription factor titration effect dictates level of gene expression[END_REF], that is when µ(x) = µ r x r . The PQ system can be written as:

ẋ = f (x) -d(x r )x, (5)
where f = ( f 1 , ..., f n ) and d(x r ) = diag(µ r x r +γ 1 , ..., µ r x r +γ n ), where diag is the diagonal matrix corresponding to the vector. First, note that the non-negative orthant remains invariant, since x i = 0 implies ẋi ≥ 0. Second, note that all solutions remain bounded, since

f i (x) -d i (x r )x i < max x∈R n ≥0 f i (x) -γ i x i
, where max x∈R n ≥0 f i (x) < ∞, by definition of f . Without loss of generality, the dynamics of the PQ subsystem can be studied in the n-dimensional state-space Ω = Ω 1 × Ω 2 × ... × Ω n , where each Ω i is defined by

Ω i = {x i ∈ R ≥0 |0 ≤ x i ≤ max i } with max i = max x∈R n ≥0 f i (x)/γ i .
The set Ω is thus invariant for system [START_REF] Carta | A class of switched piecewise quadratic systems for coupling gene expression with growth rate in bacteria[END_REF]. Each protein 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 , ..., θ p i i and let θ 0 i = 0, θ p i+1 i = max i . These thresholds partition Ω into hyper-rectangular regions called 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 |θ j i i < x i < θ j i +1 i , for j i ∈ {0, ..., p i } D i = x i ∈ Ω i |x i = θ j i i , for j i ∈ {1, ..., p i }.
Let D be the set of domains in Ω. A domain D ∈ D is called a regular domain if none of the variables x i has a threshold value in D (it is the full hyper-rectangle). In contrast, a domain D ∈ D is called a threshold domain of order k ≤ n if exactly k variables have threshold values in D (in [START_REF] Mestl | A mathematical framework for describing and analysing gene regulatory networks[END_REF] threshold domains are called switching domains, but we avoid this definition to prevent misunderstandings with switched system). The corresponding variables x i are called threshold variables in D. The two sets of domains are respectively denoted by D rg andD th .

Solutions and Stability in Regular Domains

For any regular domain D, for all x ∈ D the function f (x) ≡ f D is a constant, so (5) can be written as

ẋ = f D -d(x r )x, (6)
as done in [START_REF] Carta | A class of switched piecewise quadratic systems for coupling gene expression with growth rate in bacteria[END_REF]. Note that ( 6) is a hierarchical system, since the function f r (x) depends only on x r while f i (x) depends only on x i and x r , but not on x j , j = r, i. Thus, for any x(t 0 ) ∈ D the unique solution of (6) can be found explicitly by solving first the r-component of ( 6), which is an autonomous differential equation, and then successively solving for the i-components (i = r), having substituted x r (t) into them. The focal points of the PQ system (5) have the form:

xi = η( xr , f D i , µ r , γ i ) = f D i µ r xr + γ i , i = 1, ..., n, i = r xr = ϕ( f D r , µ r , γ r ) = -γ r + γ 2 r + 4µ r f D r 2µ r , (7)
A focal point is a regular equilibrium point provided that it belongs to its respective regular domain, i.e. Φ(D) ∈ D. We now expand some basic results from [START_REF] Carta | A class of switched piecewise quadratic systems for coupling gene expression with growth rate in bacteria[END_REF] on the stability of these points.

Theorem 1 Consider a regular domain D = D 1 × • • • × D n ∈ D rg with D i = (θ j i i , θ j i +1 i
). Assume that Φ(D) = ( x1 , ..., xn ) T ∈ D. Then Φ(D) is a locally asymptotically stable equilibrium of the system (6) with respect to its domain D. If, in addition, the following conditions are satisfied: (a)

θ j i +1 i > f D i γ i and (b) θ j i i < f D i µ r θ jr +1 r +γ i
, for all i = r, then Φ(D) is globally asymptotically stable with respect to D.

Proof:

If Φ(D) ∈ D, then θ j i i < xi < θ j i +1
i for all i. Note that variable x r is autonomous and its equation can be written ẋr = -µ r (x rxr )(x rω), where x is given by [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: The bistable switch example[END_REF] and ω < 0 is the second root of the quadratic equation. Since θ j r r < xr < θ j r +1 r , it follows immediately that ẋr > 0 at x r = θ j r r and ẋr < 0 at

x r = θ j r +1 r . Hence x r (t) → xr . For variable x i , i = r, one can write ẋi = (µ r x r + γ i ) f D i µ r x r + γ i -x i ,
where µ r x r + γ i > 0. The local stability follows by linearization around Φ(D), which yields a triangular system with negative diagonal entries. To prove global stability with respect to D, assume that (a) and (b) hold. Then

θ j i i < f D i µ r θ j r +1 r + γ i < f D i µ r x r + γ i < f D i γ i < θ j i +1 i .
Thus, x i ≥ θ j i +1 i implies ẋi < 0 and x i ≤ θ j i i implies ẋi > 0 and both imply that x i cannot leave the interval D i . Therefore, the regular domain D is a forward-invariant set for the system. Since the system is in a cascade, x r (t) → xr also implies x i (t) → xi , which concludes the first part.

If (a) or (b) fail, one can show that stability is only local. Suppose (a) is false so θ j i +1 i < f D i /γ i for some i. For any ε > 0 small enough, consider an initial condition x 0 satisfying:

x 0i = θ j i +1 i -ε < f D i γ i , f D i µ r x r0 + γ i > θ j i +1 i
Then ẋi > 0 for some time interval and, letting ε > 0 tend to zero, some trajectories will leave the domain D through the boundary x i = θ j i +1 i . Note that, for PL systems, global asymptotic stability with respect to D is guaranteed: since µ r = 0, the focal points are given by x i = f D i /γ i for all i and conditions (a)-(b) are automatically satisfied. In general, all trajectories in a regular domain D flow towards the focal point Φ(D) asymptotically unless they reach a boundary and leave the domain D. Since the step functions are not defined for x i = θ j i i , we next analyze in detail the dynamics at threshold domains.

Solutions and Stability in Threshold Domains

Since the function f is a sum of products of step functions, it is discontinuous at the threshold points. In order to provide the existence and the possibility for solutions to be continued on all domains, we have to define the righthand side of system (5) at these points. To do this, we use a construction originally proposed by Filippov [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF] and then applied to PL systems [START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF]. The method consists of extending the system (5) to a differential inclusion,

ẋ ∈ H(x), (8)
where H is a set valued function (i.e. H(x) ⊆ R n ). If D ∈ D rg , then we define H simply as

H(x) = { f D -d(x r )x}, ∀x ∈ D. (9)
If D ∈ D th , we define H as

H(x) = co({ f D -d(x r )x | D ∈ R(D)}), ∀x ∈ D, (10)
where

R(D) = {D ∈ D rg | D ⊆ ∂ D }
is the set of all regular domains with D in their boundary, and co(X) is the closed convex hull of X. For threshold domains, H(x) is typically multi-valued so solutions of the differential inclusion are defined as follows.

Definition 2 A solution of (8) 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 ].
To see how to construct H(x) at a discontinuity point, consider the case where x belongs to a threshold domain S separating two regular domains D A and D B . Hence,

H(x) = co({ f D A -d(x r )x, f D B -d(x r )x})
represents the segment joining the endpoints of the vectors g A = f D Ad(x r )x and g B = f D Bd(x r )x. Roughly speaking, trajectories can cross S if the vector fields g A and g B point in the same sense across S, slide along S if g A and g B point in opposite senses towards S and be repelled from S if g A and g B point in opposite senses away from S (see some examples of vector fields in Fig. 2). The last two cases are known as stable and unstable sliding motion in the literature [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF]. Moreover, the velocity of the sliding motion (stable or unstable) on S is given by

ẋ = α f D A + (1 -α) f D B -d(x r )x. (11)
The dynamics of PQ systems is similar to that of PL systems at threshold domains in x r , but they exhibit new behavior on the emergence of sliding motion at thresholds of the form x i , i = r (see the discussion in Sections 3.

Φ(D) = supp(D) ∩ {x : 0 ∈ H(x)}, (12)
where H(x) is defined as in [START_REF] Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF].

Hence, Φ(D) for D ∈ D th is the set containing all the equilibrium points of the differential inclusion [START_REF] Jong | Modeling and simulation of genetic regulatory systems : A literature review[END_REF], which lie on supp(D). Thus, Φ(D) can be a singleton, but more generally is a closed convex bounded set and hence is referred to as a focal set.

From now on, we make a technical assumption saying that the focal set Φ(D) does not intersect the supporting hyperplanes of threshold domains at the boundary of D:

Assumption 1 For all D ∈ D, we assume that Φ(D) ∩ supp(D ) = / 0, for all D ⊂ ∂ (D).
System (5) may also have singular equilibrium (to be studied on a case-by-case basis), defined as a point x belonging to some threshold plane and satisfying the condition 0 ∈ H( x).

Dynamics at threshold domains x r = θ r

Consider the case when two regular domains, D A and D B , share a boundary segment S with x r = θ r and their focal points are on opposite sides of this threshold domain:

sup{x r : x ∈ D A } < θ r < inf{x r : x ∈ D B } and Φ r (D A ) > θ r > Φ r (D B ).
Since the ẋr equation depends only on x r , a solution of the form [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF] can be constructed where the value of α is constant for all x ∈ S and thus a sliding motion occurs along the threshold domain, as in the case of PL systems. Moreover, the vector field on S converges towards the focal set as defined in [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF] (see [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF][START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF] for more details).

= γ i = 1, f D A r = 1.2, f D A i = 2.3, f D B r = 3.2, f D B i = 1.2, (b) the case f D A i > f D B i = 3, (c) the case f D A i = f D B i , and (d) the case f D A i < f D B i = 4.5.

Dynamics at threshold domains

x i = θ i , i = r
To give a clear idea of the different dynamical behaviors that may arise, here we will study only the 2-dimensional case, with variables (x i , x r ). Let the boundary segment S between D A and D B be of the form x i = θ i , i = r:

sup{x i : x ∈ D A } < θ i < inf{x i : x ∈ D B } (13)
and

Φ i (D A ) > θ i > Φ i (D B ), (14)
which characterizes the location of Φ i (D), the i-th coordinate of region D's focal point. Three different cases can happen, shown in Fig. 2(b)-(d), depending on the configuration of the nullclines ẋi = 0 passing through the focal points, which are of the form ν i (x r ) = f D i /(µ r x r + γ i ) given in [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: The bistable switch example[END_REF]. The three configurations are characterized by the sign of

f D A i -f D B i
which implies an ordering between the two nullclines (dashed and dot-dashed curves in Fig. 2) and introduces a subdivision of the threshold domain x i = θ i into regions where the vector fields have the same or opposite signs. Theorem 2 summarizes the main results. For comparison, the "standard" PL case is shown in Fig. 2(a).

Theorem 2 Consider a 2-dimensional system of the form (5), with coordinates (x i , x r ), and let D A and D B be two adjacent regular domains of system sharing a boundary segment S of the form x i = θ i , i = r. Assume conditions (13) and ( 14) are satisfied and that the focal set Φ(S), given by [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF], belongs to S. The dynamics of the system at the threshold domain S can then be characterized as follows:

1. If f D A i > f D B
i , then the system has a stable sliding motion on S towards the focal set Φ(S);

2. If f D A i = f D B
i , then the system has a damped oscillation on D A and D B towards the focal set Φ(S);

3. If f D A i < f D B
i , then the focal set Φ(S) is unstable, and there is an unstable sliding motion on the threshold domain.

Proof: The 2-dimensional system has the form:

ẋr = f D r -(µ r x r + γ r )x r = -µ r (x r -Φ r (D))(x r + ω D ), ẋi = f D i -(µ r x r + γ i )x i , (15)
where Φ r (D) > 0, ω D > 0 are given by

Φ r (D) = -γ r + γ 2 r + 4µ r f D r 2µ r , ω D = γ r + γ 2 r + 4µ r f D r 2µ r .
In general, let ν A (x r ) and ν B (x r ) represent the nullclines ẋi = 0 of the system passing through the focal points Φ(D A ) and Φ(D B ). Let R A and R B denote the x r coordinate of the intersection between the nullclines and the segment S:

θ i = ν A (R A ), θ i = ν B (R B ). with R D α = f Dα i µ r θ i -γ i µ r , α ∈ {A, B}. 1.
In this case, from [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: The bistable switch example[END_REF] we see that ν A (x r ) > ν B (x r ) for all x r ∈ S and also R A > R B . In particular, the focal set is between ν A (x r ) and ν B (x r ). Therefore, from the vector fields on D A and D B it becomes clear that

x r < R B : f D A i -d(x r )θ i > 0, f D B i -d(x r )θ i > 0, R B < x r < R A : f D A i -d(x r )θ i > 0 > f D B i -d(x r )θ i , x r > R A : f D A i -d(x r )θ i < 0, f D B i -d(x r )θ i < 0.
This implies that the vector fields can be normally continued "outside the nullclines", for all x r < R B (transition from D A to D B ) and x r > R A (transition from D B to D A ). Between the nullclines, a stable sliding motion can be constructed according to [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF] (see Fig. 2(b)).

2. In this case, it follows that the nullclines coincide: ν A (x r ) ≡ ν B (x r ) which is also the convex hull of the two vector fields. Therefore,

f D A i = f D B i = f i , R A = R B = R and we have: [x r -R] • [ f i -d(x r )θ i ]
< 0, so the vector fields can be continued normally from D A to D B for x r < R, and in the opposite sense for x r > R (see Fig. 2(c)). To prove that the oscillatory dynamics is converging towards Φ(S) ∩ H(x) = (θ i , R), we will translate system [START_REF] Glass | The logical analysis of continuous, non-linear biochemical control networks[END_REF] towards this point and then use Lyapunov's direct method. Thus, letting z i = x iθ i and z r = x r -R, after simplification we obtain for α ∈ {A, B}:

żr = -g α (z r ), żi = -θ i µ r z r -(µ r z r + f i θ i )z i , (16)
with

g α (z r ) := µ r (z r + R -Φ r (D α ))(z r + R + ω D α ) (17)
for z i ≥ -θ i and z r ≥ -R. Note that x ∈ D A is equivalent to z i < 0 and similarly x ∈ D B is equivalent to z i > 0. Consider the following Lyapunov function candidate for system [START_REF] Gourse | rRNA transcription and growth rate-dependent regulation of ribosome sysnthesis in Escherichia coli[END_REF]:

V (z r , z i ) = z 2 r -b A z i , z i ≤ 0 z 2 r + b B z i , z i > 0, (18)
where

b α = 2 µ r θ i |g α (0)|, α ∈ {A, B}. (19)
It is clear that V is continuous and positive definite on R 2 since V (0, 0) = 0 and V (z i , z r ) > 0 for all (z i , z r ) = (0, 0). Also, V is proper on R 2 , since the set {z ∈ R 2 : V (z) ≤ L} is compact for each L > 0. Let z * be any initial condition in D A ∩ D B and z(t) the corresponding trajectory of ( 16). The derivative of V along z(t), V (t), is given by:

V (t) = -2z r g A (z r ) + b A θ i µ r z r + b A (µ r z r + f i θ i )z i , z i ≤ 0 -2z r g B (z r ) -b B θ i µ r z r -b B (µ r z r + f i θ i )z i , z i > 0.
We need to analyze the sign of the terms in dV /dt. For the last term, note that -R ≤ z r < 0 implies :

µ r |z r | < µ r R < µ r R + γ i = f i θ i ⇒ µ r z r + f i θ i > 0.
Therefore, b α (µ r z r + f i θ i ) > 0 for all z r ≥ -R, which implies

     b A (µ r z r + f i θ i )z i ≤ 0, z i < 0 0, z i = 0 -b B (µ r z r + f i θ i )z i < 0, z i > 0. (20)
Consider now the first two terms for z i < 0:

-2z r g A (z r ) + b A θ i µ r z r = -2z r (g A (z r ) -g A (0)) < 0, (21)
for all z r ≥ -R and z r = 0. We used the fact that |g A (0)| = g A (0), due to R > Φ r (D A ), and also that g A is an increasing function for all z r ≥ -R (its derivative is given by dg A /dz r = 2(z r + R) + γ r /µ r ). Similarly, for z i > 0, we have |g B (0

)| = -g B (0), due to R < Φ r (B), hence -2z r g B (z r ) -b B θ i µ r z r = -2z r (g B (z r ) -g B (0)) < 0, (22)
for all z r ≥ -R and z r = 0, since g B is also an increasing function for all z r ≥ -R (its derivative is the same as for g A ).

Putting together the inequalities ( 20)- [START_REF] Liberzon | Basic problems in stability and design of switched systems[END_REF] we conclude that dV /dt < 0 for all (z i , z r ) = (0, 0). For points on the threshold z i = 0, although dV /dt is discontinuous, ( 21) and [START_REF] Liberzon | Basic problems in stability and design of switched systems[END_REF] show that the values at z i = 0 + and z i = 0 -are both negative, so V (z(t)) is a strictly decreasing function of time, for any trajectory z(t). Then V is a Lyapunov function in the sense of Definition 2 of [START_REF] Bacciotti | Stability and stabilization of discontinuous systems and nonsmooth Lyapunov functions[END_REF], and Theorem 3 in this same reference can be used to conclude that the point (θ i , R) is indeed an equilibrium point of system [START_REF] Gourse | rRNA transcription and growth rate-dependent regulation of ribosome sysnthesis in Escherichia coli[END_REF]. 8 3. In this case, the nullclines satisfy the opposite ordering to those in the first case, that is ν A (x r ) < ν B (x r ) for all x r ∈ S and R A < R B , which leads to:

x r < R A : f D A i -d(x r )θ i > 0, f D B i -d(x r )θ i > 0, R A < x r < R B : f D A i -d(x r )θ i < 0 < f D B i -d(x r )θ i , x r > R B : f D A i -d(x r )θ i < 0, f D B i -d(x r
)θ i < 0. This implies that the vector fields can again be normally continued "outside the nullclines", for all x r < R A (transition from D A to D B ) and x r > R B (transition from D B to D A ). However, for R A < x r < R B , the vector fields generate an unstable sliding motion on the threshold x i = θ i . Depending on the parameters, a periodic orbit may appear between the two domains (see Fig. 2(d)).

Local stability analysis of the CSPQ system

System (1) belongs to a very small family within the wide class of switching systems (see, for instance, [START_REF] Branicky | Stability of switched and hybrid systems[END_REF]) since its vector field is continuous at all points of the switching surfaces, except at the intersections with threshold domains (where the vector field is not uniquely defined). To study existence and local stability of equilibria of CSPQ systems, for each µ r > 0, we will define an open region χ r ⊂ R n , in which system ( 1) is in mode-r. Let M = {1 ≤ m ≤ n : µ m > 0} be the set of indexes of variables contributing to the growth rate and then define for all r ∈ M:

χ r = {[x 1 , ..., x n ] T ∈ R n : µ r x r < µ i x i , ∀i = r}. ( 23
)
Define also the switching surface between modes r and p:

S r,p = {[x 1 , ..., x n ] T ∈ R n : µ p x p -µ r x r = 0} . (24)
Let Ψ m and Λ m (m = 1, . . . , n) be the sets containing, respectively, all the locally stable and all the unstable points of the m-mode (as determined in Section 3). For simplicity, it is assumed that equilibria of each mode do not lie on switching surfaces S r,p .

Assumption 2 For all r, p ∈ M: Ψ m ∩ S r,p = / 0, Λ m ∩ S r,p = / 0 for m ∈ {r, p}.

Therefore, a stable (resp. unstable) point of the m-mode is also a stable (resp. unstable) point of the CSPQ (1), only if it belongs to χ m , the region in which the m-mode is active. Let

Ψ = ∪ m∈M (Ψ m ∩ χ m ) . ( 25
)
denote the set of locally stable points of (1). Prop. 1 shows that two modes cannot mutually contain each other's focal points, thus eliminating possible oscillatory behavior around the switching surfaces.

Proposition 1 Given any D ∈ D rg , let Φ P and Φ R represent its focal points in modes p and r. Then either

Φ P / ∈ χ R or Φ R / ∈ χ P .
Proof: Given regular domain D the focal points Φ P , Φ R are calculated according to (7):

Φ P p = -γ p + γ 2 p + 4µ p f D p 2µ p , Φ P r = f D r µ p Φ P p + γ r ( 26
) Φ R p = f D p µ r Φ R r + γ p , Φ R r = -γ r + γ 2 r + 4µ r f D r 2µ r (27)
To prove the statement by contradiction, assume that it is false, so both Φ P ∈ χ R and Φ R ∈ χ P . Then we can write µ p Φ P p > µ r Φ P r and µ p Φ R p < µ r Φ R r and, after some algebra, these inequalities lead to both -γ r + γ 2 r + 4µ r f D r > -γ p + γ 2 p + 4µ p f D p and its opposite. Hence the Proposition holds. Finally, consider the dynamics at the intersection of the switching surface with a regular domain S D = D∩S r,p . The continuity of the vector field at S D together with Assumption 2 and Prop. 1 imply that trajectories will eventually evolve away from S D . In particular, it rules out sliding motion on a switching surface, for the generic case. The global dynamics of the CSPQ system can clearly generate many different configurations, given the possible combinations of PQ systems. A general overview of different possibilities can be found, for instance, in [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF] (chapter 9).

Growth rate control through RNAP and ribosomes

As an application of the CSPQ formalism, we focus on the gene expression machinery of the bacterium Escherichia coli during the exponential phase [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF], with the aim of adjusting the growth rate of the cells to a desired level. We will further analyze the 2-dimensional CSPQ open loop model developed in previous work [START_REF] Carta | A class of switched piecewise quadratic systems for coupling gene expression with growth rate in bacteria[END_REF]. It is a model of a mutant E. coli inspired by the experiments in [START_REF] Tan | Emergent bistability by a growth-modulating positive feedback circuit[END_REF], where an engineered inducible-promoter is used to externally control the expression of RNAP (see parameter u in ( 28)).

Our model is based on the idea that growth rate is related to the amount of proteins in the cell. All proteins are produced in a two-step process [START_REF] George | A unified theory of gene expression[END_REF]: (i) transcription is catalyzed by RNA polymerase (RNAP), an enzyme which allows the synthesis of mRNA from DNA; (ii) during translation, the mRNA is translated into proteins by ribosomes, complex molecules composed of various proteins and rRNAs.

The aim of our study is to investigate the possibility of generating single or multiple equilibria, and under which conditions. For example, switching between different equilibria of the bacterial system might be a strategy for reducing protein production costs associated with a target protein [START_REF] Shachrai | Cost of unneeded proteins in E. coli is reduced after several generations in exponential growth[END_REF]. Thus, in Section 5.2, we focus on characterizing the most interesting stability scenarios.

CSPQ model of the RNAP-ribosomes system

The two variables of our model, x p , x r ∈ R ≥0 , are the concentrations of RNAP and ribosomes both of which play a pivotal role in protein synthesis and contribute to cell growth (Fig. 3). Using (3) for the growth rate, our model becomes:

ẋr = k 1 r s + (x p , θ 1 p ) + k 2 r s + (x p , θ 2 p ) -(µ(x) + γ r )x r ẋp = u k 0 p s + (x p , θ 1 p )s + (x r , θ 1 r ) + k 1 p s + (x p , θ 2 p )s + (x r , θ 2 r ) -(µ(x) + γ p )x p (28)
where µ(x) = min(µ p x p , µ r x r ) and u ∈ [0, 1] is a constant control parameter that modulates the expression of RNAP and eventually tunes growth rate. In synthetic biology constructions, u can be a chemical inducer such as aTc or IPTG [START_REF] Gardner | Construction of a genetic toggle switch in Escherichia coli[END_REF].

The x r equation describes synthesis of ribosomes limited by the production rate of stable-RNAs, which in turn are regulated at the level of transcription by RNAP [START_REF] Gourse | rRNA transcription and growth rate-dependent regulation of ribosome sysnthesis in Escherichia coli[END_REF]. Stable-RNAs are produced from two promoters: the house-keeping P2 (activated by a low concentration of RNAP, θ 1 p ) and a second promoter P1 (activated by a higher concentration θ 2 p ).

Table 1: The regular domains of system [START_REF] Shachrai | Cost of unneeded proteins in E. coli is reduced after several generations in exponential growth[END_REF].

D 1 = x ∈ R 2 ≥0 : 0 ≤ x r < θ 1 r , 0 ≤ x p < θ 1 p D 2 = x ∈ R 2 ≥0 : θ 1 r < x r < θ 2 r , 0 ≤ x p < θ 1 p D 3 = x ∈ R 2 ≥0 : θ 2 r < x r ≤ max r , 0 ≤ x p < θ 1 p D 4 = x ∈ R 2 ≥0 : 0 ≤ x r < θ 1 r , θ 1 p < x p < θ 2 p D 5 = x ∈ R 2 ≥0 : θ 1 r < x r < θ 2 r , θ 1 p < x p < θ 2 p D 6 = x ∈ R 2 ≥0 : θ 2 r < x r ≤ max r , θ 1 p < x p < θ 2 p D 7 = x ∈ R 2 ≥0 : 0 ≤ x r < θ 1 r , θ 2 p < x p ≤ max p D 8 = x ∈ R 2 ≥0 : θ 1 r < x r < θ 2 r , θ 2 p < x p ≤ max p D 9 = x ∈ R 2 ≥0 : θ 2 r < x r ≤ max r , θ 2 p < x p ≤ max p .
Table 2: The focal points Φ j = (Φ j,r , Φ j,p ) for each domain D j of system [START_REF] Shachrai | Cost of unneeded proteins in E. coli is reduced after several generations in exponential growth[END_REF], in terms of the functions ϕ, η defined in [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: The bistable switch example[END_REF].

Φ j mode-p mode-r Φ 1 xp = ϕ(0, µ p , γ p ) xr = ϕ(0, µ r , γ r ) Φ 2 xr = η( xp , 0, µ p , γ r ) xp = η( xr , 0, µ r , γ p ) Φ 3 Φ 4 xp = ϕ(0, µ p , γ p ) xr = ϕ(k 1 r , µ r , γ r ) xr = η( xp , k 1 r , µ p , γ r ) xp = η( xr , 0, µ r , γ p ) Φ 5 xp = ϕ(uk 0 p , µ p , γ p ) xr = ϕ(k 1 r , µ r , γ r ) Φ 6 xr = η( xp , k 1 r , µ p , γ r ) xp = η( xr , uk 0 p , µ r , γ p ) Φ 7 xp = ϕ(0, µ p , γ p ) xr = ϕ(k 1 r + k 2 r , µ r , γ r ) xr = η( xp , k 1 r + k 2 r , µ p , γ r ) xp = η( xr , 0, µ r , γ p ) Φ 8 xp = ϕ(uk 0 p , µ p , γ p ) xr = ϕ(k 1 r + k 2 r , µ r , γ r ) xr = η( xp , k 1 r + k 2 r , µ p , γ r ) xp = η( xr , uk 0 p , µ r , γ p ) Φ 9 xp = ϕ(u[k 0 p + k 1 p ], µ p , γ p ) xr = ϕ(k 1 r + k 2 r , µ r , γ r ) xr = η( xp , k 1 r + k 2 r , µ p , γ r ) xp = η( xr , u[k 0 p + k 1 p ], µ r , γ p )
The x p equation represents RNAP with both transcription and translation regulation. The synthesis of RNAP proceeds in two steps, as it is limited by the production of two sub-units (β and β ) encoded by two different genes [START_REF] Ralling | Relative activities of the transcriptional regulatory sites in the rplKAJLrpoBC gene cluster of escherichia coli[END_REF]. Thus, a lower concentration of RNAP and ribosomes (resp., θ 1 p and θ 1 r ) are required for a basal synthesis of RNAP (uk 0 p ) whereas higher concentrations (resp., θ 2 p and θ 2 r ) are needed for the main synthesis of RNAP (u[k 0 p + k 1 p ]). During exponential phase the growth rate is constant [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF] so, for each fixed u, according to (3) the system is near a steady state ( xr (u), xp (u)), with growth rate:

μ(u) = min(µ p xp (u), µ r xr (u)) (29)
Controlling bacterial growth rate is thus equivalent to controlling the steady states by means of parameter u. From the considerations above, it follows that 0 ≤ θ 1 r ≤ θ 2 r ≤ max r and 0 ≤ θ 1 p ≤ θ 2 p ≤ max p . Therefore, the state space of each of the two modes of system (28) can be partitioned into nine regular domains (see Table 1). The threshold domains are as defined in Section 3 (not shown). The focal points of the two modes for each regular domain D j ( j = 1, ..., 9) can be calculated according to [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: The bistable switch example[END_REF] (see Table 2). Higher dimensional systems may be studied through hybrid models dedicated software such as GeneticNetworkAnalyser [START_REF] Jong | Genetic network analyzer : Qualitative simulation of genetic regulatory networks[END_REF].

To illustrate the dependence of growth rate on the control parameter u, Fig. 4 depicts a bifurcation diagram with u ∈ [0, 1] and all other parameters fixed. First, for each u, the equilibria of system (28) are calculated through [START_REF] Mukherji | Synthetic biology: understanding biological design from synthetic circuits[END_REF] and then the growth rate expression is computed using [START_REF] Sontag | Some new directions in control theory inspired by systems biology[END_REF].

Asymptotic behavior and growth rate characterization

The bifurcation diagram changes as a function of the parameters. Here we chose a set of parameters to illustrate the variety of dynamical behaviors that can be obtained with system [START_REF] Shachrai | Cost of unneeded proteins in E. coli is reduced after several generations in exponential growth[END_REF]. Fig. 4 shows six different scenarios corresponding to different combinations of equilibria and hence to different growth rates. A straightforward characterization of the different scenarios can be given in terms of the formulas in Table 2.

Proposition 2 [Absence of growth] Assume the control parameter u is such that µ p xp < µ r xr , ϕ(uk 0 p , µ p , γ p ) < θ 1 p , and ϕ(u[k 0 p + k 1 p ], µ p , γ p ) < θ 2 p . Then the only equilibrium of the system is the origin. Moreover, for any u, assume that the initial condition satisfies x p (t 0 ) < θ 1 p . Then the solution converges to the origin.

From a biological point of view, the origin corresponds to absence of bacterial growth, which happens either when u is too low to induce RNAP expression (scenario-a) or when RNAP initial condition is too low to initiate gene transcription (all scenarios in Fig. 4).

From Corollary 1, it follows that system (28) eventually remains on only one of its modes. This allows for monostability (scenarios b,c in Fig. 4) or bistability (scenarios d,e,f ). In scenarios e,f, the system is mono-mode bistable in mode-r, since the two non-zero equilibria belong to the same mode. In scenario d the system is bi-mode bistable, since both modes contribute with a positive equilibrium.

Proposition 3 [Bi-mode Bistability] Assume that the equilibria are given by Φ 9 of mode-p and Φ 5 = Φ 6 of moder. Then the two growth rate expressions are μ

= µ p ϕ(u[k 0 p + k 1 p ], µ p , γ p ) (high level) and μ = µ r ϕ(k 1 r , µ r , γ r ) (low level).
Next, we investigate the occurrence of Theorem 2 configurations. It is clear that the production rates are monotone in both variables, which means that cases as in part 3 (with periodic solutions) can be ruled out. Configurations as in part 1 or 2 may happen: damped oscillations may exist around the segment S 58 = {x : θ 1 r < x r < θ 2 r , x p = θ 2 p } between domains D 5 and D 8 , as shown next and illustrated in Fig. 5.

Proposition 4 [Damped oscillations] Assume that, in mode-r, Φ 5 = Φ 6 ∈ D 8 and Φ 8 ∈ D 5 ∪D 6 . Then the system has damped oscillations converging to an equilibrium in the focal set Φ(S 58 ).

Proof: In mode-r, the vector fields for D 5 and D 8 are of the form: ). Three trajectories are shown in dash-dotted magenta, each converging to one of the equilibria points: origin, Φ 9 , or the point given by Prop. 4. For comparison, the corresponding solutions of a differentially continuous system (2) are also shown (solid magenta curves). The parameters are as in Fig. 4 with u = 0.85.

D 5 : ẋr = k 1 r -(µ r x r + γ)x r , ẋp = uk 0 p -(µ r x r + γ)x p , D 8 : ẋr = k 1 r + k 2 r -(µ r x r + γ)x r , ẋp = uk 0 p -(µ r x r + γ)x p .
For parameters satisfying:

-γ + γ 2 + 4µ r k 1 r 2µ r < θ 2 p < -γ + γ 2 + 4µ r (k 1 r + k 2 r ) 2µ r θ 1 r < uk 0 p 2 γ 2 + 4µ r k 1 r ,
the system satisfies the conditions of Theorem 2, part 2 (with x i = x p ), at the threshold domain x p = θ 2 p . This proves the statement.

This last case leads to interesting dynamics since it permits a large range of parameters to have fairly high growth rate but always with a limited concentration of ribosomes (x r ), which are thought to have a high production cost [START_REF] Shachrai | Cost of unneeded proteins in E. coli is reduced after several generations in exponential growth[END_REF]. It can also be combined with bistability, as in scenario-f of Fig. 4.

Finally, the usefulness of the abstract CSPQ formalism is illustrated by comparison with a C 1 system, where the step functions are replaced by a function of the form (2) with h(x) = b + ax 2 /(x 2 + θ 2 ) (with a, b determined by continuity and x low = θ i j -0.3 and x high = θ i j + 0.3). In Figs. 5 and6, it is clear that the solutions of both CSPQ and C 1 systems remain qualitatively close. In addition, the new behavior characterized by Theorem 2 is also observed in the "more biological" system.

Conclusion

The class of continuous-switch piecewise quadratic (CSPQ) systems was obtained from PL systems by introducing a nonlinear term to represent dilution due to growth rate, which is modeled by different expressions depending on the internal state. The nonlinear terms introduce new dynamical behavior at the threshold domains where the vector fields are discontinuous. In the case of opposing vector fields, we have characterized different configurations which include not only the possible occurrence of sliding motion, as in PL systems, but also damped oscillations or periodic behavior.

As an application, a CSPQ model was developed to investigate the control of growth rate in E. coli cells during exponential phase. The analysis revealed several interesting biological scenarios including bistability and steady states that combine high growth rate with a low amount of ribosomes (useful for lowering protein production costs). To conclude, we believe the CSPQ formalism models properties of realistic gene expression networks with dynamics dependent on dilution. This qualitative approach is a valid starting point to help guide the construction of synthetic gene networks, by allowing biologists to select experimental conditions most likely to yield successful results.

Definition 1

 1 Given a regular domain D ∈ D rg , the point Φ(D) = ( x1 , ..., xn ) T ∈ Ω (defined by f D = d(x r )x) is called the focal point for the flow in D.

Definition 3

 3 3, 3.4 below). It is useful to extend the concept of focal point to threshold domains by defining the focal set: Let supp(D) be the (nk)-dimensional hyperplane supporting D. If D is a threshold domain of order k, then its focal set Φ(D) is

Figure 2 :

 2 Figure 2: Different configurations of dynamics at a threshold domain, x i = θ i with i = r, for a 2-dimensional PQ system. Note that x r is represented on the horizontal axis and x i on the vertical axis. The figures represent two regular domains D A and D B , separated by a threshold domain (dashed black horizontal line). Each domain contains the other's focal point (diamond symbols, Φ A and Φ B ). The dashed (resp. dot-dashed) blue curve represents the nullcline ẋi = 0 for domain D A (resp. D B ) through the respective focal point, while the solid blue curve represents a linear combination of the two vector fields: (α f D A i + (1α) f D B i )/γ i . The dashed (resp. dot-dashed) arrows represent the vector field in domain D A (resp. D B ). One trajectory is depicted for each case (solid red curves). The common parameters for (b)-(d) are: µ r = 0.5, γ r = γ i = 1, θ i = 2.1, f D A r = 2, f D A i = 3.75, f D B r = 7. The plots represent: (a) the PL case, µ r = 0, γ r = γ i = 1, f D A r = 1.2, f D A

Figure 3 :Corollary 1

 31 Figure 3: Regulatory network of the CSPQ model[START_REF] Shachrai | Cost of unneeded proteins in E. coli is reduced after several generations in exponential growth[END_REF]. The system consists of RNAP and ribosomes (RIB), encoded by their respective proxy genes rnap and rib (in fact, RNAP and RIB are made up of different subunits but here, for simplicity, we consider two proxy gene classes). The synthesis of RNAP is positively regulated by RIB, RNAP and the control parameter u while RIB is positively regulated by RNAP from the promoters P1 and P2. Corollary 1 Consider system (1) with assumption 2. Let D ∈ D rg and assume it intersects a switching surface, S D = D ∩ S r,p = / 0. Then S D contains no stable equilibria.

  μ

Figure 4 :

 4 Figure 4: A growth rate bifurcation diagram as function of the control parameter u. Black stars represent the zero steady state (present for any value of u). Blue circles (resp., red squares) represent steady states corresponding to mode-p (resp., mode-r). Red triangles represent a steady state for mode-r in the case of damped oscillations between domains D 5 and D 8 (see Prop. 4). Parameter values used: θ 1 r = 1, θ 2 r = 5, θ 1 p = 1, θ 2 p = 3, k 1 r = 6, k 2 r = 17, k 0 p = 10,k 1 p = 15, γ r = γ p = 1, µ r = 0.5, and µ p = 0.8.

Φ 8 Φ 5 Φ 9 Figure 5 :

 8595 Figure 5: Phase plane of the CSPQ system (28), corresponding to scenario-f of Fig. 4. The dashed horizontal and vertical lines represent the thresholds (θ i j ) of the system and the domains are labeled D k , k ∈ {1, . . . , 9}. The mode switching line is in solid black. The dashed red curve represents the nullcline ẋp = 0 in mode-r and also corresponds to the focal set Φ(S 58 ) (see Prop.4). Three trajectories are shown in dash-dotted magenta, each converging to one of the equilibria points: origin, Φ 9 , or the point given by Prop. 4. For comparison, the corresponding solutions of a differentially continuous system (2) are also shown (solid magenta curves). The parameters are as in Fig.4with u = 0.85.

  Figure 5: Phase plane of the CSPQ system (28), corresponding to scenario-f of Fig. 4. The dashed horizontal and vertical lines represent the thresholds (θ i j ) of the system and the domains are labeled D k , k ∈ {1, . . . , 9}. The mode switching line is in solid black. The dashed red curve represents the nullcline ẋp = 0 in mode-r and also corresponds to the focal set Φ(S 58 ) (see Prop.4). Three trajectories are shown in dash-dotted magenta, each converging to one of the equilibria points: origin, Φ 9 , or the point given by Prop. 4. For comparison, the corresponding solutions of a differentially continuous system (2) are also shown (solid magenta curves). The parameters are as in Fig.4with u = 0.85.

Figure 6 :

 6 Figure 6: Comparison of the CSPQ system (dash-dotted curves) with an analogous continuously differentiable system (2) (solid curves). The solutions as a function of time correspond to the damped oscillatory trajectories in Fig. 5.

Acknowledgements

We thank our biologist research partners D. Ropers and J. Geiselmann (IBIS team, Inria) for many fruitful discussions and their help in the development of the E. coli model.

This work was supported in part by the projects GeMCo (ANR 2010 BLANC020101), ColAge (Inria-Inserm Large Scale Initiative Action), RESET (Investissements d'Avenir, Bioinformatique), and by the LABEX SIG-NALIFE (ANR-11-LABX-0028-01). The material in this paper was in part presented at the 9th IFAC Symp. Nonlinear Control Systems (NOLCOS), 2013, Toulouse, France.