Differential algebra and QSSA methods in biochemistry

Proceedings of the 15th IFAC Symposium on System Identification Saint-Malo, France, July 6-8, 2009

Differential algebra and QSSA methods in biochemistry ⋆ Fran¸ cois Boulier ∗ Fran¸ cois Lemaire ∗∗ ∗ ∗∗

USTL - LIFL, 59655 V. d’Ascq France ([email protected]). USTL - LIFL, 59655 V. d’Ascq France ([email protected]).

Abstract: Differential algebra is an algebraic theory for differential equations (ordinary or with partial derivatives, linear or non linear). It involves an elimination theory which is completely algorithmic since the mid 1990s. In addition to the numerous applications which have been developed in control theory for more than twenty years, it was recently shown that differential elimination makes algorithmic a rigorous quasi-steady state approximation method for generalized chemical reaction systems. This result was recently succesfully applied to model reduction problems in cellular modeling. This paper presents this recent advance in connection with a new MAPLE package (with new functionalities). This package is built over the BLAD libraries, which are LGPL-licensed open source, written in the C programming language. Keywords: differential algebra ; elimination theory ; quasi-steady state approximation ; chemical reaction systems ; software 1. INTRODUCTION

2. ON DIFFERENTIAL ELIMINATION

Differential algebra has been applied in control theory for more than twenty years. The early papers are more conceptual in nature. Later, the developments of the characteristic sets theory permitted to tackle applications. See Fliess (1989); Ollivier (1990); Diop and Fliess (1991); Ljung and Glad (1994); Audoly et al. (2001); Denis-Vidal et al. (2003); Sedoglavic (2002) and references therein.

The differential elimination theory is a subtheory of the differential algebra (Ritt, 1950; Kolchin, 1973). 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 derivatives 1 and a ranking. 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.

This paper presents two applications and illustrates them using the new DifferentialAlgebra package. The applications are not new and the paragraphs which present them are borrowed from Boulier (2007); Boulier et al. (2007); Boulier and Lemaire (2007); Boulier et al. (2009). However, their illustrations based on DifferentialAlgebra are new. Indeed, it is the first time that some subalgorithms such as PARDI and NF get available for practitioners. The DifferentialAlgebra package is still under development and is thus not yet available. However, it is only composed of functions which implement a user interface for the BLAD libraries, which are available at (Boulier, 2004). The BLAD libraries are open source software libraries, written in the C programming language and licensed under the Lesser General Public License.

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 a˙ the derivative of a): ˙ ˙ d \ (1) b) = a˙ b + a b˙ . (a + b) = a˙ + 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(k1 , . . . , kr ) where the ki denote parameters that would be assumed to be algebraically

The BLAD libraries are not very convenient to use for a casual reader since one needs to write a C program to use them. This reader may then want to use the diffalg MAPLE package, whose first version was developed in 1995 by the first author. It provides an alternative way to perform some of the computations with a less powerful interface and fewer functionnalities.

⋆ This work was partially supported by the PPF Bioinfo of the USTL.

978-3-902661-47-0/09/$20.00 © 2009 IFAC

1

33

This paper is only concerned by the ordinary case.

10.3182/20090706-3-FR-2004.0261

15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 independent) and U is a finite set of dependent variables 2 . 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.

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 (Boulier and Lemaire, 2007, algorithm NF).

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 ap ∈ 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. The next theorem (Nullstellensatz) establishes a relationship between the radical differential ideal generated by a system and its solutions. Theorem 2. 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 .

The Rosenfeld-Gr¨ obner algorithm gathers as input a finite system F of differential polynomials and a ranking. It returns a finite family (possibly empty) C1 , . . . , Cr of finite subsets of K{U } \ K, called regular differential chains. Each system Ci defines a differential ideal Ci (it turns out to be a characteristic set of Ci ) in the sense that, for any f ∈ K{U }, we have f ∈ Ci iff NF(f, Ci ) = 0 . (3) The relationship with the radical A of the differential ideal generated by F is the following: A = C1 ∩ · · · ∩ Cr . (4) 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, Ci ) = 0 , 1 ≤ i ≤ r . (5) The differential ideals Ci do not need to be prime. They are however necessarily radical. The NF( · , Ci ) function permits to compute canonical representatives of the residue classes of the differential ring R/Ci . It plays the same role, in the setting of polynomial differential equations as the normal form algorithm of the Gr¨ obner bases theory, in the setting of polynomial equations. The next theorem is a variant of well-known elimination theorems. Theorem 4. Fix a ranking and split ΘU into two disjoint subsets X and Y such that any element of X is less than any element of Y w.r.t. the ranking. Let C be a regular differential chain, defining a differential ideal C. Then, for any f ∈ K[X] we have: NF(f, C) = NF(f, C ∩ K[X]) (6)

