domain: info.info-au math spi sdv

text:
Stefano Casagranda
email: stefano.casagranda@inria.fr

Delphine Ropers
email: delphine.ropers@inria.fr

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

Model reduction and process analysis of biological models

Keywords: Dynamical Systems, Biological Networks, Model Reduction, Process Analysis 1

Understanding the dynamical behavior of biological networks is complicated due to their large number of components and interactions. We present a method to analyse key processes for the system behavior, based on the a priori knowledge of the system trajectory and the simplification of mathematical models of these networks. The method consists of the model decomposition into biologically meaningful processes, whose activity or inactivity is evaluated during the time evolution of the system. The structure of the model is reduced to the core mechanisms involving active processes only. We assess the quality of the reduction by means of global relative errors and apply our method to two models of the circadian rhythm in Drosophila and the influence of RKIP on the ERK signaling pathway.

I. INTRODUCTION

Mathematical modeling and simulation tools are determinant in understanding how adaptation of living organisms to environmental cues results from large networks of metabolites, RNAs, proteins, and their mutual interactions. Larger and larger kinetic models of cellular networks are nowadays published, as a result of decades of work in biology, recent advances in high throughput technologies, and progress in modeling and parameter estimation approaches (see [START_REF] Bettenbrock | A quantitative approach to catabolite repression in escherichia coli[END_REF] and [START_REF] Kuepfer | Ensemble modeling for analysis of cell signaling dynamics[END_REF] for an example). The large size of these models and their non linearity due to complex feedback loops make their dynamical analysis rather difficult. More specifically, it is extremely difficult to relate the global behaviour of the system to the functioning of specific cellular processes (e.g. RNA transcription, protein phosphorylation, or complex formation), while this gain of knowledge is crucial to identify what are the key cellular processes for the environmental adaptation and when they are at play. Efforts in this direction rely on the reduction of the system complexity, as reduced models can be more easily analysed. Beside testing model robustness, sensitivity analysis is hence used to detect model parameters that are not influential for the system dynamics and that can be removed (i.e. [START_REF] Degenring | Sensitivity analysis for the reduction of complex metabolism models[END_REF]). When different time scales are present in a model, quasisteady-state approximation is often applied to reduce the system dimension and replace the fast variables by algebraic expressions. Nevertheless, the resulting differential algebraic system can be still complex to analyse [START_REF] Segel | The quasi-steady-state assumption: a case study in perturbation[END_REF]. Finally, phase plane analysis allows analysing the asymptotic behaviour of kinetic models, but its use is restricted to small size systems [5, p. 21]. Alternative approaches consist in simplifying the mathematical functions describing the molecular processes. For instance, piece-wise affine differential equations approximate by step functions the sigmoidal functions describing cooperativity in the regulation of gene expression. The dynamics of this simplified system can be easily analysed by means of state transition graphs [START_REF] Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF]. However, these simplifications are generally restricted to models of gene expression and are more difficult to apply to model other types of network [START_REF] Baldazzi | The carbon assimilation network in escherichia coli is densely connected and largely sign-determined by directions of metabolic fluxes[END_REF]. We tackle the problem in this study, by developing a mathematical and numerical approach for analysing the contribution of cellular processes to the dynamics of an ordinary differential equation system. The method allows i) to identify, at different times of the dynamics, the mechanisms that are active and contribute most to the system dynamics and ii) to reduce the kinetic model into a series of new models reduced to the core mechanisms, neglecting the inactive processes. The method is general and can be easily applied to any ODE model of biological system. To some extent, our work is similar to the model reduction approach of L. Petzold and W. Zhu [START_REF] Petzold | Model reduction for chemical kinetics: An optimization approach[END_REF] applied to chemical kinetics, which consists in decomposing ODE models into a sum of elementary reactions and finding simplified models with a reduced number of reactions. Removing part of the network is also considered by Apri et al. in [START_REF] Apri | Complexity reduction preserving dynamical behavior of biochemical networks[END_REF], where parameters that are not important for the system dynamics are removed from the model. We apply our method to two models with different dynamical behaviors. In the model of the circadian rhythm in Drosophila [START_REF] Leloup | A model for circadian rhythms in drosophila incorporating the formation of a complex between the per and tim proteins[END_REF], the system dynamics reaches a regime of sustained periodic oscillations, whereas this dynamics is simpler in the model of the influence of Raf Kinase Inhibitor Protein on the Extra-cellular Signal Regulated Kinase [START_REF] Kwang-Hyun | Mathematical modeling of the influence of rkip on the erk signaling pathway[END_REF]. The solutions quickly reach a steady-state. Section 2 of this paper describes the methodology of the model simplification. The method is applied to the circadian and ERK models in Section 3 and 4, respectively. We conclude in Section 5, by giving some future perspectives.

II. METHODOLOGY

We evaluate the different processes along the trajectory of the model, in order to study their relative influence during the evolution of the system. As such, this technique takes naturally into account the order of magnitude of the variables. Consider a general ODE model of biological network:

ẋ = f (x, p) (1)
and an ODE system which approximates the previous model:

ẋr = g (x r , p r ) (2)
where x = (x 1 , x 2 , ..., x n ) R n and x r = (x r 1 , x r 2 , ..., x r n ) R n represent, respectively, the vectors of concentration of components in the original and reduced models. p R m and p r R k with k ≤ m are respectively the vectors of parameters of the original and reduced models, where each parameter has a precise and known value. The vector of initial conditions is x 0 = x r 0 , where x 0 R n is also known. The time interval on which we want to study the processes is [0, T ] (it can be either a short or a long period). The basic idea of our method is based on the following theorem: if the vector fields of two systems are close (f (x) ≈ g(x)), then the solutions of the original and approximated systems are close during some time interval under the assumptions on the Lipschitz conditions listed in [START_REF] Khalil | Nonlinear systems[END_REF]p. 79,Th. 2.5]. We first identify the processes in the model equations.

ẋi = j f ij (x, p) (3)
where f ij represents the j th process involved in the evolution of the i th coordinate where i = 1, ..., n. Generally a biological model is already decomposed in terms of a linear combination of processes: they evolve dynamically over time and contribute to the time evolution of one or more variables. Let us take as an example the 8 th coordinate of the Goldbeter and Leloup's model [START_REF] Leloup | A model for circadian rhythms in drosophila incorporating the formation of a complex between the per and tim proteins[END_REF] (see Section 3 and the model equations in Appendix).

dT 2 dt = V 3T T 1 K 3T + T 1 -V 4T T 2 K 4T + T 2 -k 3 P 2 T 2 (4) +k 4 C -v dt T 2 K dT + T 2 -k d T 2
The equation represents the variation of the concentration of the double phosphorylated timeless protein. It contains six different processes, each of which with a specific biological meaning. They can be positive or negative.

ẋ8 = f 8,1 -f 8,2 -f 8,3 + f 8,4 -f 8,5 -f 8,6 (5)
where

f 8,1 = V 3T T1 K 3T +T1 , ..., f 8,6 = k d T 2 .
In order to compare the influence of the different processes f ij in the time evolution of each variable x i , we have to find a criterion to compare the processes between themselves. It is possible, for example, to compare their absolute function (|f ij (x, p)|), to scale it by the i th initial condition (

|fij (x(t),p)| x0i
) or to scale it by the solution of the

i th ODE ( |fij (x(t),p)| x(t)i
). In this work we choose to associate a relative weight to each process with the purpose to make it dimensionless:

W ij (t, p) = |f ij (x(t), p)| j |f ij (x(t), p)| (6)
where 0 ≤ W ij (t, p) ≤ 1 and j W ij (t, p) = 1. Definition: Let the continuous function f ij (x(t), p) be the j th process of ẋ(t) i in t [0, T ] and let the threshold

θ ⊂ R θ [0,1]. We call a process f ij (x(t), p) always inactive when W ij (t, p) < θ ∀ t [0,T]. We call a process f ij (x(t), p) inactive at time t when W ij (t, p) < θ. We call a process f ij (x(t), p) active at time t when W ij (t, p) ≥ θ.
The switching time for a process f ij (x(t), p) is the time t s ij when W ij (t, p) = θ. A process can have 0, 1, ..., s switching times. The switching time set S i for the i th variable contains all the switching times t s ij where j = 1, .., k. The global switching time set S is the union of all S i . The relative weights evolve with time: for example a process can be active at the beginning of the system evolution and becomes inactive at the end. Fig. 1(A) shows the dynamics of the absolute values of the processes in (4) for t [0,24] (hours). In Fig. 1(B) is shown the dynamics of the absolute weights. First step: eliminate the always inactive processes. We are looking for a function g(x r ) which approximates the function f (x) and with N g P < N f P , where N P is the number of processes in the ODE system. The threshold θ value must be chosen between the range [0,1]: a low threshold should be chosen to avoid neglecting too many processes and changing significantly the dynamics of the new system with respect to the original one. In our work a threshold θ = 0.1 is set. After having applied this procedure for every kinetic equation, we obtain the reduced model where:

if W ij (x(t), p) < θ ∀ t [0,T] then f ij (x(t), p) = 0; if not, g ij (x(t), p) = f ij (x(t), p).
To test the quality of the reduced model g(x r ), we numerically compute the global relative error between the original and the reduced models for each variable: it is defined for the i th variable as

