François Boulier
email: francois.boulier@lifl.fr

François Lemaire
email: francois.lemaire@lifl.fr

Differential algebra and system modeling in cellular biology

Among all the modeling approaches dedicated to cellular biology, differential algebra is particularly related to the well-established one based on nonlinear differential equations. In this paper, it is shown that differential algebra makes both simple and algorithmic one of the model reduction methods, the quasi-steady state approximation theory, in the particular setting of generalized chemical reactions systems. This recent breakthrough may suggest some evolution of modeling techniques based on nonlinear differential equations, by incorporating the reduction hypotheses in the models. Potential improvements of parameters fitting methods are discussed too.

Introduction

Among all the modeling approaches dedicated to cellular biology, differential algebra is particularly related to the well-established one based on nonlinear differential equations [START_REF] Leloup | Modeling the molecular regulatory mechanism of circadian rhythms in Drosophila[END_REF][2][START_REF] Conrad | Modeling Molecular Interaction Networks with Nonlinear Differential Equations[END_REF].

Differential equations apply to rather small models. There are many subapproaches, restricted to (sometimes piecewise) linear differential equations, qualitative simulations . . . [START_REF] De | Modeling and Simulation of Genetic Regulatory Systems: A Literature Review[END_REF][START_REF] De Jong | Genetic Network Analyzer: qualitative simulation of genetic regulatory networks[END_REF][START_REF] De Jong | Qualitative Approaches to the Analysis of Genetic Regulatory Networks[END_REF]. The approach based on nonlinear differential equations was however successfully applied for the analysis of many genetic regulatory circuits, including those that work during complex animal embryogenesis [START_REF] Saka | A mechanism for the sharp transition of morphogen gradient interpretation in Xenopus[END_REF]. One of the major successes of this approach is certainly the modeling of the segment polarity network of the drosophila [START_REF] Von Dassow | The segment polarity network is a robust developmental module[END_REF][START_REF] Von Dassow | Exploring modularity with dynamical models of gene networks[END_REF]. This paper is mostly concerned by the model reduction problem for chemical reactions systems based on the generalized mass action law [START_REF] Horn | General mass action kinetics[END_REF]. Among all the approximation techniques available for such systems [11, lumping, sensitivity analysis, . . . ], this paper is exclusively concerned by the quasi-steady state approximation theory.

The model reduction problem for systems of differential equations broadly consists in simplifying the given system, by means of some simplification hypotheses. It aims to get a tractable reduced system, while preserving some properties of interest (e.g. presence of oscillations, number of equilibria).

Modeling approaches based on nonlinear differential equations do not all directly rely on the formalism of chemical reactions systems. Quite often, more sophisticated functions are used such as the Henri-Michaelis-Menten factors and the Hill functions. It is worth noticing that these sophisticated functions can quite often be deduced from generalized chemical reactions systems through some model reduction. However, practitioners using these functions do not always explicitly formulate the simplifying assumptions which justify these approximations. Thereby, they do not formulate the domains of validity of their models as precisely as they could.

Differential algebra makes both simple and algorithmic the quasi-steady state approximation methods in the particular setting of systems of generalized chemical reactions systems [START_REF] Boulier | Model Reduction of Chemical Reaction Systems using Elimination[END_REF]. This fact thus suggests to widen the use of these systems for systems modeling, and to incorporate in the models the hypotheses which lead to the more sophisticated formulas. One should thereby expect to obtain models with more ascertained domains of validity and to estimate the practical relevance of some of the simplification hypotheses.

In addition to the above issues, the differential algebra and the quasi-steady state approximation theory also suggest improvements of the parameters fitting problem, which is an important concern in systems modeling. It has been known for a few years that differential elimination somehow transforms nonlinear least squares into linear ones in the parameters estimation problem, by suggesting a starting point for the Newton-like estimation methods [START_REF] Denis-Vidal | System identifiability (symbolic computation) and parameter estimation (numerical computation)[END_REF][START_REF] Boulier | proceedings of the ICPSS conference[END_REF]. By incorporating moreover, in these methods, the recent progress obtained in the quasi-steady state approximation theory, one may also expect to reduce the stiffness of the differential systems that they need to numerically integrate, thereby speeding up the overall optimisation processes.

Differential elimination

The content of this section owes a lot to [START_REF] Boulier | Differential Elimination and Biological Modelling[END_REF]. We refer to that paper for more details.