Proof. (Ritt, 1950, chap. II, §7, 11) or Boulier and Lemaire (2007). The Rosenfeld-Gr¨ obner algorithm Boulier et al. (1995) solves the membership problem to radical differential ideals. To present it, one defines the concept of ranking. Definition 3. 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 < a˙ and a < b ⇒ 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 .

Proof. Sketched. The NF algorithm computes a sequence of rational fractions fi /gi starting at f1 /g1 = f /1. The (i + 1)th rational fraction is computed from the ith by applying either Ritt’s reduction or an algebraic inverse computation. The key arguments of the proof are as follows:

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 .

(1) if the ith rational fraction lies in K(X) then the computation of the next one can only involve an element of ΘC whose leading derivative lies in X ; (2) since the ranking eliminates Y w.r.t. X, if the leading derivative of some element of ΘC lies in X, then this element itself lies in K[X] ; (3) rewriting an element of K(X) by an element of K[X] can only produce an element of K(X).

Rankings permit to define leaders. Leaders permit to use differential polynomial as rewrite (substitution) rules. Assume that f = ad v d + · · · + a1 v + a0 is a differential polynomial with leader v (the coefficients ai are themselves differential polynomials). Then the equation f = 0 can be written: ad−1 v d−1 + · · · + a1 v + a0 (2) v d −→ − · ad

3. APPLICATION TO THE QSSA FOR CHEMICAL REACTION SYSTEMS

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.

The quasi-steady state approximation (QSSA) is one of the model reduction methods widely applied by modelers Okino and Mavrovouniotis (1998); Kokotovic et al. (1999). In the general setting, there does not exist any algorithm to compute the change of variables which is necessary to

34

15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 convert a given ODE system into the standard singular perturbation form (Kokotovic et al., 1999, §1.2). Algorithm 1 DifferentialModelReduction(X˙ = N V ) Input: The initial parametric ODE system X˙ = N V derived from a generalized chemical reaction system involving n chemical species and p reactions. The stoichiometry matrix N has dimension n × p and integer entries. The vector X of dependent variables has dimension n. The vector of reaction rates V has dimension n. Its entries are power products of system parameters and dependent variables. Each reaction rate of V (hence each column of N ) is assumed to be tagged “fast” or “slow”. Output: a list of dynamical systems in the dependent variables X obtained by quasi-steady state approximation from the initial system or Fail. 1:

2:

3: 4:

5: 6: 7:

8:

Split N into two matrices Nf (columns of the fast reactions) and Ns (columns of the slow reactions). Split V into Vf and Vs so that the initial ODE model writes X˙ = Ns Vs + Nf Vf . By (say) Gaussian elimination, determine a maximal linearly independent set of columns of Nf and remove the other ones, giving N f . Update the vector of reaction rates Vf , giving a new vector V f , such that the initial ODE model writes as follows. The entries of V f are linear combinations of elements of Vf with rational number coefficients. X˙ = Ns Vs + N f V f . Build a vector F of new dependent variables, having the same dimension as V f . Build the following DAE of the differential polynomial ring K{X ∪ F } where K is the field of the rational fractions in the system parameters: X˙ = Ns Vs + N f F , V f = 0 . (7) R := a ranking F ≫ X eliminating F w.r.t. X. [C1 , . . . , Ct ] := RosenfeldGroebner((7), R) If there exists some dependent variable Xi and some regular differential chain Ck such that NF(X˙ i , Ck ) is not a rational fraction in the variables X then Fail. ˙ Ck ) for 1 ≤ k ≤ t. for Return the list X˙ = NF(X, 1 ≤ k ≤ t.

However, in the particular setting of generalized chemical reactions systems, a general method is available, provided that each reaction is tagged “fast” or “slow” Van Breusegem and Bastin (1991); Vora and Daoutidis (2001); Bennet et al. (2007). One step of this method amounts to simplify a differential algebraic equations system (DAE). It turns out that this simplification can easily be carried out by means of differential elimination as shown in Boulier et al. (2007). The DifferentialModelReduction algorithm, summarized in Algorithm 1, is borrowed from Boulier et al. (2009). It is illustrated over a famous example: the beginning of the Henri-Michaelis-Menten reduction of the basic enzymatic reaction system:

( k1 E S − k2 C ) = 0 .

The computation can be perfomed using DifferentialAlgebra as follows. The DifferentialRing and RosenfeldGroebner functions return MAPLE tables. The first one permits to

k1

k3 −−−→ E+S − ←−−−− C −−−−→ E + P .

The key assumption is that the two first reactions (on the left) are faster than the third one (on the right). The initial system of ODE writes: X˙ = N V i.e.     E˙ ! −1 1 1 k1 E S  C˙   1 −1 −1   = k2 C (9) ·  S˙  −1 1 0  k3 C 0 0 1 P˙ 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-hand side minus the number of times it occurs on the reaction left-hand side). The rate of a reaction is the product of the left-hand side species (with multiplicities) times the reaction rate constant (the parameter over the arrow). Split the stoichiometry matrix N into two matrices Nf and Ns putting the columns which correspond to fast reactions in Nf and the ones which correspond to slow reactions in Ns . Split accordingly the rows of the vector V into two vectors Vf and Vs . One gets a formula X˙ = Ns Vs + Nf Vf . Over system (9), one gets:       E˙ −1 1  1   C˙   −1   1 −1  k1 E S  = k C . (10) · ·( )+    3  S˙  −1 1 k2 C 0 0 0 1 P˙ Determine a maximal linearly independent set of columns of Nf (i.e. a basis of that matrix) and remove the other ones, giving a new matrix N f . Update the vector of reaction rates Vf , giving a new vector V f such that Nf Vf = N f V f . Over the example, removing the second column, one gets a new formula X˙ = Ns Vs + N f V f which is equivalent to formula (9):       E˙ 1 −1  C˙   −1   1  = · ( k1 E S − k2 C ) . · ( k3 C ) +   S˙  0 −1  0 1 P˙ (11) Replace the vector V f by a vector F of new dependent variables Fi . The slow variety 3 is defined by letting the entries of V f all equal to zero. The DAE to be considered for quasi-steady state approximation is (12) Vf = 0. X˙ = Ns Vs + N f F , Over the example, one gets:       E˙ 1 −1  C˙   −1   1  = · ( k3 C ) +  · ( F1 ) ,  S˙  0 −1  (13) 1 0 P˙

3

More precisely, its approximation M0 , following the notations of (Kokotovic et al., 1999, Sect. 1.4).

(8)

k2

35

15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 describe the ranking, the used notation and the system parameters. From a mathematical point of view, parameters are just viewed as dependent variables whose derivatives are zero. The second one receives four arguments: the system to be simplified ; a list of derivatives (the list of parameters) considered as algebraically independent (this tells the software not to generate algebraic relations among parameters) ; a list of polynomials which must considered as non zero (this is another way to avoid considering degenerate cases) ; and the table R. The call to RewriteRules then permits to extract, among the differential regular chain ˙ displayed equations, the one whose leading derivative is S, as a rewrite rule i.e. solved w.r.t. its leading rank.

> hypotheses := subs (P[0]=0, C[0]=0, conservation_laws); > R := DifferentialRing (blocks = [F[1], species, params], notation = diff, derivations = [t], parameters = params): > ideal := RosenfeldGroebner ([ op(sys), op(hypotheses) ], algindep = params, R); > formula := RewriteRules (ideal[1], leader=diff(S(t),t)); d formula := -- S(t) = dt 2 2 S(t) k[1] k[3] E[0] + S(t) k[1] k[2] k[3] E[0] - ------------------------------------------------------2 2 2 S(t) k[1] + 2 S(t) k[1] k[2] + k[1] k[2] E[0] + k[2]

> with (DifferentialAlgebra): > sys := [ diff(E(t),t) = -F[1](t) + k[3]*C(t), diff(S(t),t) = - F[1](t), diff (C(t),t) = F[1](t) - k[3]*C(t), diff (P(t),t) = k[3]*C(t), k[1]*E(t)*S(t) - k[2]*C(t) = 0 ]:

Then denote K = k2 /k1 and Vmax = k3 E0 and neglect the term K E0 (assuming S0 ≫ E0 ). One gets:

> species := [C,E,P,S]; > params := [k[1],k[2],k[3]]; > R := DifferentialRing (blocks = [F[1], species, params], notation = diff, derivations = [t], parameters = params):

> > > >

formula formula formula formula

:= := := :=

subs (k[2]=K*k[1], formula): normal (formula): algsubs (k[3]*E[0]=V[max], formula): normal (subs (K*E[0]=0, formula));

d V[max] S(t) formula := -- S(t) = - ----------dt S(t) + K

R := DRing

The above method applies to larger examples such as the ones featured in Boulier et al. (2008). However, specially dedicated implementations of Algorithm 1, are today necessary to solve the largest cases.

> ideal := RosenfeldGroebner (sys, algindep = params, nonzero = [C(t)], R); ideal := [regchain] > RewriteRules (ideal[1], leader=diff(S(t),t)); d -- S(t) = dt E(t) S(t) k[1] k[3] + E(t) S(t) k[1] k[2] k[3] - -----------------------------------------------2 E(t) k[1] k[2] + S(t) k[1] k[2] + k[2]

4. APPLICATION TO THE PARAMETERS ESTIMATION PROBLEM This section describes an application of differential elimination and, more precisely, an application of algorithms which perform changes of rankings over regular differential chains. The principle of this application was designed Denis-Vidal et al. (2003).The addressed problem is this one: estimate parameters values of parametric ordinary differential systems the dependent variables of which are not all observed. The work of Joly-Blanchard, Denis-Vidal and Noiret is strongly related to the problem of the identifiability study of differential systems, for which a huge literature is available. The method of Joly-Blanchard, Denis-Vidal and Noiret is original for two reasons: it relies on rigorous differential elimination methods and it carries out the study of real examples up to the final numerical treatment. It mixes symbolics and numerics.

˙ It The above rewrite rule provides the normal form of S. is free of F1 . Let us check moreover that the normal forms of the derivatives of the other dependent variables are also free of F1 . > NF := NormalForm ([diff(E(t),t), diff(C(t),t), diff(P(t),t)], ideal[1]): > Indets (NF, R); {k[1], k[2], k[3], E(t), S(t)}

To recover the well-known Henri-Michaelis-Menten formula, some extra assumptions, which have nothing to do with QSSA, must be taken into account. One uses conservation laws (that can be deduced from the stoichiometry matrix) and some hypotheses on initial values: P0 = C0 = 0. For simplicity, the whole process (including QSSA) is run again:

4.1 Statement of the Problem over an Example Figure 1 represents a compartmental model. The two compartments represent the blood and some organ. A medical product is injected in the blood at t = 0. It can go from the blood to the organ and conversely. It may also get degraded and exit from the system. In order to write the corresponding differential system, some hypotheses must be made on the nature of the exchanges: exchanges between the two compartments are assumed to be linear.

> conservation_laws := [E(t) + C(t) = E[0] + C[0], S(t) + C(t) + P(t) = S[0] + C[0] + P[0]]; > params := [k[1],k[2],k[3],E[0],C[0],S[0],P[0]];

36

15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 afterwards by the purely numerical method as a starting point.

The degradation is assumed to follow a Michaelis-Menten law. (linear exchange) (Michaelis–Menten exchange) Ve , ke

The idea here consists in eliminating the non observed variables of the model. In other words, the idea consists in computing a differential polynomial which lies in the differential ideal generated by the model equations and which only involves the observed variable x1 , its derivatives up to any order and the model parameters. Let us show how to do this with DifferentialAlgebra.

k12

compartment 1

compartment 2 organ

blood

k21 (linear exchange)

The differential system actually is a regular differential chain w.r.t. some orderly ordering. The next call to RosenfeldGroebner thus does not perform any relevant computation. Observe that the first ODE is a rational fraction. Its numerator is kept among the equations to be processed while its denominator is automatically inserted in the nonzero list.

Fig. 1. Compartmental model From Figure 1, it is possible to derive a system of parametric ordinary differential equations. One associates to compartments 1 and 2, dependent variables x1 and x2 which represent the concentrations of product present in these compartments. For simplicity, it is assumed here that both compartments have a unitary volume. One gets the following differential system. Ve x1 x˙ 1 = −k12 x1 + k21 x2 − , (14) ke + x1 x˙ 2 = k12 x1 − k21 x2 .

> with (DifferentialAlgebra); > pars := [k[e], V[e], k[12], k[21]]; > R := DifferentialRing (derivations = [t], blocks = [[x[1],x[2]], pars], parameters = pars, notation = diff); R := DRing