|x i (t) -x r i (t)|dt |x i (t)|dt (7)
where x i (t) and x r i (t) are respectively the solutions of the original and reduced systems. Taking (4) as an example, we see that the processes f 8,3 , f 8,4 , f 8,6 are always under θ and that they can be neglected. The processes f 8,1 and f 8,5 are always over θ, while f 8,2 crosses the threshold twice and its switching times (t 1 8,2 = 0.2 h and t 2 8,2 = 19.8 h) are collected in the set S 8 . Second step: visualize the dynamical activity of the remaining processes. We qualitatively study the evolution of the intensity of the process activity in the reduced model g(x r ). We use in this paper two graphic tools to summarize this information in one picture: a three dimensional bar graph and a heat map. Third step: create a succession of sub-models. The time period in which the system evolves can be split into time intervals using the switching times t i previously grouped in the set S, sorted in ascending order: a succession of sub-models can be created, each one valid for a certain time window (∆t i = t i -t i-1 ), which contains the core mechanisms in that period of time. To avoid having a large chain of sub-models, we reduce the number of time windows, by compacting with a clustering method the switching times that are close. We obtain a reduced number of time windows based on the reduced number of switching times: ∆t r i = t r i -t r i-1 where N ∆t r i << N ∆ti . A good method to compact them in few clusters is the k-means clustering [START_REF] Kanungo | An efficient k-means clustering algorithm: Analysis and implementation[END_REF], which is really easy to compute in one dimension (time axis). Given our global switching time set (sorted in ascending order) S = [t 1 , t 2 , ..., t n ], this method leads to group the n switching times into k (≤ n) clusters C={C 1 , C 2 , ..., C n } so as to minimize the within-cluster sum of squares (WCSS):

argmin C k i=1 t Ci ||t -µ i || 2 (8)
where µ i is the mean of the switching times in C i . Following this operation, the cluster C i is replaced by the mean µ i , which is going to be the new switching time t r i , and we assume that all the processes that were going to switch their state in the cluster C i are going to switch together at time t r i . There is no precise rule to choose k, but its value can be related with the difference between the maximum and minimum number of active processes during the time evolution of the system: if the difference is low, k should be chosen low as well. We finally end up with a small chain of sub-models that can be biologically interpreted.

III. MODEL FOR CIRCADIAN RHYTHMS IN DROSOPHILA

A. Description

We consider the model in [START_REF] Leloup | A model for circadian rhythms in drosophila incorporating the formation of a complex between the per and tim proteins[END_REF] given by an ODE system of dimension 10. The model describes the circadian oscillations of the proteins timeless (TIM) and period (PER), which involve a negative feedback loop: the double phosphorylated forms of these proteins can be degraded in the cytoplasm or form the PER-TIM complex which, following its transport to the nucleus, inhibits the transcription of the tim and per genes and the subsequent accumulation of their mRNAs and proteins. The light sets the period of the oscillations to 24 hours precisely, by increasing the velocity of the degradation of the double phosphorylated form of TIM (from 2 nM/h to 4 nM/h in the model). Since the same oscillatory behavior is repeated every day, we focus our analysis on 24 hours only and use the model parameters given in Figures 2 and4 of [START_REF] Leloup | A model for circadian rhythms in drosophila incorporating the formation of a complex between the per and tim proteins[END_REF].

B. Model reduction

After having decomposed every kinetic equation of the model into a set of processes that we consider to be biologically relevant, we simulate their relative weight as shown before for (4) in Fig. 1. Every time that the weight of a process crosses the threshold θ, the corresponding process changes its state. The timing of the process activation/inactivation in each equation and the active/inactive processes along time is conveniently shown in a Temporal Process Map, as displayed in Fig. 2. From this analysis, we obtain a reduced model by neglecting the processes that are always inactive (under the threshold θ = 0.1 during the whole time). We hence neglect 18 of 44 processes: f 1,3 , f 2,3 , f 2,4 , f 3,2 , f 3,4 , f 3,5 , f 4,2 , f 4,6 , f 5,3 , f 6,3 , f 6,4 , f 7,2 , f 7,5 , f 8,3 , f 8,4 , f 8,6 , f 9,5 , f 10,3 . To test the quality of the reduced model, we calculate the global relative error between the original and the reduced models for each variable using [START_REF] Baldazzi | The carbon assimilation network in escherichia coli is densely connected and largely sign-determined by directions of metabolic fluxes[END_REF]. The results in Table I show a good match between the original model and the reduced one.

C. Qualitative tool: Heat Map

In the previous analysis we used a Boolean approach to study the process activity or inactivity by means of the Temporal Process Map. We also have considered alternative graphical representations in this paper, since we also want to qualitatively understand, in one graph, the change over time of the intensity of the activity/inactivity of important biological processes. For the circadian model, we chose an Heat Map as a qualitative approach, where the individual values contained in a matrix are represented using colors. In the Heat Map that we generate, values along the rows are standardized, so that processes assume a red color (active) if their value is above the mean, black if their value is equal to the mean and green (inactive) if their value is below the mean of a column across all rows. If the red (or green) color is lighter, it means that the process is more active (or inactive). The Heat Map for the reduced model g(x r ) is shown in Fig. 3.

D. Creation of sub-models based on time windows

The Temporal Process Map is not only useful to identify negligible processes for model reduction, it can be used also to study the evolution of important (non-negligible) processes and create a succession of sub-models containing the core mechanisms, each one being valid for a certain time window. the number of time windows is 10. We use the k-means clustering to compact together more switching times: we choose k = 2 because the difference between the maximum and minimum number of active processes is five (low value). The steps are shown in Fig. 5. We have thus created 2 submodels: the first one (valid from 0 to 1.96 h and from 17.8 h to 24 h) coincides with the reduced model and the second one (valid from 1.96 h to 17.8 h) suppresses also the processes f 4,3 , f 4,4 , f 7,4 , f 8,2 , f 9,1 . In Fig. 6 is shown a Flux Map for the two sub-models, where red lines represent active processes in that time-window, while inactive ones are represented in black. The first submodel essentially corresponds to a situation of night time and the second one, of day time. The results in Table II show a

IV. MODEL FOR THE INFLUENCE OF RKIP ON THE ERK SIGNALING PATHWAY A. Description

We consider the model in [START_REF] Kwang-Hyun | Mathematical modeling of the influence of rkip on the erk signaling pathway[END_REF] of dimension 11. The ERK signaling pathway controls important cellular phenomena like proliferation or differentiation. The model describes the inhibition of the activation of RAF by RKIP, which regulates the ERK signaling pathway. We use in our analysis the parameters of [START_REF] Kwang-Hyun | Mathematical modeling of the influence of rkip on the erk signaling pathway[END_REF] and initial values in [START_REF] Petrov | Reduction of nonlinear dynamic systems with an application to signal transduction pathways[END_REF].

Since no switches occur during the simulation (due to the strong activity/inactivity boundary separation and the absence of process weights crossing the threshold θ), we will show neither the Temporal Process Map nor the creation of submodels.

B. Model reduction

We simulate the relative weight of the processes for every ODE as we did for the previous model. We are able to neglect 12 of 34 processes:

f 1,2 , f 2,2 , f 3,2 , f 3,4 , f 4,2 , f 5,3 , f 6,3 , f 7,2 , f 8,2 , f 9,2 , f 10,2 , f 11,2 .
The results in Table III show a really good match between the original model and the reduced one with a low global error. In addition to our study, other model reduction approaches with different goals have been applied to the ERK model in [START_REF] Kwang-Hyun | Mathematical modeling of the influence of rkip on the erk signaling pathway[END_REF]: the quasi-steady-stateapproximation used by Petrov et al. in [START_REF] Petrov | Reduction of nonlinear dynamic systems with an application to signal transduction pathways[END_REF] and the automatic complexity analysis by Lebiedz et al. in [START_REF] Lebiedz | Automatic complexity analysis and model reduction of nonlinear biochemical systems[END_REF]. The work of [START_REF] Petrov | Reduction of nonlinear dynamic systems with an application to signal transduction pathways[END_REF] is concerned with the separation of variables with fast dynamics from those with slow dynamics with respect to a time scale and uses a mathematical scaling. The components of the ODEs related to the fast variables are expressed by algebraic equations. Analysing the resulting reduced model allowed the authors of the study to conclude that the variable m 4 (the complex Raf-1*/RKIP/ERK-PP) has the biggest influence on the system when it approaches its quasi-steady state. The work of [START_REF] Lebiedz | Automatic complexity analysis and model reduction of nonlinear biochemical systems[END_REF] uses two different methods: the first one combines dynamic sensitivity analysis with singular value decomposition to find a minimal dimension of the model and the second one permits to reduce actually the dimension of the model and determine the variables which contribute more to the full dynamics of the system (variables m 5 , m 8 , m 11 ). Our analysis is different in the sense that we are not interested in which variable gives the bigger contribution to the dynamics of the full system but which mechanisms (processes) give the highest contribution to the dynamics of the variable: we study the influence of the processes. It is a reduction method that does not change in general the dimension and the biological structure of the system. For instance, the model in [START_REF] Kwang-Hyun | Mathematical modeling of the influence of rkip on the erk signaling pathway[END_REF] is composed by two types of processes: association processes, where two or more proteins combine together to form a complex, and dissociation processes which correspond to the reverse mechanism. While the association processes are always active during the time evolution of the system, most of the dissociation processes can be considered negligible: it means that association processes play the bigger role in the dynamic for every variable of the system.

C. Qualitative tool:3-D bar Graph

We chose this time a three-dimensional bar graph as a qualitative approach to show the evolution of the intensity of the activity for every process (Fig. 7). The height of the column gives the intensity of each process activity, while the color code indicates the weight of the process. In this graph we put the processes of the reduced model g(x r ) during a one-second dynamics discretized in six time intervals. Fig. 7 shows that the most active process during the system dynamics is the seventh process, the association of Raf-1*/RKIP and ERK-PP to form the Raf-1*/RKIP/ERK-PP 3D bar graph: the x-axis represents the processes (the black arrow starts from the first process and ends with the last one), the y-axis corresponds to the time evolution and the z-axis, to the value of the weight. complex. This process strongly influences the variable m 4 , which was found by Petrov et al. to be the most important one for the system dynamics.

V. CONCLUSIONS