The differential elimination theory is a subtheory of the differential algebra [START_REF] Ritt | Differential Algebra[END_REF][START_REF] Kolchin | Differential Algebra and Algebraic Groups[END_REF]. See also [START_REF] Wang | Elimination Practice: Software Tools and Applications[END_REF]. The differential elimination processes that are presented in this paper take as input two parameters: a system of polynomial (thus nonlinear) differential equations, ordinary or with partial derivatives1 and a ranking2 . They produce on the output an equivalent finite set of polynomial differential systems, which are simpler, in the sense that they involve some differential equations which are consequences of the input system but were somehow hidden. The output may consist of more than one differential system because the differential elimination process may need to split cases. The set of the differential equations which are consequences of the input system forms a so-called differential ideal of some polynomial differential ring. Since this ideal is an infinite set, a natural question arises: how does the process select the finitely many differential equations which appear in the output system ? This is indeed the role of the rankings.

An example

The following example, borrowed from [19, Chapter VII, page 454], is a socalled differential algebraic equations system, since it mixes differential and non differential equations. There are three unknown time varying functions (three dependent variables) x, y and z. The dot over a variable denotes its derivative w.r.t. the independent variable, which is assumed to be the time, throughout this paper : ẋ = 0.7 y + sin(2.5 z), ẏ = 1.4 x + cos(2.5 z),

1 = x 2 + y 2 .
Let us try to apply an explicit numerical integration scheme such as Euler's one over this system, for some initial conditions. For each dependent variable (say) x, the scheme transforms the explicit differential equation ẋ = f (x, y, z) as a recurrence formula x n+1 = x n +h f (x n , y n , z n ), where h denotes some stepsize, and computes the terms x n the ones after the others. However, the absence of a differential equation of the following form seems to forbid the application of Euler's scheme:

ż = f (x, y, z). (1)
Indeed, this equation is not missing. It is a "hidden" consequence of the input system: it can be automatically extracted from it by differential elimination. In order to show how to proceed with the help of the diffalg package of MAPLE, one first converts the system as a polynomial differential system. For this, one denotes s the sine, c the cosine and one introduces a few more equations. The following differential polynomial system is equivalent to the above one.

ẋ = 0.7 y + s, ẏ = 1.4 x + c, 1 = x 2 + y 2 . ṡ = 2.5 ż c, ċ = -2.5 ż s, 1 = s 2 + c 2 .
In order to compute equation (1) using diffalg, one stores the differential polynomial system in the variable syst, converting floating point numbers as rational numbers. Then the variable R is assigned the context of the computation: the only derivation is taken with respect to the time, the notation is the standard diff notation of MAPLE and the ranking is provided. According to the ranking notation of the diffalg package, this is the orderly ranking such that s > c > x > y > z.

The properties of that ranking ensure that, if there exists a differential equation of the form (1) in the radical of the differential ideal generated by syst then the differential elimination process will find it. Next the Rosenfeld-Gröbner function [START_REF] Boulier | Representation for the radical of a finitely generated differential ideal[END_REF] is applied to syst and R. It returns a list of MAPLE tables. Each table provides a regular differential chain defining some differential ideal. The list should be understood as an intersection. Over the example, the list only involves one entry so that the regular differential chain does represent the radical differential ideal generated by the input system. The desired equation stands on the second place of the chain (only the two first equations are displayed). Enlarging the input system with this equation, it is now easy to perform any numerical integration method and our problem is solved.

Differential algebra

A differential ring (resp. field) is a ring (resp. field) R endowed with a derivation (this paper is restricted to the case of a single derivation but the theory is more general) i.e. a unitary mapping R → R such that (denoting ȧ the derivative of a):

˙ (a + b) = ȧ + ḃ, ˙ (a b) = ȧ b + a ḃ.
Observe that, theoretically, the derivation is an abstract operation. For legibility, one views it as the derivation w.r.t. the time t. Algorithmically, one is led to manipulate finite subsets of some differential polynomial ring R = K{U } where K is the differential field of coefficients (in practice, K = Q, Q(t) or Q(k 1 , . . . , k r ) where the k i denote parameters that would be assumed to be algebraically independent) and U is a finite set of dependent variables 3 . The elements of R, the differential polynomials are just polynomials in the usual sense, built over the infinite set, denoted ΘU , of all the derivatives of the dependent variables. Definition 1. A differential ideal of a differential ring R is an ideal of R, stable under the action of the derivation.

Let F be a finite subset of a differential ring R. The set of all the finite linear combinations of various orders derivatives of elements of F , with elements of R for coefficients, is a differential ideal. It is called the differential ideal generated by F . An ideal A is said to be radical if a ∈ A whenever there exists some nonnegative integer p such that a p ∈ A. The radical of an ideal A is the set of all the ring elements a power of which belongs to A. The radical of a (differential) ideal is a radical (differential) ideal.

Theorem 1. Let R be a differential polynomial ring and F be a finite subset of R. A differential polynomial p of R lies in the radical of the differential ideal generated by F if and only if it vanishes over every analytic solution of F .