Let us consider now some instance of the above model and assume that some extra information is available: parameters k12 , k21 and Ve are assumed to be unknown while ke is assumed to be fixed. Compartment 1 is assumed to be observed i.e. a file of measures is assumed to be available for x1 . Compartment 2 is assumed to be non observed. One just assumes x2 (0) = 0 i.e. that no product is initially present in the organ. We are now ready to state the problem over this example: given the system of parametric ordinary differential equations, the file of measures and the extra information, estimate the values of the three unknown parameters: Ve , k12 and k21 .

> edoA := diff(x[1](t),t) = -k[12]*x[1](t) + k[21]*x[2](t) - (V[e]*x[1](t))/(k[e]+x[1](t)); > edoB := diff(x[2](t),t) = k[12]*x[1](t) - k[21]*x[2](t); edoA := d V[e] x[1](t) -- x[1](t) = -k[12] x[1](t) + k[21] x[2](t) - -------------dt k[e] + x[1](t) d edoB := -- x[2](t) = k[12] x[1](t) - k[21] x[2](t) dt > ideal := RosenfeldGroebner([edoA,edoB], algindep=pars, R);

4.2 The Numerical Method ideal := [regchain]

There exists a purely numerical method to solve this problem. It is a non linear least squares solving method i.e. a Newton method. Precisely, a Levenberg-Marquardt solver is called. 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 Levenberg-Marquardt method 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 any Newton-like method, the Levenberg-Marquardt is very likely to end up in a local minimum and produce wrong parameters values.

Eliminating the non observed variable could be performed by RosenfeldGroebner. However, it is more efficient to apply PARDI (Boulier et al., 2001) which, following Ollivier (1990), takes advantage from the fact that the input equations form a regular differential chain. PARDI takes two arguments: an existing regular differential chain and a target ranking. It returns a regular differential chain w.r.t. the target ranking, defining the same differential ideal as the first argument. > idealES := Pardi (ideal [1], ordering (derivations = [t], blocks = [x[2],x[1],pars])); idealES := regchain > Equations (idealES, leader=derivative(x[1](t)));

4.3 The Symbolic-Numeric Method

The equation returned by the last command above can be pretty-printed as follows:

The previous method has a drawback: it relies on nonlinear least squares which require the a priori knowledge of a good approximation of the parameters values. Thanks to differential elimination and to linear least squares, it is possible to estimate a first approximation of the parameters values. This first approximation may be used

x ¨1 (x1 + ke )2 + [k12 + k21 ] x˙ 1 (x1 + ke )2 + [Ve ] x˙ 1 ke + [k21 Ve ] x1 (x1 + ke ) = 0.

(15)

Since the target ranking eliminates x2 w.r.t. x1 , the equation whose leader is a derivative of x1 is necessarily free

37

15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 of x2 , which is the case for equation (15). The expressions enclosed between square brackets are called parameters blocks. Equation (15) tells us that the model is globally identifiable i.e. that, given a function x1 and a parameter value ke , the three unknown parameters are uniquely defined. Indeed, assume that the function x1 is known. Then so are its derivatives x˙ 1 and x ¨1 . These three functions can therefore be evaluated for three different values of the time t. The known parameter ke 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 k12 , k21 and Ve also are uniquely defined. QED. In practice, the function x1 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. One can then recover the parameters values from the blocks values. These values, which are not very accurate, can however be used as a starting point for the purely numerical method. 4.4 Issues In general, there is no guarantee that the first estimation provided by the symbolic-numeric method leads the purely numerical method in the global minimum. Estimating parameters only makes sense for models at least locally identifiable. However, testing this property does not raise any difficulty. Some seminumerical algorithms are available Sedoglavic (2002). These are probabilistic algorithms for which the failure probability is known and can be decreased up to any value. Numerically estimating the derivatives raises an important difficulty. To overcome it, a good method consists in converting the differential equations as integral equations. Under some conditions, integral equations are less sensitive to the noise than differential equations. There exists another important difficulty: 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: [Ve ], [k21 ], [Ve k21 ]. 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. REFERENCES Audoly, S., Bellu, G., D’Angio, L., Saccomani, M.P., and Cobelli, C. (2001). Global Identifiability of Nonlinear Models of Biological Systems. IEEE Transactions on Biomedical Engineering, 48(1), 55–65.