We have presented in this paper a method that allows to reduce biological models of high dimension and to analyze their core mechanisms and how these evolve dynamically, provided that their initial conditions and parameter values are known. Reduced models and sub-models obtained by this approach can be refined, by calibrating the reduced vector of parameters so as to obtain models that better approximate the original one. In a future step, Parameter Sensitivity Analysis will be applied to the reduced models and sub-models to test their robustness. Biological models of higher dimension with some parameter uncertainty (e.g. in some intervals) will be also analyzed in future works to study the consequences of this uncertainty on the dynamics of the sub-models.

APPENDIX

Full Drosophila model (see [START_REF] Leloup | A model for circadian rhythms in drosophila incorporating the formation of a complex between the per and tim proteins[END_REF])

dM P dt = v sP K n IP K n IP + C n N -v mP M P K mP + M P -k d M P dP 0 dt = k sP M P -V 1P P 0 K 1P + P 0 + V 2P P 1 K 2P + P 1 -k d P 0
= V 1P P 0 K 1P + P 0 -V 2P P 1 K 2P + P 1 -V 3P P 1 K 3P + P 1 +V 4P P 2 K 4P + P 2 -k d P 1 dP 2 dt = V 3P P 1 K 3P + P 1 -V 4P P 2 K 4P + P 2 -k 3 P 2 T 2 +k 4 C-V dP P 2 K dP + P 2 -k d P 2 dM T dt = v sP K n IT K n IT + C n N -v mT M T K mT + M T -k d M T dT 0 dt = k sT M T -V 1T T 0 K 1T + T 0 + V 2T T 1 K 2T + T 1 -k d T 0 dT 1 dt = V 1T T 0 K 1T + T 0 -V 2T T 1 K 2T + T 1 -V 3T T 1 K 3T + T 1 +V 4T T 2 K 4T + T 2 -k d T 1 dT 2 dt = V 3T T 1 K 3T + T 1 -V 4T T 2 K 4T + T 2 -k 3 P 2 T 2 +k 4 C -V dT T 2 K dT + T 2 -k d T 2 dC dt = k 3 P 2 T 2 -k 4 C -k 1 C + k 2 C N -k dC C dC N dt = k 1 C -k 2 C N -k dN C N
Reduced Drosophila model and first sub-model

K 1T + T 0 dT 1 dt = V 1T T 0 K 1T + T 0 -V 3T T 1 K 3T + T 1 + V 4T T 2 K 4T + T 2 dT 2 dt = V 3T T 1 K 3T + T 1 -V 4T T 2 K 4T + T 2 -V
K 1T + T 0 dT 1 dt = V 1T T 0 K 1T + T 0 -V 3T T 1 K 3T + T 1 dT 2 dt = V 3T T 1 K 3T + T 1 -V dT T 2 K dT + T 2 dC dt = -k 4 C -k 1 C + k 2 C N dC N dt = k 1 C -k 2 C N
Full model ERK model (see [START_REF] Kwang-Hyun | Mathematical modeling of the influence of rkip on the erk signaling pathway[END_REF])

Fig. 1 .

 1 Fig. 1. (A) Example of simulation of the processes of the 8th coordinate of the ODE system in [10] and (B) of the relative weights with the threshold θ set at 0.1. Three processes are always inactive during the time evolution of the system.

Fig. 2 .

 2 Fig. 2. Temporal Process Map of the process activity for the ODE model in Drosophila [10]. Black: active state. White: inactive state. Times t i are, in the order, 0.2, 0.5, 1.3, 1.8, 6, 12.3, 15.8, 19.6, 19.8, and 21.5 h. The total duration of the simulation is 24 h.