Proof. [START_REF] Ritt | Differential Algebra[END_REF]chap. II,[START_REF] Saka | A mechanism for the sharp transition of morphogen gradient interpretation in Xenopus[END_REF][START_REF] Okino | Simplification of Mathematical Models of Chemical Reaction Systems[END_REF] or [START_REF] Boulier | A computer scientist point of view on Hilbert's differential theorem of zeros[END_REF].

The Rosenfeld-Gröbner algorithm [START_REF] Boulier | Representation for the radical of a finitely generated differential ideal[END_REF] solves the membership problem to radical differential ideals. To present it, one needs to define the concept of ranking. Definition 2. If U is a finite set of dependent variables, a ranking over U is a total ordering over the set ΘU of all the derivatives of the elements of U which satisfies: a < ȧ and a < b ⇒ ȧ < ḃ for all a, b ∈ ΘU .

Let U be a finite set of dependent variables. A ranking such that, for every u, v ∈ U , the ith derivative of u is greater than the jth derivative of v whenever i > j is said to be orderly. If U and V are two finite sets of dependent variables, one denotes U ≫ V every ranking such that any derivative of any element of U is greater than any derivative of any element of V . Such rankings are said to eliminate U w.r.t. V .

Assume that some ranking is fixed. Then one may associate with any differential polynomial f ∈ K{U } \ K the greatest (w.r.t. the given ranking) derivative v ∈ ΘU such that deg(f, v) > 0. This derivative is called the leading derivative or the leader of f . Rankings permit to define leaders. Leaders permit to use differential polynomial as rewrite (substitution) rules. Assume that f = a d v d + • • • + a 1 v + a 0 is a differential polynomial with leader v (the coefficients a i are themselves differential polynomials). Then the equation f = 0 can be written (as the rewrite rules function of diffalg presented in section 2.1 does, though it uses an equal sign instead of an arrow):

v d -→ - a d-1 v d-1 + • • • + a 1 v + a 0 a d •
It can be used afterwards as a rule to simplify any differential polynomial g such that deg(g, v) ≥ d or deg(g, v (k) ) > 0 where v (k) denotes any proper derivative of v. There are precise algorithms for performing these sorts of substitution by finite sets of rewrite rules: Ritt's reduction algorithm or the normal form algorithm [START_REF] Boulier | A computer scientist point of view on Hilbert's differential theorem of zeros[END_REF]algorithm NF].

The Rosenfeld-Gröbner algorithm gathers as input a finite system F of differential polynomials and a ranking. It returns a finite family (possibly empty) C 1 , . . . , C r of finite subsets of K{U } \ K, called regular differential chains. Each system C i defines a differential ideal C i (it is a characteristic set of C i ) in the sense that, for any f ∈ K{U }, we have

f ∈ C i iff NF(f, C i ) = 0.
The relationship with the radical A of the differential ideal generated by F is the following:

A = C 1 ∩ • • • ∩ C r .
When r = 0 we have A = K{U }. Combining both relations, one gets an algorithm to decide membership in A. Indeed, given any f ∈ K{U } we have:

f ∈ A iff NF(f, C i ) = 0, 1 ≤ i ≤ r.
The differential ideals C i do not need to be prime. They are however necessarily radical. The NF( • , C i ) function permits to compute canonical representatives of the residue classes of the differential ring R/C i .

3 Quasi-steady state approximation for generalized chemical reactions systems Differential elimination makes the quasi-steady state approximation (QSSA) theory both simple and algorithmic in the special setting of generalized4 chemical reactions systems as shown by [START_REF] Boulier | Model Reduction of Chemical Reaction Systems using Elimination[END_REF].

QSSA in general: fast and slow variables

In principle, QSSA applies to systems under the two time-scales standard form. Consider the following system in two dependent variables x and y (for legibility, but there may be more than two variables) and assume that ε is a small positive number.

ẋ = f (x, y), ε ẏ = g(x, y). (2)
On a random point (x, y) ∈ R 2 , the speed of y is high and thus, under some general conditions, rapidly approaches an area where g(x, y) ≃ 0. The variable x is said to be slow while the variable y is said to be fast. The QSSA amounts to approximate such a system by the following one:

ẋ = f (x, y), 0 = g(x, y) (3)
which mixes differential and algebraic equations. Note that this approximation is only valid under some conditions (e.g. stability of the differential system in the neighborhhod of g(x, y) = 0) given in the Tikhonov theorem [START_REF] Kokotovic | Singular Perturbation Methods in Control: Analysis and Design[END_REF] and after the transient step is elapsed. Performing a QSSA over a differential system presents two advantages: it reduces the number of ODE and it tends to transform stiff systems into nonstiff ones, much easier to solve numerically.

However, the QSSA is not proven algorithmic in general. The issue relies in the fact that the fast and slow variables, if they exist, may only be obtained through a change of coordinates and there does not seem to exist any algorithm which decides, given a differential system, if such a change of coordinates does exist.