38

Bennet, M.R., Volfson, D., Tsimring, L., and Hasty, J. (2007). Transient Dynamics of Genetic Regulatory Networks. Biophysical Journal, 92, 3501–3512. Boulier, F. (2004). The BLAD libraries. http://www.lifl.fr/ ∼boulier/BLAD. Boulier, F. (2007). Differential Elimination and Biological Modelling. Radon Series on Computational and Applied Mathematics (Gr¨ obner Bases in Symbolic Analysis), 2, 111–139. http://hal. archives-ouvertes.fr/hal-00139364. Boulier, F. and Lemaire, F. (2007). A computer scientist point of view on Hilbert’s differential theorem of zeros. Submitted to Applicable Algebra in Engineering, Communication and Computing. http://hal.archives-ouvertes.fr/hal-00170091. Boulier, F., Lazard, D., Ollivier, F., and Petitot, M. (1995). Representation for the radical of a finitely generated differential ideal. In ISSAC’95: Proceedings of the 1995 international symposium on Symbolic and algebraic computation, 158–166. ACM Press, New York, NY, USA. http://hal.archives-ouvertes.fr/ hal-00138020. Boulier, F., Lefranc, M., Lemaire, F., and Morant, P.E. (2007). Model Reduction of Chemical Reaction Systems using Elimination. Presented at the international conference MACIS 2007, http://hal.archives-ouvertes.fr/hal-00184558, submitted to Mathematics in Computer Science. Boulier, F., Lefranc, M., Lemaire, F., and Morant, P.E. (2008). Applying a rigorous quasi-steady state approximation method for proving the absence of oscillations in models of genetic circuits. In K.H. et al. (ed.), Proceedings of Algebraic Biology 2008, number 5147 in LNCS, 56–64. Springer Verlag Berlin Heidelberg. Boulier, F., Lemaire, F., and Maza, M.M. (2001). PARDI! In ISSAC’01: Proceedings of the 2001 international symposium on Symbolic and algebraic computation, 38–47. ACM Press, New York, NY, USA. http://hal.archives-ouvertes.fr/ hal-00139354. ¨ upl¨ Boulier, F., Lemaire, F., Urg¨ u, A., and Sedoglavic, A. (2009). Towards an automated reduction method for polynomial ODE models in cellular biology. Mathematics in Computer Science, Special Issue on Symbolic Computation in Biology. To appear. Denis-Vidal, L., Joly-Blanchard, G., and Noiret, C. (2003). System identifiability (symbolic computation) and parameter estimation (numerical computation). In Numerical Algorithms, volume 34, 282–292. Diop, S. and Fliess, M. (1991). Nonlinear observability, identifiability, and persistent trajectories. In Proc. 30th CDC, 714–719. Brighton. Fliess, M. (1989). Automatique et corps diff´ erentiels. Forum Math., 1, 227–238. Kokotovic, P., Khalil, H.K., and O’Reilly, J. (1999). Singular Perturbation Methods in Control: Analysis and Design. Classics in Applied Mathematics 25. SIAM. Kolchin, E.R. (1973). Differential Algebra and Algebraic Groups. Academic Press, New York. Ljung, L. and Glad, S.T. (1994). On global identifiability for arbitrary model parametrisations. Automatica, 30, 265–276. Okino, M.S. and Mavrovouniotis, M.L. (1998). Simplification of Mathematical Models of Chemical Reaction Systems. Chemical Reviews, 98(2), 391–408. Ollivier, F. (1990). Le probl` eme de l’identifiabilit´ e structurelle globale : approche th´ eorique, m´ ethodes effectives et bornes de ´ complexit´ e. Ph.D. thesis, Ecole Polytechnique, Palaiseau, France. Ritt, J.F. (1950). Differential Algebra. Dover Publications Inc., New York. http://www.ams.org/online bks/coll33. Sedoglavic, A. (2002). A Probabilistic Algorithm to Test Local Algebraic Observability in Polynomial Time. Journal of Symb. Comp., 33(5), 735–755. Van Breusegem, V. and Bastin, G. (1991). Reduced order dynamical modelling of reaction systems: a singular perturbation approach. In Proceedings of the 30th IEEE Conference on Decision and Control, 1049–1054. Brighton, England. Vora, N. and Daoutidis, P. (2001). Nonlinear model reduction of chemical reaction systems. AIChE Journal, 47(10), 2320–2332.