Fig. 5 (

 5 Fig. 3. Heat Map applied to the Drosophila Model.

Fig. 4 .

 4 Fig. 4. Evolution of the number of active processes as a function of time. The function increases or decreases at the switching time t i in the time interval of the system evolution.

Fig. 5 .

 5 Fig. 5. (A)shows the based-event grid: every time window is delimited by the switching time t i in the time interval of the system evolution (upper picture). In (B) the algorithm subsequently chooses the membership of every t i in the Cluster C 1 or C 2 to minimize the WCCS expressed by (8) (middle picture). In (C) the cluster is then replaced with its centroid (in this case the mean) that will be the new approximate switching time. (D) shows the approximation of Fig.4. The number of time windows becomes 2.

Fig. 6 .

 6 Fig. 6. Upper panel: Flux Map for the sub-model valid from 0 to 1.96 h and 17.8 h to 24 h. Lower panel: Flux Map for the sub-model valid from 1.96 h to 17.8 h. Black: inactive processes. Red: active processes.

P 1 K 3P + P 1 dP 2 dt = V 3P P 1 K 2 K dP + P 2

 112122 3P + P 1 -k 3 P 2 T 2 + k 4 C -V dP P = k sT M T -V 1T T 0

dT T 2 K dT + T 2 dC dt = k 3 P 2 T 2 -k 4 C -k 1 C + k 2

 2322412 = k sT M T -V 1T T 0

dm 1 dt = -k 1 m 1 m 2 + k 2 m 3 + k 5 m 4 dm 2 dt = -k 1 m 1 m 2 + k 2 m 3 + k 11 m 11 dm 3 dt = k 1 m 1 m 2 -k 2 m 3 -k 3 m 3 m 9 + k 4 m 4 dm 4 dt = k 3 m 3 m 9 -k 4 m 4 -k 5 m 4 dm 5 dt = k 5 m 4 -k 6 m 5 m 7 + k 7 m 8 dm 6 dt = k 5 m 4 -k 9 m 6 m 10 + k 10 m 11 dm 7 dt = -k 6 m 5 m 7 + k 7 m 8 + k 8 m 8 dm 8 dt = k 6 m 5 m 7 -k 7 m 8 -k 8 m 8 dm 9 dt = -k 3 m 3 m 9 + k 4 m 4 + k 8 m 8 dm 10 dt = -k 9 m 6 m 10 + k 10 m 11 + k 11 m 11 dm 11 dt = k 9 m 6 m 10 -k 10 m 11 -k 11 m 11 Reduced ERK model dm 1 dt = -k 1 m 1 m 2 + k 5 m 4 dm 2 dt = -k 1 m 1 m 2 + k 11 m 11 dm 3 dt = k 1 m 1 m 2 -k 3 m 3 m 9 dm 4 dt = k 3 m 3 m 9 -k 5 m 4 dm 5 dt = k 5 m 4 -k 6 m 5 m 7 dm 6 dt = k 5 m 4 -k 9 m 6 m 10 dm 7 dt = -k 6 m 5 m 7 + k 8 m 8 dm 8 dt = k 6 m 5 m 7 -k 8 m 8 dm 9 dt = -k 3 m 3 m 9 + k 8 m 8 dm 10 dt = -k 9 m 6 m 10 + k 11 m 11 dm 11 dt = k 9 m 6 m 10 -k 11 m 11

 12322332394945476410778878994101011111011111222324954647787991010111011

TABLE I GLOBAL

 I RELATIVE ERRORS FOR THE REDUCED DROSOPHILA MODEL

	Variable	G. Rel. Err. (%)
	Period mRNA	10.82
	Total Period Protein	3.70
	Timeless mRNA	7.54
	Total Timeless Protein	5.80
	Complex	6.08
	Nuclear Complex	4.56

TABLE III GLOBAL

 III RELATIVE ERRORS FOR THE REDUCED ERK MODEL

	Variable	G. Rel. Err. (%)
	Raf-1*	0.0012
	RKIP	0.106
	Raf-1*/RKIP	6.172
	Raf-1*/RKIP/ERK-PP	0.0036
	ERK-P	9.07
	RKIP-P	0.138
	MEK-PP	0.0056
	MEK-PP/ERK	0.0148
	ERK-PP	0.0013
	RP	0.00049
	RKIP-P	0.03451
	Fig. 7.

ACKNOWLEDGMENT

We acknowledge Conseil Régional PACA and the Investissements d'Avenir Bio-informatique programme under project RESET (ANR-11-BINF-0005) for funding the PhD thesis of S. Casagranda. We also thank the research program LABEX SIGNALIFE.


A Merlen
email: merlen@univ-tln.fr

J Plathier

A Ruediger
email: ruediger@emt.inrs.ca

Near Field optical image of a gold surface: a luminescence study

This paper addresses recent experimental findings about luminescence of a gold tip in near-field interaction with a gold surface. Our electrochemically etched gold tips show a typical, intrinsic luminescence that we exploit to track the plasmon resonance modeled by a Lorentzian oscillator. Our experimental device is based on a spectrometer optically coupled to an atomic force microscope used in tuning fork mode. Our measurements provide evidence of a strong optical coupling between the tip and the surface. We demonstrate that this coupling strongly affects the luminescence (Intensity, wavelength and FHWM) as a function of the tip position in 2D maps. The fluctuation of these parameters is directly related to the plasmonic properties of the gold surface and is used to qualify the optical near field enhancement (which subsequently plays the predominant role in surface enhanced spectroscopies) with a very high spatial resolution (typically around 20 nm). We compare these findings to the independently recorded near-field scattered elastic Rayleigh signal.

Introduction

Noble metal surfaces play a predominant role in plasmon-based surface enhanced spectroscopy techniques. Their investigation requires characterization tools with a very high spatial resolution, traditionally reserved to electron microscopy techniques. In the particular case of metallic nanostructures the strong near field electromagnetic enhancement has a typical range of a few nanometers [START_REF] Novotny | Principles of Nano-Optics[END_REF] . As it plays a key role in surface enhanced spectroscopy [START_REF] Moskovits | Surface Enhaced SPectroscopy[END_REF][START_REF] Moskovits | Surface-Enhanced Raman Spectroscopy: A Brief Retrospective[END_REF] (SERS, SEF SEIRA etc.), the determination of its spatial distribution around metallic nanostructures is required to achieve a better understanding of the mechanisms involved in any of these techniques. For this reason a variety of experimental methods [START_REF] Vogelgesang | Real-Space Imaging of Nanoplasmonic Resonances[END_REF] have been adapted for the investigation of the near field. Some use electronic microscopy (e.g. EELS [START_REF] Liang | Asymmetric Silver "Nanocarrot" Structures: Solution Synthesis and Their Asymmetric Plasmonic Resonances[END_REF] , Cathodolumuniscence [START_REF] Gómez-Medina | Mapping Plasmons in Nanoantennas via Cathodoluminescence[END_REF] ). This technique has demonstrated several benefits (high resolution for both spectral and spatial measurements) but it has also some major limitations such as the requirement for conductive samples, high vacuum environments, transparency for electron beam for EELS measurements.In addition the exact nature of the interaction between electrons and plasmons is still under debate even if its understanding has made substantial progresses recently [START_REF] Kociak | Mapping Plasmons at the Nanometer Scale in an Electron Microscope[END_REF] . For these reasons, optical measurements are generally privileged. However, it was not until the end of the 1990s that optical scanning probe techniques bypassed the diffraction limit of classical optical measurements. Scanning near field Optical Microscopy (SNOM) is divided in two main approaches: aperture and apertureless. Aperture SNOM has already demonstrated its efficiency for the study of metallic plasmonic structures [START_REF] Imura | Photoluminescence from Gold Nanoplates Induced by near-FIeld Two-Photon Absorptio[END_REF] but this technique has a weak optical collection leading to a tradeoff between signal strength and spatial resolution. For this reason, apertureless SNOM (sometimes referred to as aSNOM) has been proposed: the probe is a tip (AFM or STM tip) whose role is to transform an undetectable near field signal into a measurable far field one. For the study of plasmonic nanostructures specific tips have been developed [START_REF] Koglin | Material Contrast in Scanning near-FIeld Optical Microscopy at 1-10 Nm Resolutio[END_REF][START_REF] Heimel | SNOM/STM Using a Tetrahedral Tip and a Sensitive Current-to-Voltage Converter[END_REF] . The nature of the optical interaction between the tip and the sample is the key point in aSNOM [START_REF] Esteban | Tip-Substrate Interaction in Optical near-Field Microscopy[END_REF][START_REF] Porto | Resonance Shift Effects in Apertureless Scanning near-Field Optical Microscopy[END_REF] . To ensure a strong far field signal, the coupling between the tip and the sample needs to be strong. This is the case with plasmonic tips: the tip-sample interaction can generate highly enhanced and confined (a few nm) near fields leading to a measurable far field signal above the confocal background. High spatial resolution can thus be obtained and measurement down to even single orbital level has already been reported [START_REF] Hamann | Enhanced Sensitivity near-Field Scanning Optical Microscopy at High Spatial Resolution[END_REF][START_REF] Zhang | Chemical Mapping of a Single Molecule by Plasmon-Enhanced Raman Scattering[END_REF] . aSNOM has been used to perform measurements of elastic-scattering phenomena and also inelastic scattering such as Raman (Tip Enhanced Raman Spectroscopy [START_REF] Stadler | Developments in and Practical Guidelines for Tip-Enhanced Raman Spectroscopy[END_REF] ), fluorescence [START_REF] Mauser | Tip-Enhanced near-Field Optical Microscopy[END_REF] , second harmonic detection [START_REF] Zayats | Apertureless Scanning near-Field Second-Harmonic Microscopy[END_REF] , etc. The high spatial resolution of aSNOM is thus steadily transforming into a table-top alternative to electron microscopy in the field of plasmionics with the additional virtue of flexible environments. Technically, even liquid media are accessible.

Even aSNOM techniques have several remaining challenges. While the alignment of the tip and the laser focus has been improved through persistent progress in instrumentation, the signal to noise ratio relies on the availability of excellent near field antennae. This is particularly true when measurements are performed on non-transparent samples, which is the case for nanoelectronics and metallic surface. For these samples, optical side illumination and collection is ideal [START_REF] Nicklaus | Tip Enhanced Raman Spectroscopy with Objective Scanner on Opaque Samples[END_REF] and requires a dedicated AFM head to accommodate a long working distance objective with high optical aperture while at the same time minimizing mechanical loops. For chemical and structural mapping, i.e. the acquisition of hyperspectral images, the total time is a relevant factor and scales linearly with the quality factor of the near field antenna. For instance, for a 128 × 128 pixel image, an acquisition of 1 s (which is already rather low for many optical signals, in particular inelastic scattering) it will take more than 4 hours to get the full mapping. At this time-scale, thermal drift of piezoelectric systems needs to be taken into consideration. This drift still remains a limitation in the final possible spatial resolution in SNOM measurements: even if the theoretical resolution can be very high, the final experimental resolution will be fixed by the number of points to limit the total acquisition time.

In most aSNOM studies, the signal originates from the sample itself and is enhanced by the tip. In this configuration, the tip acts like a nanoantenna, which enhances the optical response of the sample. We now introduce the inverse configuration: the optical signal from the tip itself is enhanced by the plasmonic properties of the sample to provide a position-dependent signal: the intensity of the signal is used to map the plasmonic properties of the sample. In such a configuration, the spatial resolution can be better than 10 nm and relies on the radius of the tip. The main challenge of this approach is the realization of a tip with a measureable optical signal. Functionalized tip with fluorescent or Raman active [START_REF] Bortchagovsky | Possibilities of Functionalized Probes in Optical near-FIeld Microscop[END_REF] molecules have been proposed [START_REF] Aigouy | Field Distribution on Metallic and Dielectric Nanoparticles Observed with a Fluorescent near-Field Optical Probe[END_REF] but their realization remains problematic and their final size can strongly affects the spatial resolution.

In this article, we therefore introduce an alternative approach: the luminescence from the gold tip itself is strongly enhanced by the interaction with a gold surface. We have developed a reproducible experimental procedure for the preparation of luminescent gold tips with a very small radius (typically a few nm). In such a configuration, a strong signal is measured with short integration times and without compromising the spatial resolution. The final contrast in the optical mappings originates from the variation of the tip surface coupling due to the plasmonic properties of the substrate. Such an approach offers different advantages: a relatively simple optical device without complex post-treatment analysis, high spatial resolution, short acquisition times and limited topographic artifacts compared to the elastically scattered signal.

Method

The optical near filed measurements are performed with an AIST-NT OmegaScope 1000 using gold tips and He-Ne laser excitation (632.8 nm, TEM00). A schematic view of this experimental device can be seen in Figure 1. The beam power was set to 300 µW and its polarization is aligned parallel to the tip axis. The tip is illuminated using a 0.7 NA objective (Mitutoyo M Plan Apo 100x) mounted on a piezoelectric stage oriented at 65° form the tip axis. The piezoelectric objective scanner has a range of 15 μm along the optical axis and 30 μm on the other directions. That allows for a precise alignment of the tip apex with the focal spot of the objective by directly tracking the luminescence emission of the tip. Optical and topographic maps were performed using the tip in a shear force microscopy mode. The optical signal was collected by the same objective. Elastically back-scattered light is simultaneously detected with a Photomultiplier whereas luminescence is analyzed using a 600 l/mm grating coupled to a cooled CCD detector (Andor DU420A-BU). In this configuration, only the sample moves during the measurements in order to assure the same illumination condition at each points and the same alignment of the tip apex with the focal spot of the objective. The main advantage of this experimental device is that it allows measurements with any kind of samples (transparent or not) compatible with AFM. Tips were fabricated by an electrochemical process. A 100 µm gold wire (Goodfellow, Purity 99.99+ Annealed) was immersed by 1 mm into a hydrochloric acid solution (Sigma Aldrich, 37%) and connected as anode to an arbitrary waveform generator (Keithley 3390). The cathode consisted of a gold ring. A pulsed voltage (7.5 V peak to peak, 30 µs pulses, 3 kHz) was applied to the tip with an offset of 0.5 V. Then, the tip was cut to a length of about 300 µm and glued to a tuning fork (Abracon Corporation, AB38T-32.768KHZ) using a UV-curable glue. This process gives smooth tip with a curvature radius down to 10 nm and an aspect ratio estimated to 30:1when the process is highly optimized (see electronic supplementary information for a typical SEM image of a tip).Their final optical properties strongly depend on their exact shape and final size and even a slight variation of these parameters can affect luminescence. Some tips give rise to a very low luminescence whereas some others provide a very intense signal. For our measurement we have privileged the tip with the most intense signal.

Near field measurements were performed on a gold surface (Arrandee TM ) previously prepared by flame annealing. The optical near field study of such a sample appears as a major challenge: it is not transparent and it has a rather low intrinsic signal compared for instance with plasmonic metallic nanostructures [START_REF] Imura | Near-Field Two-Photon-Induced Photoluminescence from Single Gold Nanorods and Imaging of Plasmon Modes[END_REF][START_REF] Imura | Imaging of Surface Plasmon and Ultrafast Dynamics in Gold Nanorods by Near-Field Microscopy[END_REF] .

The luminescence peak analysis was performed using homemade software. It allows performing complex mathematical operations such as non-linear regression (using the ALGLIB library) of peaks on complete mapping or algebraic operations on and between maps and spectra.

Results and Discussion

In Figure 2a we show the luminescence spectra from the apex of the tip when the tip is not in contact with the gold surface and when it is in contact. The intensity of this luminescence is clearly higher at the apex of the tip as can be seen in the mapping shown in Figure 2b andc. A strong luminescence signal has already been reported in previous tip enhanced spectroscopy experiments and was attributed to inelastic light scattering from the adsorbate-substrate system [START_REF] Pettinger | Direct Monitoring of Plasmon Resonances in a Tip-Surface Gap of Varying Width[END_REF] . If it was the case here we should have this signal whatever the tip, and this is not what we observed. We therefore attribute this signal to gold interband transitions enhanced by the tip-sample coupling. In aSNOM the tip can be approximated as a nanorod. The luminescence of gold nanorods has already been studied [START_REF] Mohamed | The "lightning" Gold Nanorods: Fluorescence Enhancement of over a Million Compared to the Gold Metal[END_REF][START_REF] Zijlstra | Single Metal Nanoparticles: Optical Detection, Spectroscopy and Applications[END_REF] : it can be extremely intense compared to a gold surface. The final enhancement and the emission wavelength strongly depend on the shape of the rod; this is surely the reason why luminescence is not always detected even with our tips. In any case, the final intensity of this luminescence is a result of the strong near field enhancement of the incoming and outgoing light via the excitation of the surface plasmon resonance of the tip itself. Several numerical studies have shown that this enhancement can be affected by different parameters, notably the environment of the plasmonic nanostructure. In Figure 2a it appears clearly that luminescence of the tip is enhanced when it is in contact with the gold surface. This is induced by the strong tip-surface coupling [START_REF] Behr | Optical Antenna Properties of Scanning Probe Tips: Plasmonic Light Scattering, Tip-Sample Coupling, and Near-Field Enhancement[END_REF] . Numerical simulations have conclusively shown that when a plasmonic tip is very close to a gold surface a tremendous near field enhancement is observed near the apex [START_REF] Yang | Electromagnetic Field Enhancement in TERS Configurations[END_REF][START_REF] Kazemi-Zanjani | Localized Enhancement of Electric Field in Tip-Enhanced Raman Spectroscopy Using Radially and Linearly Polarized Light[END_REF] . This very enhancement is used in tip-enhanced spectroscopy. In our case, this enhancement increases the intensity of the luminescence. This feature clearly confirms that there is a strong optical coupling between the tip and the gold surface.

b) and c) Mapping of the luminescence intensity of the tip when it is respectively in or not in contact with the gold surface. Insets represent the tip (top) sample (bottom) configuration seen by the scanning objective.