QSSA for chemical reactions systems: fast and slow reactions

For differential systems arising from generalized chemical reactions systems, there exists a standard way to perform the QSSA, provided that the set of chemical reactions is divided in two parts: the fast ones and the slow ones.

As far as we know, the first clear relationship between this method and the Tikhonov theorem was established in [START_REF] Van Breusegem | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF]. Afterwards, close variants of the same method were rediscovered more or less independently [START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF][START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF]. Though all these papers present methods, none of them is fully presented in an algorithmic manner. This may at least partly be due to the fact that some steps of the methods require the inversion of a matrix over a residue class ring, a non obvious task which may imply splitting cases.

Indeed, it turns out that the whole method is equivalent to a differential elimination process, as shown for the first time in [START_REF] Boulier | Model Reduction of Chemical Reaction Systems using Elimination[END_REF]. For a general presentation of the method, one refers to [START_REF] Boulier | Model Reduction of Chemical Reaction Systems using Elimination[END_REF]. In this paper, one only presents the method over a famous example: the Henri-Michaelis-Menten reduction of the basic enzymatic reaction system:

E + S k1 ----→ ← ---- k2 C k3 ----→ E + P (4)
The initial system of ODE writes:

Ẋ = N V i.e.     Ė Ċ Ṡ Ṗ     =     -1 1 1 1 -1 -1 -1 1 0 0 0 1     •   k 1 E S k 2 C k 3 C   ( 5
)
where X is the vector of the chemical species, N is the system stoichiometry matrix and V is the vector of the reaction rates. The stoichiometry matrix is built as follows: it involves one row per species and one column per reaction. The entry at row r, column c is the number of molecules of species r produced by the reaction c (i.e. the number of times species r occurs on the reaction right handside minus the number of times it occurs on the reaction left handside). The rate of a reaction is the product of the left handside species (with multiplicities) times the reaction rate (the parameter over the arrow). Expanding the formula, one gets:

Ė = -k 1 E S + (k 2 + k 3 ) C, Ṡ = -k 1 E S + k 2 C, Ċ = k 1 E S -(k 2 + k 3 ) C, Ṗ = k 3 C. (6)
Among all the assumptions leading to the Henri-Michaelis-Menten formula:

Ṡ = - V max S K + S , (7)
the only one concerning the QSSA is that k 1 , k 2 ≫ k 3 i.e. that the two leftmost reactions of the chemical system (4) are fast while the rightmost one is slow.

The other assumptions will be given later.

In order to perform the QSSA over the above system, one builds a DAE system from system (6) by replacing the contribution of the fast reactions by a new dependent variable, F 1 , and by inserting an algebraic equation, defining the "slow variety", stating that the fast reactions are at quasi-equilibrium.

Ė = -F 1 + k 3 C Ṡ = -F 1 Ċ = F 1 -k 3 C Ṗ = k 3 C 0 = k 1 E S -k 2 C. (8)
In general there may be many different new dependent variables F i and many different algebraic equations. Observe that, over the example, the two fast reactions are considered as one reversible reaction and associated to only one dependent variable: F 1 . The general process is described in section 3.3. Differential elimination is performed below using diffalg. The parameters are put in the base field F of the equations to avoid discussions over their values (they are then considered as algebraically independent). The inequation C(t) = 0 is inserted to avoid considering this degenerate case. The output is pretty printed.

F 1 = k 3 k 1 E S (k 1 S + k 2 ) k 2 (k 1 S + k 1 E + k 2 ) , Ė = k 2 1 E 2 k 3 S k 2 (k 1 S + k 1 E + k 2 ) , Ṗ = k 3 k 1 E S k 2 , Ṡ = - k 3 k 1 E S (k 1 S + k 2 ) k 2 (k 1 S + k 1 E + k 2 ) , C = k 1 E S k 2
Of course, one does not recognize formula [START_REF] Saka | A mechanism for the sharp transition of morphogen gradient interpretation in Xenopus[END_REF] in this regular differential chain: the reduction is incomplete since the extra assumptions made in the Henri-Michaelis-Menten reduction have not yet been taken into account. These assumptions are: S(0) ≫ E(0), P ≃ 0 and C(0) = 0. One refers to [START_REF] Boulier | Model Reduction of Chemical Reaction Systems using Elimination[END_REF] for a complete reduction.

Construction of the DAE to consider for the QSSA