The most commonly used analytical model for the description of the tip sample interaction is the dipole model. The tip is approximated by a sphere (with a radius equal to the curvature radius of the tip). The tip-sample optical coupling is seen as the interaction between this sphere and its image due to the substrate [START_REF] Formanek | Analysis of the Measured Signals in Apertureless near-Field Optical Microscopy[END_REF] . More sophisticated and realistic models have been developed since [START_REF] Esteban | Tip-Substrate Interaction in Optical near-Field Microscopy[END_REF] and the full interpretation of the tip-surface coupling must take into account several parameters: retardation, multipoles excitation, multiple scattering etc. [START_REF] Porto | Resonance Shift Effects in Apertureless Scanning near-Field Optical Microscopy[END_REF] In any case, numerical studies have shown that tipsample coupling sensitively depends on both the tip and sample materials as well as the ambient medium. For linear light scattering the optical image is based on a dielectric contrast [START_REF] Raschke | Apertureless near-Field Optical Microscopy: Tip-Sample Coupling in Elastic Light Scattering[END_REF] . In Figure 3a and b we show the topography and the corresponding elastic light scattered mapping of a gold surface. Even if some similarities are clearly visible, a comparison of both images reveals different information. In the optical image, some details with a typical dimension of a few 10 nm are visible. This confirms the high spatial resolution of optical near field measurements. The final contrast of this image is a complex convolution of topography and optical properties including local and propagating plasmons of the surface. It needs to be said that several studies have shown that the interpretation of a SNOM image is indeed far from obvious and that knowledge of even hardly accessible local materials parameters is required to avoid ambiguities [START_REF] Kohlgraf-Owens | Optical-Force-Induced Artifacts in Scanning Probe Microscopy[END_REF][START_REF] Hecht | Facts and Artifacts in near-Field Optical Microscopy[END_REF] . The most common artifact is induced by the strong dependence of the signal on the tip-sample distance. This explains why the final image is a complex convolution of topography and optical properties of the surface. The deconvolution between those two contributions is beyond our present knowledge of local information. In addition, in most of SNOM studies based on the elastic scattered signal the final contrast is very weak and different experimental techniques must be used to improve it (harmonic demodulation and/or interferometric measurements). In our experimental device, none of these methods is required. But it explains why the contrast in Figure 3b is relatively weak and that some details are not visible. To improve these limitations we have recorded the same optical image using the luminescence signal. The corresponding mapping can be seen in Figure 3c. We would like to point out that, as the luminescence signal is very intense the acquisition time for each point has been limited to 100 ms, the time cost for acquiring a 128 × 128 image is thus less than one hour. It offers an excellent compromise between spatial resolution and limited drift of the piezoelectric controller. The obtained image offers strong similarities with both the topography (Figure 3a) and optical properties (Figure 3b). Some details not visible in Figure 3b can now be clearly distinguished. The very high spatial resolution of the optical image is remarkable. Some spatial details that are hardly visible with the AFM topography image appear clearly on the optical image. This result confirms that the final resolution of a SNOM image can be higher than the curvature radius of the probe [START_REF] Konopsky | Operation of Scanning Plasmon near-Field Microscope with Gold and Silver Tips in Tapping Mode: Demonstration of Subtip Resolution[END_REF] . Cross sections (See electronic supplementary information) reveal that the spatial resolution of the optical image is typically around 20 nm. As evocated above the Figure 3b is dominated by topographic artifacts and surely contains little optical information. The fact that Figure 3c is different suggests that this image is now dominated by the optical contribution. To confirm this hypothesis we show in Figure 4a and b respectively the position of the luminescence maximum and the full width half maximum (FWHM). Both parameters were obtained fitting the curves of each point with a lorentzian with our homemade program. The luminescence of a nanorod, as any plasmonic nanostructure, is strongly influenced by the optical properties of its environment. As a consequence, we derive that a combined evolution of intensity, position and FWHM of luminescence is a sign of an optical interaction between the probe and the optical properties of the sample. This is exactly what we observe here: Figure 3c, Figure 4a andb clearly show strong similarities. Another experimental result demonstrates that those images are dominated by the optical contribution. In Figure 4 the region between y = 0.15 m and 0.25 m is clearly brighter than the rest of the image. It appears that the luminescent wavelength is slightly higher (the typical shift is a few nanometers) and the FWHM has increased. We have observed that for this region some Raman modes are present (see Figure 5). They originate from organic impurities that were temporarily adsorbed at the tip apex and therefore modify the dielectric medium around the tip. No Raman modes are observed outside this region. The Raman modes of these parasitic molecules are also strongly enhanced by the near field antenna as for any TERS measurement. Under ambient conditions and considering that the tip represents an optical tweezer, it is almost impossible to prevent such pollution. Among all the parameters evocated above that can affect the nanorod luminescence only one can explain this result: a change in the optical index in the medium surrounding the nanorod. As the amount of deposited molecules is very low the induced index change is certainly very weak, yet it is clearly visible in Figure 4. This result demonstrates the high optical sensitivity of the luminescence and that the corresponding image is dominated by the optical contribution compared to topographic artifacts. As for any SNOM measurements, this optical signal is the result of two different contributions: the optical properties of the sample itself, that is to say a plasmonic surface, and the optical coupling between the tip and the sample. It is extremely difficult to determine the exact contribution between those two features. Pettinger et al. [START_REF] Pettinger | Direct Monitoring of Plasmon Resonances in a Tip-Surface Gap of Varying Width[END_REF] have studied the evolution of luminescence background in TERS measurements. They reported a blue shift and an intensity decrease when the tip-sample distance increases and attributed it to a decreasing electrodynamic coupling between the excited surface plasmons of the tip and the gold surface. Even if we cannot totally exclude this contribution, we do not consider it to be dominant in our measurements for a simple technical reason: each spectrum is acquired during 0.1 s, a very long time compared to what is required for the feedback loop to insure a constant tipsample distance. As a consequence we assume that the main features observed in both Figure 3 and Figure 4 are a sign of the optical properties of the gold surface. Those properties are dominated by the field distribution around plasmonic nanostructures. Such nanostructures exhibit strong near field enhancement at their edge. This is what we observe in Figure 3c: the luminescence intensity is higher at the edge of the gold nanostructures. As a consequence, it appears that this image is the distribution of the electromagnetic field on the gold plasmonic surface.

As a conclusion we suggest that, even if the presence of a tip can affect the optical properties of the sample [START_REF] Bozhevolnyi | Near-Field Optical Microscopy of Localized Excitations on Rough Surfaces: Influence of a Probe[END_REF] , our luminescence-based near field mapping technique is a simple and very efficient tool for the study of the optical properties of surfaces. Our measurements have evidenced several advantages of this method: a very high signal to noise ratio leading to reasonable acquisition time compatible with high resolution scanning, a very high sensitivity of the luminescence parameters with the optical environment of the tip, a simple experimental device without complex post analysis and a technique that can be used with any kind of sample compatible with AFM measurements. Our study has demonstrated that this approach is extremely suitable for the study of plasmonic surfaces. We suggest that it should be applied for the optical study of nanoantennas [START_REF] Novotny | Nano-Optics: Optical Antennas Tuned to Pitch[END_REF] , an extremely promising and growing field in nano-optics and surface enhanced spectroscopy.

Figure 1 :

 1 Figure 1: Experimental setup of the near field scanning optical microscope. It is based on a Tuning Fork AFM optically coupled to a laser excitation and a spectrometer through a side illumination.

Figure 2 :

 2 Figure 2: a) Luminescence from the apex of the gold tip. b) and c) Mapping of the luminescence intensity of the tip when it is respectively in or not in contact with the gold surface. Insets represent the tip (top) sample (bottom) configuration seen by the scanning objective.

Figure 3 :

 3 Figure 3:1 µm × 1 µm scanning images of the gold surface a) Topography b) Intensity of the elastically scattered light c) Intensity of the luminescence.

Figure 4 :

 4 Figure 4: 1 µm × 1 µm scanning images of the gold surface. The luminescence curve is fitted with a lorentzian and the following parameters are reported a) Position of the maximum b) Full Width Half Maximum map of the fitted curve.

Figure 5 :

 5 Figure 5: Luminescence for the region at y = 0.25 µm. Some unidentified Raman modes (sharp peaks) are clearly present. The offset between the four spectra is for clarification only.

Acknowledgements: A.M. acknowledges the Agence Nationale de la Recherche Scientifique (ANR) for financial support under the CARIOCA project (2010-JCJC-918-01). A.R. is grateful for a CFI leaders opportunity fund for his infrastructure as well as an NSERC discovery grant and an FRQNT team grant.


text:
Francis Mairet
email: fran-cis.mairet@inria.fr

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

Hybrid Control of a Bioreactor with Quantized Measurements

Keywords: Hybrid systems, bioreactor, differential inclusions, quantized output, process control

We consider the problem of global stabilization of an unstable bioreactor model (e.g. for anaerobic digestion), when the measurements are discrete and in finite number ("quantized"), with control of the dilution rate. The measurements define regions in the state space, and they can be perfect or uncertain (i.e. without or with overlaps). We show that, under appropriate assumptions, a quantized control may lead to global stabilization: trajectories have to follow some transitions between the regions, until the final region where they converge toward the reference equilibrium. On the boundary between regions, the solutions are defined as a Filippov differential inclusion. If the assumptions are not fulfilled, sliding modes may appear, and the transition graphs are not deterministic.

I. INTRODUCTION

Classical control methods are often based on the complete knowledge of some outputs y(t) of the system [START_REF] Sontag | Mathematical control theory: deterministic finite dimensional systems[END_REF]. By complete, we mean that any output y i is a real number, possibly measured with some noise δ i . The control is then built with this (noisy) measurement. These tools have been successfully applied in many domains of science and engineering, e.g. in the domains of biosystems and bioreactors [START_REF] Dochain | Automatic control of bioprocesses[END_REF]. However, in these domains, detailed quantitative measurements are often difficult, too expensive, or even impossible. A striking example is the measurements of gene expression by DNA-chips, giving only a Boolean measure equal to on (gene is expressed) or off (not expressed). In the domain of bioprocesses, it frequently happens that only a limited number or level of measurements are available (e.g. low, high, very high ...) because the devices only give a discretized semi-quantitative or qualitative measurement [START_REF] Bernard | Non-linear qualitative signal processing for biological systems: application to the algal growth in bioreactors[END_REF]. For this case of quantized outputs, the problem of control has to be considered in a non-classical way: the control cannot be a function of the full continuous state variables anymore, and most likely will change only when the quantized measurement changes. This framework has been considered by numerous works, having their own specificity: quantized output and control with adjustable zoom [START_REF] Nesic | A unified framework for design and analysis of networked and quantized control systems[END_REF], hybrid systems abstracting continuous ones (cf. [START_REF] Lunze | Handbook of hybrid systems control: theory, tools, applications[END_REF] for many examples).