First split the stoichiometry matrix N into two matrices N f and N s putting the columns which correspond to fast reactions in N f and the ones which correspond to slow reactions in N s . Split accordingly the rows of the vector V into two vectors V f and V s . One gets a formula Ẋ = N s V s + N f V f . Over system (5), one gets:

    Ė Ċ Ṡ Ṗ     =     1 -1 0 1     • k 3 C +     -1 1 1 -1 -1 1 0 0     • k 1 E S k 2 C . ( 9
)
Determine a maximal linearly independent set of columns of N f (i.e. a basis of that matrix) and remove the other ones, giving a new matrix N f . Update the vector of reaction rates V f , giving a new vector

V f such that N f V f = N f V f .
Over the example, removing the second column, one gets a new formula Ẋ = N s V s + N f V f giving system [START_REF] De Jong | Qualitative Approaches to the Analysis of Genetic Regulatory Networks[END_REF]:

    Ė Ċ Ṡ Ṗ     =     1 -1 0 1     • k 3 C +     -1 1 -1 0     • k 1 E S -k 2 C . (10)
Replace the vector V f by a vector F of new dependent variables F i . The slow variety is defined by letting the entries of V f all equal to zero. The DAE to be considered for quasi-steady state approximation is

Ẋ = N s V s + N f F, V f = 0.
Over the example, one gets a formula giving system (8):

    Ė Ċ Ṡ Ṗ     =     1 -1 0 1     • k 3 C +     -1 1 -1 0     • F 1 , k 1 E S -k 2 C = 0. ( 11
)

Limits and generalizations of the method

In some cases, the F i variables cannot be all eliminated. In that case, the method does not apply. In [START_REF] Boulier | Model Reduction of Chemical Reaction Systems using Elimination[END_REF], this is checked by testing whether NF( Ẋi , C) only depends on the X j (with order 0), for each dependent variable X i and each regular differential chain C produced by differential elimination. A simple example is given by the generalized chemical reactions system, assuming k 1 , k 2 are both fast:

A k1 ----→ B, A k2 ----→ C.
If one applies the method sketched above, one gets the following differentialalgebraic system:

Ȧ = -F 1 -F 2 , Ḃ = F 1 , Ċ = F 2 , 0 = k 1 A, 0 = k 2 A
which simplifies to the regular differential chain:

F 1 = -Ċ, F 2 = Ċ, Ḃ = -Ċ, A = 0.
The normal form of Ḃ is -Ċ and Ċ is equal to its own normal form: the method failed.

According to the Tikhonov theorem, there are some extra conditions to check for the approximation to be valid [23, conditions C3 and C4]. In particular, the slow variety defined by the algebraic equations must be attractive.

In some cases, there exists a better slow variety than the one provided by the method. To compute a reduced system w.r.t. this different variety, one just needs to change the set of algebraic equations of the differential-algebraic system and to run the differential elimination process. Such a situation is illustrated by the next example borrowed from [26] where

k 2 ≫ k 1 : A k1 ----→ B k2 ----→ ∅
If one applies the method, one gets the following differential-algebraic system:

Ȧ = -k 1 A, Ḃ = k 1 A -F 1 , 0 = k 2 B
which simplifies to the following regular differential chain:

F 1 = k 1 A, Ȧ = -k 1 A, B = 0.
However, the slow variety k 1 A-k 2 B = 0 is better than B = 0. A better reduced system is thus obtained by performing differential elimination over the following differential-algebraic system:

Ȧ = -k 1 A, Ḃ = k 1 A -F 1 , 0 = k 1 A -k 2 B
which simplifies to the following regular differential chain:

F 1 = (k 1 + k 2 ) B, Ḃ = -k 1 B, A = k 2 k 1 B.
This section aims at establishing some relationship between the quasi-steady state approximation method sketched in section 3 and system modeling in cellular biology. For this purpose, an application of this technique to the analysis of a genetic regulatory circuit made of a single gene, regulated by a polymer of its own protein [START_REF] Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF][START_REF] Boulier | Applying a rigorous quasisteady state approximation method for proving the absence of oscillations in models of genetic circuits[END_REF] is described. One considers the genetic circuit depicted in Figure 1. The single gene is regulated by an order n polymer of its own protein. The integer number n is a parameter of the system. This study was motivated by the activity of a working group aiming at modeling the circadian clock of the green alga ostreococcus tauri. The addressed question was: does there exist biologically meaningful parameters values which make this circuit oscillate ? More technically: does there exist biologically meaningful parameters values which make a Poincaré-Andronov-Hopf bifurcation occur ? One refers to [START_REF] Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF] for a more detailed motivation of the addressed question and to [START_REF] Niu | Algebraic Approaches to Stability Analysis of Biological Systems[END_REF] for a related work.

There are many different ways to derive a system of ordinary differential equations from the considered circuit but one of the simplest schemes consists in first translating it as a system of generalized chemical reactions (observe that transcription and translation are not balanced reactions). The variables G and H represent the state of the gene. The mRNA concentration and the concentration of the protein translated from the mRNA are represented by M and P . The n types of polymers of P are denoted by P = P 1 , P 2 , . . . , P n . Greek letters and

k - i , k + i (1 ≤ i ≤ n -1) represent parameters: G + P n α ----→ ← ---- θ H, G ρ f ----→ G + M, H ρ b ----→ H + M, M β ----→ M + P, M δM ----→ ∅, P δP ----→ ∅, P i + P k + i ----→ ← ---- k - i P i+1 (1 ≤ i ≤ n -1).
This generalized chemical reactions system can now be canonically translated as a system of parametric ordinary differential equations, denoting

A i = (k - i P i+1 - k + i P i P ).
Variables G, H, M, P = P 1 , . . . , P n are dependent variables. They all represent species concentrations except G and H, which should rather be viewed as "random variables".

Ġ = θ H -α G P n , Ḣ = -θ H + α G P n , Ṁ = ρ f G + ρ b H -δ M M, Ṗ = β M -δ P P + 2 A 1 + A 2 + • • • + A n-1 , Ṗi = -A i-1 + A i (2 ≤ i ≤ n -1), Ṗn = -A n-1 + θ H -α G P n (12)
This system of n + 3 differential equations depending on 2 n + 5 parameters is actually much too large for any further symbolic analysis. It needs to be reduced.

In order to apply a quasi-steady state approximation, it is assumed that the n -1 chemical reactions describing the polymerization of the protein are fast compared to the other ones.

Then, according to the technique sketched in section 3, one gets an approximation of system ( 12) by replacing each expression A i by a new dependent variable F i (1 ≤ i ≤ n -1) and by augmenting this system by the n -1 following algebraic equations:

0 = k + i P P i -k - i P i+1 , ( 1
≤ i ≤ n -1)
It is now sufficient to eliminate the F i from the so obtained differential-algebraic system. Unfortunately, this cannot be performed by a standard differential elimination algorithm since the number of equations depends on the parameter n. However, such an algorithm can be applied for many different values of n and the general formula can be inferred.

Ġ = θ H -α K n-1 P n G, Ḣ = -θ H + α K n-1 P n G, Ṁ = ρ b H + ρ f G -δ M M, Ṗ = n θ H -n α K n-1 P n G -δ P P + β M n-1 i=0 (i + 1) 2 K i P i (13)
where

K i = k + 1 • • • k + i k - 1 • • • k - i with the convention K 0 = 1.
The redundant equation describing the evolution of H can be removed and H can be replaced by γ 0 -G in the three remaining equations, where γ 0 is some new parameter denoting the total quantity of gene. Some further exact reduction of the parameters set can moreover be performed and one is led to the three parametric ODE model [28, system (3)].

This three parametric ODE reduced system is simple enough to be tackled by the Hopf criterion (a variant of the Routh-Hurwitz one) and it was proven in [START_REF] Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF][START_REF] Boulier | Applying a rigorous quasisteady state approximation method for proving the absence of oscillations in models of genetic circuits[END_REF] that a Poincaré-Andronov-Hopf bifurcation occurs for biologically meaningful values of the parameters if and only if n > 8. One refers to those papers for a biological interpretation of this result.

Let us conclude this section by a few comments on the involved quasi-steady state approximation methods. When [START_REF] Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF] was written, its authors did not know the particular way to perform the quasi-steady state approximation over an ODE system derived from a generalized chemical reactions system. Thus they applied a "handmade" technique inspired from the singular perturbation theory yielding a three ODE model and "solved" it. Integral curves obtained by numerical integration showed differences between the initial and the reduced system but this was expected. Later, the same authors understood the particular technique for generalized chemical reactions system and showed how it is related to differential elimination [START_REF] Boulier | Model Reduction of Chemical Reaction Systems using Elimination[END_REF]. This technique, applied over the very same system in [START_REF] Boulier | Applying a rigorous quasisteady state approximation method for proving the absence of oscillations in models of genetic circuits[END_REF], produced (after a suitable change of coordinates) a slightly different reduced system: one gets the old [START_REF] Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF] reduced system by clearing the denominator of the last ODE of system [START_REF] Denis-Vidal | System identifiability (symbolic computation) and parameter estimation (numerical computation)[END_REF]. By some reduction argument, the proof of the old paper could be applied to the new model. This somehow confirmed that the handmade reduction technique was good. However, comparisons of integral curves obtained by numerical integration showed that the new model is more accurate, in the sense that it has a wider domain of validity, than the old one.

Parameters estimation

The problem which is addressed in this section can be stated as follows: given a parametric ordinary differential system and files of measures for some of the dependent variables, estimate the values of the unknown parameters. In this context, differential elimination somehow transforms a nonlinear least squares problem into a linear one by guessing a starting point for a Newton-like method.