In this paper, we consider a classical problem in the field of bioprocesses: the stabilization of an unstable bioreactor model, representing for example anaerobic digestion, towards a working set point. Anaerobic digestion is one of the most employed process for (liquid) waste treatment [START_REF] Dochain | Automatic control of bioprocesses[END_REF]. Considering a simplified model with two state variables (substrate and biomass), the system has two stable equilibria (and an unstable one, with a separatrix between the two basins of attractions of the respective stable equilibria), one being the (undesirable) washout of the culture [START_REF] Hess | Design and study of a risk management criterion for an unstable anaerobic wastewater treatment process[END_REF]. The goal is to globally stabilize the process toward the other locally stable reference equilibrium. The (classical) output is the biomass growth (through gaseous production), the control is the dilution rate (see [START_REF] Steyer | Lessons learnt from 15 years of ICA in anaerobic digesters[END_REF] for a review of control strategies). There exists many approaches for scalar continuous output [START_REF] Mailleret | Nonlinear adaptive control for bioreactors with unknown kinetics[END_REF], [START_REF] Antonelli | Set-point regulation of an anaerobic digestion process with bounded output feedback[END_REF], based on well-accepted models [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF]. Here, we suppose that the outputs are discrete or quantized: there are available in the form of finite discrete measurements. Moreover, we introduce an uncertain model where the discrete measurements may overlap, and the true value is at the intersection between two quantized outputs. For this problem, the general approaches described above do not apply, and we have to turn to more tailored methods, often coming from the theory of hybrid systems, or quantized feedbacks. We here develop our adapted hybrid approach. It has also some relations with the fuzzy modeling and control approach: see e.g. in a similar bioreactor process the paper [START_REF] Estaben | Fuzzy control for an anaerobic digester[END_REF]. We provide here a more analytic approach, and prove our results of stability with techniques coming from differential inclusions and hybrid systems theory [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF]. Our work has some relations with qualitative control techniques proposed in [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: the bistable switch example[END_REF], and with the domain approaches used in hybrid systems theory, where there are some (controlled) transitions between regions, forming a transition graph [START_REF] Belta | Controlling a class of nonlinear systems on rectangles[END_REF], [START_REF] Habets | A control problem for affine dynamical systems on a full-dimensional polytope[END_REF].

The paper is organized as follows: Section II describes the bioreactor model and the measurements models. The next section is devoted to the model analysis in open-loop (with a constant dilution rate). In Section IV, we propose a control law and show its global stability through the analysis of the transitions between regions. Finally, in Section V, we explain how to choose the dilution rates by a graphical approach, and we end by giving some simulations.

II. FRAMEWORK A. Model presentation

In a perfectly mixed continuous reactor, the growth of biomass x limited by a substrate s can be described by the following system (see [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF], [START_REF] Dochain | Automatic control of bioprocesses[END_REF]):

ṡ = u(t)(s in -s) -kµ(s)x ẋ = (µ(s) -u(t))x (1)
where s in is the input substrate concentration, u(t) the dilution rate, k the pseudo yield coefficient, and µ(s) the specific growth rate. Given ξ = (s, x), let us rewrite System (1) as ξ = f (ξ, u(t)), where the dilution rate u(t) is the manipulated input. The specific growth rate µ(s) is assumed to be a Haldane function (i.e. with substrate inhibition) [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF]:

µ(s) = μ s k S + s + s 2 /k I (2)
where μ, k S , k I are positive parameters. This function admits a maximum for a substrate concentration s = √ k S k I := s, and we will assume s < s in .

Lemma 1. The solutions of System (1) with initial conditions in the positive orthant are positive and bounded.

Proof. See [START_REF] Mairet | Hybrid Control of a Bioreactor with Quantized Measurements: Extended Version[END_REF].

In the following, we will assume initial conditions within the interior of the positive orthant.

B. Quantized measurements

We consider that a proxy of biomass growth y(ξ) = αµ(s)x is monitored (e.g. through gas production), but in a quantized way, in the form of a more or less qualitative measure: it can be levels (high, medium, low...) or discrete measures. Finally, we only know that y(ξ) is in a given range, or equivalently that ξ is in a given region (parameter α is a positive yield coefficient): n, where 0 = y 1 < y 2 < . . . < y n and y 1 < y 2 < . . . < y n = +∞. We will consider two cases: (A1) Perfect quantized measurements:

Y i = {ξ ∈ R 2 + : y i ≤ y(ξ) ≤ y i }, i = 1, . . . ,
y i = y i+1 , ∀i ∈ {1, . . . , n -1}.
This corresponds to the case where there is no overlap between regions. The boundaries are perfectly defined and measured.

0 y i = y i-1 y i = y i+1 Yi-1 Yi Yi+1 y (A2) Uncertain quantized measurements: y i < y i-1 < y i+1 , ∀i ∈ {2, . . . , n -1}.
In this case, we have overlaps between the regions. In these overlaps, the measure is not deterministic, and may be any of the two values.

0 y i y i-1 y i+1 y i Yi-1 Yi Yi+1 y
For both cases, we define (open) regular domains Ỹi :=

Y i \ (Y i-1 ∪ Y i+1 ), and (closed) switching domains Y i|i+1 := Y i ∩ Y i+1 where the measurement is undetermined, i.e. if ξ ∈ Y i|i+1 , then either ξ ∈ Y i or ξ ∈ Y i+1 .
For perfect measurements (A1), we have Ỹi = intY i , and the switching domains Y i|i+1 correspond to the lines y(ξ) = y i = y i+1 . For uncertain measurements (A2), the switching domains Y i|i+1 become the regions {ξ ∈ R 2 + :

y i+1 ≤ y(ξ) ≤ y i }.
Unless otherwise specified, we consider in the following uncertain measurements.

C. Quantized control

Given the risk of washout, our objective is to design a feedback controller that globally stabilizes System (1) towards a set-point. Given that measurements are quantized, the controller should be defined with respect to each region:

ξ(t) ∈ Y i ⇔ u(t) = D i , i = 1, ..., n. (3)
Here D i is the positive dilution rate in region i. This control scheme leads to discontinuities in the vector fields. Moreover, in the switching domains, the control is undetermined. Thus, solutions of System (1) under Control law (3) are defined in the sense of Filippov, as the solutions of the differential inclusion [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF] ξ ∈ H(ξ), where H(ξ) is defined on regular domains Ỹi as the ordinary function H(ξ) = f (ξ, D i ), and on switching domains Y i|i+1 as the closed convex hull of the two vector fields in the two domains i and i + 1:

H(ξ) = co{f (ξ, D i ), f (ξ, D i+1 )}.
Following [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF], [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: the bistable switch example[END_REF], a solution of System (1) under Control law We first consider System (1) when a constant dilution rate D is applied in the whole space (i.e. u(t) = D, ∀t ≥ 0). In this case, the system is a classical ordinary differential equation which can present bistability, with a risk of washout. Let us denote s a (D) and s b (D) the two solutions, for D ∈ (0, µ(s)), of the equation µ(s) = D, with 0 < s a (D) < s < s b (D). For the Haldane growth rate defined by (2), we have:

(3) on [0, T ] is an absolutely continuous (w.r.t. t) function ξ(t, ξ 0 ) such that ξ(0, ξ 0 ) = ξ 0 and ξ ∈ H(ξ) for almost all t ∈ [0, T ].

III. MODEL ANALYSIS WITH A CONSTANT DILUTION

s a,b (D) = k I 2 μ D -1 ∓ k I 2 μ D -1 2 -k S k I
The asymptotic behavior of the system can be summarized as follows: Proof. See [START_REF] Hess | Design and study of a risk management criterion for an unstable anaerobic wastewater treatment process[END_REF].

For D ∈ (0, µ(s)), let us define y j (D) = αD k [s in -s j (D)], for j = a, b. y a (D) and y b (D) are the growth proxy obtained respectively at the equilibria ξ a (D) and ξ b (D) (if it exists) 1 .

In order to design our control law, we need to provide some further properties of the system dynamics. In particular, we need to characterize ẏ(ξ), the time derivative of y(ξ) along a trajectory of System (1) with a constant dilution rate u(t) = D:

ẏ(ξ) = α [D(s in -s) -kµ(s)x] µ ′ (s)x + αµ(s)(µ(s) -D)x.
Let consider the following functions:

g D : s -→ µ(s) -D kµ ′ (s) + D(s in -s) kµ(s) , h j D : s -→ D(s in -s j (D)) kµ(s) , j = a, b,
defined respectively on (0, s) ∪ (s, +∞) and (0, +∞). In the Proof. See [START_REF] Mairet | Hybrid Control of a Bioreactor with Quantized Measurements: Extended Version[END_REF].

(s, x) plane, g D (s), h a D (s)
1 if D < µ(s in ), ξ b (D) does not exist and y b (D) < 0.
This allows us to determine the monotonicity of y(ξ) in a region of interest (for the design of the control law). Lemma 3. Consider System (1) with a constant dilution rate u(t) = D. For ξ ∈ R2 + such that y b (D) < y(ξ) < y a (D), we have ẏ(ξ) > 0.

Proof. See [START_REF] Mairet | Hybrid Control of a Bioreactor with Quantized Measurements: Extended Version[END_REF].

IV. CONTROL WITH QUANTIZED MEASUREMENTS

A. Control design

Our goal is to globally stabilize the system towards an equilibrium with a high productivity (the productivity is the output αµ(s)x), corresponding to a high dilution rate D n where there is bistability in open loop (case (ii) of Proposition 1). Number n is the number of measurements (see section II-B) and the control is such that D 1 < D 2 < . . . < D n . We consider the following control law, based on the quantized measurements y(ξ), and constant within a given region:

∀t ≥ 0, ξ(t) ∈ Y i ⇔ u(t) = D i , (4)
given that the following conditions are fulfilled:

y b (D i ) < y i i = 1, . . . , n, (5)
y a (D i ) > y i i = 1, . . . , n -1, (6)
y a (D n ) > y n-1 . (7)
Thus, for each region Y i , i < n, these conditions guarantee that the stable operating equilibrium for the dilution rate D i (see Proposition 1) is located in an upper region Y j , j > i, while the saddle point is located in a lower region Y k , k < i. These conditions make the equilibrium ξ a (D n ) globally stable, as we will see below by focusing on the transitions between regions [START_REF] Belta | Controlling a class of nonlinear systems on rectangles[END_REF].

B. Transition between two regions

We first characterize the transition between regions: Lemma 4. For any i ∈ {1, ..., n -1}, consider System (1) under Control law (4) with Conditions (5-6) for i, i + 1 2 . All trajectories with initial conditions in Ỹi ∪ Y i|i+1 enter the regular domain Ỹi+1 .

Proof. First, we consider i = n -1. We will apply LaSalle theorem on a domain

Ω := {ξ ∈ intR 2 + | x ≤ max(x(0) + s(0)/k, s in /k), s ≤ max(s(0), s in ), y i-1 < y(ξ) < y † }
for any y † ∈ Ỹi+1 . Let us consider the function V (ξ) = y a (D n ) -y(ξ). Given Conditions (5-6), we get y b (D i ) < y i < y(ξ) < y i < y a (D i ), ∀ξ ∈ Y i . If a constant dilution rate u(t) = D i is applied on Y i , we can apply Lemma 3 to conclude that ẏ(ξ) > 0, so V (ξ) is decreasing. Similarly, V (ξ) is also decreasing on Y i+1 whenever u(t) = D i+1 . Now under Control law (4), we have shown that V (ξ) is decreasing on the regular domains Ỹi and Ỹi+1 . V (ξ) is a regular C 1 function, and can be differentiated along the differential inclusion. On the switching domains Y i|i+1 , we have:

V (ξ) ∈ co ∂V ∂x f (ξ, D i ), ∂V ∂x f (ξ, D i+1 ) < 0.
Thus, V (ξ) is decreasing on Ω. Recalling that the trajectories are also bounded (Lemma 1), we can apply LaSalle invariance theorem [START_REF] Lasalle | The Stability of Dynamical Systems[END_REF] on the domain Ω. Given that the set of all the points in Ω where V (ξ) = 0 is empty, any trajectory starting in Ω will leave this region. The boundaries x = max(x(0) + s(0)/k, s in /k) and s = max(s(0), s in ) are repulsive (see Proof of Lemma 1). Finally, the boundary y(ξ) = y i corresponds to the maximum of V (ξ) on Ω, so every trajectory will reach the boundary y(ξ) = y † , i.e. it will enter Ỹi+1 . For i = n -1, taking any y † ∈ (y n-1 , y a (D n )), we can show similarly that V (ξ) is decreasing on Ω so every trajectory will enter the region Ỹn .

Following the same proof, we can show that the reverse path is not possible, in particular for the last region: Lemma 5. Consider System (1) under Control law (4) with Conditions (5) for n -1, n, Condition [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF] for n -1, and Condition [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: the bistable switch example[END_REF]. The regular domain Ỹn is positively invariant.

C. Global dynamics

Now we are in a position to present the main result of the paper:

Proposition 2. Control law (4) under Conditions (5-7) with perfect or uncertain measurements (A1 or A2) globally stabilizes System (1) towards the point ξ a (D n ).

Proof. From Lemmas 4 and 5, we can deduce that every trajectory will enter the regular domain Ỹn , and that this domain is positively invariant. System (1) under a constant control u(t) = D n has two non-trivial equilibria (see Proposition 1): ξ a (D n ), and ξ b (D n ). The growth proxy at these two points satisfy y a (D n ) > y n-1 and y b (D n ) < y n (Conditions (5,7)), so there is only one equilibrium in Ỹn : ξ a (D n ). Moreover, it is easy to check that Ỹn is in the basin of attraction of ξ a (D n ), therefore all trajectories will converge toward this equilibrium.

V. IMPLEMENTATION OF THE CONTROL LAW

A. How to fulfill Conditions (5-7)

The global stability of the control law is based on Conditions [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF][START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF][START_REF] Chaves | Exact control of genetic networks in a qualitative framework: the bistable switch example[END_REF]. We now wonder how to easily check if these conditions hold, or how to choose the dilution rates and/or to define the regions in order to fulfill these conditions. In this purpose, a graphical approach can be used. As an example, we will consider the case where the regions are imposed (by technical constraints) and we want to find the different dilution rates D i such that Conditions (5-7) hold. Our objective is to globally stabilize the equilibrium point ξ a (D * ) ∈ Ỹn , with µ(s in ) < D * < µ(s ⋄ ) (s ⋄ is defined just after). Let φ(s) := α k µ(s)(s in -s), which represents the steady state productivity. On [0, s in ], φ(s) admits a maximum for

s ⋄ := s in 1 + 1 + sin kS 1 + sin kI < s.
Note that we impose D * < µ(s ⋄ ) given that for any D * > µ(s ⋄ ), the same productivity can be achieved with a smaller dilution rate, leading to a reduced risk of instability. Let us denote s c (y) and s d (y) the two solutions, for y ∈ (0, φ(s ⋄ )), of the equation φ(s) = y, with 0 < s c (y) < s ⋄ < s d (y) < s in . For the Haldane growth rate, we have:

s c,d (y) = μs in -ky α ∓ ky α -μs in 2 -4k S ky α ( ky αkI + μ) 2 ky αkI + μ .
How to choose D i in one region: For given lower and upper bounds y i < y i < φ(s ⋄ ), we define D(y i ) := min(µ(s ⋄ ), (µ • s d )(y i )) and D(y i ) := (µ • s c )(y i ). This can be done analytically or graphically, as shown on Figure 2A. Whenever we choose D i ∈ (D(y i ), D(y i )), we have 6)). so Conditions (5-6) for i hold. If D(y i ) ≥ D(y i ) (see Figure 2A, dashed lines), it is not possible to fulfill the conditions and thus to implement the control law with this measurement range.

• s ⋄ < s d (y i ) < s b (D i ), so y b (D i ) < y i (Condition (5)), • s c (y i ) < s a (D i ) < s ⋄ , so y a (D i ) > y i (Condition (
How to choose all the D i : The procedure proposed in the previous subsection should be repeated for all the regions. We can depict two particular cases:

• For the first region Y 1 , given y 1 = 0, we actually impose D 1 < µ(s in ). • For the last region Y n , one should check that D(y n-1 ) < D n = D * < D(y n ). This approach is illustrated in [START_REF] Mairet | Hybrid Control of a Bioreactor with Quantized Measurements: Extended Version[END_REF].

Increasing measurement resolution: We have seen that for a given region, it is not always possible to fulfill Conditions [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF][START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF][START_REF] Chaves | Exact control of genetic networks in a qualitative framework: the bistable switch example[END_REF]. This gives rise to a question: is there any constraint on the measurement that guarantees the possibility to implement the control law? For perfect measurements (A1) with equidistribution, we will show that it will always be possible to implement the control law increasing the number of regions.

First, we can arbitrarily define the lower bound of the last region: y n ∈ (y b (D * ), y a (D * )), recalling nonetheless that n is unknown. Now, we will determine the limit range of a measurement region. Let ŷ such that D(ŷ) = µ(s in ). For y i < y i < ŷ, we have D(y i ) < µ(s in ) < D(y i ) so it is always possible to choose a D i in order to fulfill Conditions [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF][START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF].

For y ∈ [ŷ, y n ], we define ψ(y

) := (φ • s b • µ • s c )(y)
, see Figure 2B. It gives a lower bound on y i : whenever one choose y i > ψ(y i ), then we have D(y i ) > D(ψ(y i )) = D(y i ). Thus, one can choose any D i ∈ (D(y i ), D(y i ) and Conditions (5-6) for i will hold.

We can now define the mapping ∆ : y → y -ψ(y) on [ŷ, y n ]. ∆(y) defines the maximal range of the region with upper bound y (in order to be able to find a dilution rate such that Conditions (5-6) hold). ∆ admits a minimum ∆(y) > 0 on [ŷ, y n ]. Thus, for n > y n ∆m + 1, the regions Y i defined by y i = y i+1 = i n-1 y n , i = 1, ..., n -1 allow the implementation of the control law. In conclusions, whenever the measurement resolution is good enough (i.e. the number of regions is high enough), it is always possible to find a set of dilution rates D i such that Conditions (5-7) hold in the case of perfect measurements (A1) with equidistribution. For uncertain measurements (A2), the same result holds, but the proof is omitted for sake of brevity.

B. Simulations

As an example, we consider the anaerobic digestion process, where the methane production rate is measured. More details on parameters, set-point and measurements can be found

Y 2|3 Ỹ4 Ỹ2 Y 1|2 Ỹ1 Y 3|4

Ỹ3

(a) Conditions (5-7) are fulfilled (see Fig. 3) 5)-2 does not hold (see Fig. 5) 6)-2 does not hold (see Fig. 5) Fig. 4. Transitions between regions. Ellipses: regular domains; rectangles: switching domains. White regions are transient, while grey regions have a stable equilibrium. Ellipses with a thick line are positively invariant.

Y 2|3 Ỹ3 Ỹ2 Y 1|2 Ỹ1 (b) Condition (
Y 2|3 Ỹ3 Ỹ2 Y 1|2 Ỹ1 (c) Condition (
in [START_REF] Mairet | Hybrid Control of a Bioreactor with Quantized Measurements: Extended Version[END_REF]. For uncertain measurements, we use discrete time simulation. At each time step t k (with ∆t = 0.05d), when ξ(t k ) is in a switching region Y i|i+1 , we choose randomly the control u(t k ) between D i and D i+1 . In this case, we perform various simulations for a same initial condition. With four regions (n = 4), we can define a set of dilution rates such that Conditions (5-7) are fulfilled. Trajectories for various initial conditions are represented in the phase portrait for uncertain measurements, see Fig. 3. In accordance with Proposition 2, all the trajectories converge towards the set-point. Thus, the transition graph is deterministic (there is only one transition from a region to the upper one), see Fig. 4a. If the number of regions is reduced (three regions only), it is not possible in this example to choose dilution rates such that Conditions (5-7) hold for all i. In this case, some trajectories do not converge towards the set-point. This aspect will be further discussed in the next subsection.

C. When conditions are not verified: risk of failure

We here detail what happens if Conditions (5-7) are not fulfilled, and in particular if there is a risk of washout. First, given the previous analysis of the system, one can easily see that only the condition y b (D 1 ) < y 1 , i.e. D 1 < µ(s in ) is necessary to prevent a washout, so D 1 can be chosen with a safety margin in order to avoid such situation. Now, if Condition (5) is not fulfilled for some i > 1, the unstable equilibrium ξ b (D i ) will be located in the region Y i . Thus, the region Ỹi have transitions towards the lower region, and a trajectory can stay in the switching domain Y i-1|i . Given that z = s + kx converges towards s in , such trajectory will converge towards the intersection between the switching domain and the invariant manifold z = s in (see Fig. 4b and Fig. 5):

• For perfect measurements (A1), this gives rise to a sliding mode and the convergence towards a singular equilibrium point. 5) or ( 6) are not fulfilled . Same legend as Fig. 3. Left: Condition (6)-2 is not fulfilled. Some trajectories converge towards a singular equilibrium point (black diamond) with a sliding mode. Right: Conditions (5)-2 is not fulfilled. Some trajectories (left of the figure) converge towards an equilibrium point in the region Ỹ2 . See also the associated transition graph in Fig. 4b and4c.

• For uncertain measurements (A2), all the trajectories converge towards a line segment. In our simulation, they actually also converge towards a singular equilibrium point.

This situation can be detected by the incessant switches between two regions. In such case, the dilution rate D i should be slightly decreased.

On the other hand, if Condition (6) does not hold for some i, the stable equilibrium ξ a (D i ) will be located in the region Y i , so some trajectories can converge towards this point instead of going to the next region, see Fig. 4c and Fig. 5. To detect such situation, a mean escape time for each region can be estimated using model simulations. A trajectory which stay much more than the escape time in one region should have reached an undesirable equilibrium. In this case, the dilution rate D i should be slightly increased.

In all the cases, the trajectories converge towards a point or a line segment. Although it is not desired, this behavior is particularly safe (given that there is theoretically no risk of washout). Moreover, as explained above, these situations can easily be detected and the dilution rates can be changed accordingly (manually or through a supervision algorithm).

VI. CONCLUSION

Given the quantized measurements, we were able to design (under some conditions) a control based on regions and transition between regions. Moreover, we have seen that for some undesirable cases, singular behaviors (sliding modes) are possible on the boundaries between regions. We think that this kind of control on domains, and the design of the resulting transition graph, is a promising approach, that we want to deepen in future works. This approach could be generalized to other classical systems, e.g. in mathematical ecology.

Fig. 1 .

 1 Fig. 1. Phase portrait of System (1) with a constant dilution rate u(t) = D ∈ (µ(s in ), µ(s)) (case ii of Proposition 1). Magenta lines: trajectories, cyan dashed lines: nullcline ẏ(ξ) = 0, green dash dotted lines: isolines y(ξ) = ya(D) and y(ξ) = y b (D), red dotted line: separatrix, black vertical line: s = s, dark circles: stable equilibria, open circle: unstable equilibrium.

Proposition 1 .

 1 Consider System (1) with a constant dilution rate u(t) = D and initial conditions in the interior of the positive orthant. (i) If D < µ(s in ), the system admits a globally exponentially stable equilibrium ξ a (D) = s a (D), sin-sa(D) k . (ii) If µ(s in ) < D < µ(s), the system admits two locally exponentially stable equilibria, a working point ξ a (D) = s a (D), sin-sa(D) k and the washout ξ 0 = (s in , 0), and a saddle point ξ b (D) = s b (D), sin-s b (D) k , see Figure 1. (iii) If D > µ(s), the washout ξ 0 = (s in , 0) is globally exponentially stable.

1 . 2 .

 12 and h b D (s) represent respectively the nullcline ẏ(ξ) = 0 and the isolines y(ξ) = y a (D) and y(ξ) = y b (D) (i.e. passing through the equilibria ξ a (D) and ξ b (D)), see Figure 1. Knowing that the nullcline ẏ(ξ) = 0 is tangent to the isoline y(ξ) = y a (D) (resp. y(ξ) = y b (D)) at the equilibrium point ξ a (D) (resp. ξ b (D)), we will determine in the next lemma the relative positions of these curves, see Fig. Lemma Consider System (1) with a constant dilution rate u(t) = D. (i) For s ∈ (0, s), we have g D (s) ≥ h a D (s): the nullcline ẏ(ξ) = 0 is above the isoline y(ξ) = y a (D). (ii) For s ∈ (s, s in ), we have g D (s) ≤ h b D (s): the nullcline ẏ(ξ) = 0 is below the isoline y(ξ) = y b (D).

Fig. 2 .

 2 Fig. 2. How to choose D in one region (see Section V-A). A, thin lines: one can choose any D ∈ (D(y), D(y)); dashed lines: it is impossible to fulfill Conditions (5-6), given that D(y) > D(y)); B: ψ(y) represents the lower bound for y (limit case).

Fig. 3 .

 3 Fig. 3. Trajectories (magenta lines) with Control law (4) under uncertain measurements (A2) for various initial conditions in the phase portrait. Conditions (5-7) are fulfilled, so all the trajectories converge towards the set-point (dark circle), see Proposition 2. Open circle: washout. The frontiers are represented by the green dashed lines, switching regions are colored in gray.

Fig. 5 .

 5 Fig. 5. Trajectories with Control law (4) when Conditions (5) or (6) are not fulfilled . Same legend as Fig.3. Left: Condition (6)-2 is not fulfilled. Some trajectories converge towards a singular equilibrium point (black diamond) with a sliding mode. Right: Conditions (5)-2 is not fulfilled. Some trajectories (left of the figure) converge towards an equilibrium point in the region Ỹ2 . See also the associated transition graph in Fig.4b and 4c.

or if i = n -1, Conditions (5) for n -1, n, Condition (6) for n -1, and Condition[START_REF] Chaves | Exact control of genetic networks in a qualitative framework: the bistable switch example[END_REF].

This work was supported by Phycover (ANR-14-CE04-0011), PEPS BMI-Pectolyse, and Investissement d'avenir Reset projects.