The approach described in this section was published in [START_REF] Noiret | Utilisation du calcul formel pour l'identifiabilité de modèles paramétriques et nouveaux algorithmes en estimation de paramètres[END_REF][START_REF] Denis-Vidal | System identifiability (symbolic computation) and parameter estimation (numerical computation)[END_REF]. It was implemented in a software based on the BLAD libraries [START_REF] Boulier | The BLAD libraries[END_REF][START_REF] Boulier | proceedings of the ICPSS conference[END_REF]. It is strongly related to the study of the identifiability of parametric differential systems, for which a huge literature is available [START_REF] Walter | Identifiability of State Space Models[END_REF][START_REF] Ollivier | Le problème de l'identifiabilité structurelle globale : approche théorique, méthodes effectives et bornes de complexité[END_REF][START_REF] Diop | Nonlinear observability, identifiability, and persistent trajectories[END_REF][START_REF] Ljung | On global identifiability for arbitrary model parametrisations[END_REF][START_REF] Sedoglavic | A Probabilistic Algorithm to Test Local Algebraic Observability in Polynomial Time[END_REF].

To be honest, the method is not very convenient for system modeling in cellular biology for it requires pretty precise model equations as well as pretty accurate measures. In cellular biology, one usually does not know the model equations though the situation might be better for in vitro systems or synthetic genes. But accurate measures are definitely not available in 2008. However, we chose to present the method anyway in this paper for two reasons: first, the use of differential elimination always improves the classical numerical one ; second, the quasi-steady state approximation technique presented in section 3 suggests another, new, improvement.

The content of this section owes a lot to [START_REF] Boulier | Differential Elimination and Biological Modelling[END_REF] and one refers to this paper for more detailed explanations.

The problem is stated over an academic example. Consider the following system of parametric ordinary differential equations. There are two dependent variables x 1 , x 2 and four parameters k 12 , k 21 , k e , V e . It could be derived from a chemical reactions system since it involves two linear exchanges and a Henri-Michaelis-Menten degradation term.

ẋ1 = -k 12 x 1 + k 21 x 2 - V e x 1 k e + x 1 , ẋ2 = k 12 x 1 -k 21 x 2 . ( 14
)
Assume that x 1 is observed i.e. that a file of measures is available for this dependent variable but that x 2 is not. For the sake of simplicity, assume moreover that x 2 (0) = 0 and k e = 7 are known. Issue: estimate the values of the three unknown parameters k 12 , k 21 , V e . There exists a purely numerical method to solve this problem. It is based on a nonlinear least squares method i.e. a Newton method. The idea is simple: pick random values for the three unknown parameters. Integrate numerically the differential system w.r.t. these values and compare the curve obtained by simulation with the file of measures. The error is defined as the sum, for all abscissas, of the squares of the ordinates differences between the two curves. The Newton method (a Levenberg-Marquardt scheme was applied in [START_REF] Boulier | proceedings of the ICPSS conference[END_REF]) updates the values of the three unknown parameters if the error is considered as too large. It stops either if the error is small enough of if a stationary point is reached. However, as every Newton method, nonlinear least squares are very sensitive to the choice of the starting point (the initial random values) and are very likely to end up in a local minimum, with wrong parameters values.

By means of differential elimination, it is sometimes possible to compute a first estimate of the unknown parameters. This estimate is usually not very precise but can be used as a starting point for the Newton method.

The idea consists in eliminating the non observed variables in order to get a relation which only involves the observed variable x 1 , its derivatives up to any order and the parameters. Let us proceed with the help of diffalg. The Rosenfeld-Gröbner algorithm 5 is applied over the model equations. The ranking x 2 ≫ x 1 eliminates x 2 w.r.t. x 1 . The right-hand side of the first model equation is a rational fraction. It is decomposed as a numerator and a denominator. The numerator is stored in the list of the equations (first parameter to Rosenfeld-Gröbner). The denominator is stored in the list of the inequations. To avoid splitting cases on parameters values, one views them as algebraically independent elements of the base field of the differential polynomials. The characteristic set ideal involves two polynomials. The one which does not involve x 2 is the second one, which is displayed below, slightly pretty printed. The expressions enclosed between square brackets are called "parameters blocks".

ẍ1 (x 1 + k e ) 2 + [k 12 + k 21 ] ẋ1 (x 1 + k e ) 2 + [V e ] ẋ1 k e + [k 21 V e ] x 1 (x 1 + k e ) = 0.
This equation tells us that the systems parameters are in principle uniquely defined. Indeed, assume that the function x 1 is known. Then so are its derivatives ẋ1 and ẍ1 . These three functions can therefore be evaluated for three different values of the time t. The known parameter k e can be replaced by its value. One thereby gets an exactly determined system of three linear equations whose unknowns are the parameters blocks. This system admits a unique solution. The values of the parameters blocks being fixed, it is obvious (over this example !) that the values of k 12 , k 21 and V e also are uniquely defined. QED.

In practice, the function x 1 is known from a file of measures and one can try to numerically estimate the values of its first and its second derivative. If the measures are free of noise, the first derivative can be quite accurately estimated but this is usually not the case for the second derivative. To overcome these difficulties due to numerical approximations, one builds an overdetermined linear system that one solves by means of linear least squares.

The values of the blocks of parameters being known, one still has to recover the values of the parameters by solving the above algebraic system. Over this example, it is very easy. The obtained values can then be used as a starting point for the Newton method. Observe that one cannot guarantee that this starting point does actually lead to the global minimum.

There are however two more important difficulties to overcome. There may exist algebraic relations between the parameters blocks. There is no such relation over the example. But assume, for the sake of the explanation, that the computed differential polynomial involves the three following blocks of parameters so that the third block is the product of the two first ones:

[V e ], [k 21 ], [V e k 21 ].
There is no doubt that the numerical values produced during the resolution of the linear overdetermined system would not satisfy this relation. This would imply that the final algebraic system to solve in order to get the values of the parameters would be inconsistent.

The differential systems that the Newton method needs to numerically integrate at runtime may become stiff, forcing the numerical integrators to choose very small step sizes, thus slowing down the whole optimisation process. Indeed, in these nonlinear fitting methods, the parameters are the variables. The sets of large and small parameters have to change during the process and stiffness is often caused by the presence of different time scales in the differential systems.

A promising idea would consist, in the context of systems evolving from generalized chemical reactions systems, in applying the quasi-steady state approximation technique described in section 3. The differential system to be integrated could be replaced by a reduced one. The error computation of the Newton method could then be performed over the reduced system, speeding up the overall process.

Conclusion

As a conclusion, let us give some elements of answer to a quotation of Eric Ponder (a director of Long Island Biology Association), borrowed from [START_REF] Mishra | Algebra, Automata, Algorithms & Beyond[END_REF]Postscript]: "work on the mathematics [for biology] seems to me to have developed along two equally unprofitable lines".

The use of computer algebra methods for performing differential elimination seems very promising for enhancing (at least) software for system modeling in cellular biology. The recent progresses in the quasi-steady state approximation theory suggest to widen the use of generalized chemical reactions systems and make more explicit the simplification hypotheses which would otherwise justify the use of Henri-Michaelis-Menten or Hill terms. Already available standalone C libraries for differential elimination [START_REF] Boulier | The BLAD libraries[END_REF] should still simplify this evolution.

In turn, the search of applications in the field of system modeling in cellular biology also led to improvements in differential algebra. It is this work which pushed the authors to study of the quasi-steady state approximation theory in the context of generalized chemical reactions systems and which has also suggested a new improvement to the nonlinear least squares methods.

Last, these works are far from over. For instance, an important theoretical challenge consists in designing an algorithmic method for determining good slow varieties of generalized chemical reactions systems. The authors are involved in another, more concrete challenge: to develop a modeling software based on these ideas and to use it to design and analyze synthetic gene networks.

R

  := differential_ring (derivations = [t], notation = diff, ranking = [[s, c, x, y, z]]):

2 Fig. 1 .

 21 Fig.1. A single gene regulated by a polymer of its own protein.

  eq1 := diff (x1(t),t) + k12*x1(t) -k21*x2(t) + Ve*x1(t)/(ke + x1(t)); eq2 := diff (x2(t),t) -k12*x1(t) + k21*x2(t); K := field_extension (transcendental_elements = [k21, k12, ke, Ve]): R := differential_ring (derivations = [t], notation = diff, field_of_constants = K, ranking = [x2, x1]): ideal := Rosenfeld_Groebner ([numer (eq1), eq2], [denom (eq1)], R); ideal := [characterizable]

This paper is only concerned by the ordinary case.

Emphasized words have a technical meaning, which will be defined in section 2.2.

In the differential algebra theory, the terminology differential indeterminates is preferred to dependent variables for derivations are abstract and differential indeterminates are not even assumed to correspond to functions. In order not to mix different expressions in this paper, the second expression, which seems to be more widely known, was chosen.

Chemical reactions systems are said to be generalized when their elementary reactions are not required to be balanced. See[START_REF] Horn | General mass action kinetics[END_REF].

There are more efficient algorithms than Rosenfeld-Gröbner for performing this elimination since the input system already is a characteristic set w.r.t. some orderly ranking: one could apply a change of ranking algorithm[START_REF] Ollivier | Le problème de l'identifiabilité structurelle globale : approche théorique, méthodes effectives et bornes de complexité[END_REF][START_REF] Boulier | Computing differential characteristic sets by change of ordering[END_REF] which would avoid splitting cases.