Application of Differential Algebra to the Quasi-Steady State - - PowerPoint PPT Presentation

Application of Differential Algebra to the Quasi-Steady State Approximation in Biology and Physics

François Lemaire

Université de Lille I (France) Symbolic Computation Team (Boulier, Oussous, Petitot, Sedoglavic)

Differential Algebra and Related Topics 2010, Beijing Supported by the ANR LÉDA (Logistique des Équations Différentielles Algébriques)

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 1 / 49

Plan

1

Background Differential elimination Slow/fast dynamics Tikhonov theorem

2

Slow/fast chemical reaction systems reduction Michaelis Menten example General method

3

Application to physics Communicating vessels Diffusion Pendulum Other examples ...

4

MABSys

5

Conclusion

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 2 / 49

Summary

1

Background Differential elimination Slow/fast dynamics Tikhonov theorem

2

Slow/fast chemical reaction systems reduction Michaelis Menten example General method

3

Application to physics Communicating vessels Diffusion Pendulum Other examples ...

4

MABSys

5

Conclusion

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 3 / 49

Differential Algebra

founders of differential algebra: Ritt (50), Kolchin (73) many people in this room have contributed to differential algebra key notions:

ranking → an ordering on the derivatives characteristic sets (regular differential chains, characterizable sets) → membership test to an ideal

package for computing the characteristic sets: diffalg package (Boulier, Hubert, Lemaire) evolution of diffalg: → Differential Algebra (Boulier, Maple 14), based on BLAD (Boulier, C library, GPL)

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 4 / 49

Differential elimination

System (diff notation)

dC(t) dt = F1(t) − k2C(t) dS(t) dt = −F1(t) dE(t) dt = −F1(t) + k2C(t) dP(t) dt = k2C(t) k1E(t)S(t) = k−1C(t)

Symbols

differential indeterminates: C(t), S(t), E(t), P(t), F1(t) parameters (=constants): k1, k−1, k2

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 5 / 49

Differential elimination

System (dotted notation)

˙ C = F1 − k2C ˙ S = −F1 ˙ E = −F1 + k2C ˙ P = k2C k1E S = k−1C

Symbols

differential indeterminates: C, S, E, P, F1 derivatives: C, ˙ C, S, ˙ S, E, ˙ E, P, ˙ P, F1 ranking [F1] ≫ [C, E, P, S]: any derivative of F1 is greater than any derivative of C, S, E, P. → elimination ranking (which eliminates F1) parameters (=constants): k1, k−1, k2 In this talk, all equations are ordinary (i.e. no partial derivatives)

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 6 / 49

Rosenfeld-Gröbner algorithm

compute a list of r.d.c. C1, . . . , Cs from a input system Σ p [Σ] = ∩s

i=1Sat(Ci) with Sat(Ci) = [Ci] : HCi ∞

each r.d.c. yields a rewritting system such that p ∈ Sat(Ci) ⇐ ⇒ p

∗

− →

Ci

C → E S k1 k−1 ˙ S → − k2k1ES(k1S + k−1) k−1(k−1 + k1S + k1E) ˙ P → . . . ˙ E → . . . F1 → k2k1ES(k1S + k−1) k−1(k−1 + k1S + k1E) ranking: · · · > ˙ F1 > F1 > · · · > ˙ C > ˙ E > ˙ P > ˙ S > C > E > P > S

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 7 / 49

Normal forms

• ne has a normal form for polynomials (BL Issac01):

p = q mod Sat(Ci) ⇐ ⇒ NF(p, Ci) = NF(q, Ci) the normal form of a polynomial is a rational fractions f/g recent paper : A Normal Form Algorithm For Regular Differential Chains, → Boulier Lemaire, AADIOS09MCS

• ne can also consider normal forms of a rational fraction (provided its denominator

is not a zero divisor).

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 8 / 49

Timescales

Different timescales

in physics, biology, . . . : phenomena can have very different timescales exemple: communicating vessels with an input

x2 x1 u V1,S1 V2,S2

Two phenomena: input: u(t) (m3/s) exchange between 1 and 2 Variables: volumes: V1(t), V2(t) (in m3) sections: S1, S2 (in m2) water heights: x1(t), x2(t) (in m)

Objective

Get rid of the fast timescale

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 9 / 49

Removing the fast timescale (1/2)

First case: assume the input is much slower than the exchange

the two vessels balance instantly

• ne can assume x1(t) = x2(t) (after some transient time)

the input is split between the two vessels w.r.t. their surfaces

x2 x1 u V1,S1 V2,S2

If one assumes it starts balanced: ˙ x1 = 1 S1 + S2 u ˙ x2 = 1 S1 + S2 u x1(0) = x2(0)

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 10 / 49

Removing the fast timescale (2/2)

Second case: assume the input is much faster than the exchange

the two vessels have no time for balancing the input only goes in vessel 1

x2 x1 u V1,S1 V2,S2

the vessel 1 will be full in a very short time removing the fast timescale here makes no sense indeed, the fast phenomena does not reach a equilibrium

• ne expects the fast phenomena to

balance

Remark

One can easily remove the slow timescale: ˙ x1 = u/S1, ˙ x2 = 0

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 11 / 49

Intuition of the Tikhonov theorem

Σ  ˙ x = f(x, y) ˙ y =

1 εg(x, y)

ε is small x is the slow variable y is the fast variable (because of 1/ε) a solution starting from (x0, y0) behaves in two steps: fast transient x(t) does not change y(t) quickly reaches the curve g(x0, y) = 0 slow the solution slowly drifts on g(x, y) = 0 following ˙ x = f(x, y). The variety g(x, y) = 0 is called the slow variety

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 12 / 49

Example of the Tikhonov theorem

˙ x = x + y ˙ y =

1 ε(1 − xy)

0.5 1 1.5 2 2.5 3 y 2 4 6 8 10 x

for ε = 0.01 and three different initial conditions (x0, y0) = (2, 3), (3, 1) and (1, 0.5). The transient is (almost) vertical, and solutions slides along y = 1/x.

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 13 / 49

Tikhonov theorem

Idea of the theorem

it is a limit theorem: the limits of the solutions (when ε → 0) are the solutions when ε = 0 in the system.

The theorem (rough version)

Σ  ˙ x = f(x, y) ˙ y =

1 εg(x, y)

If for any (x0, y0), the solution of ˙ y = 1

εg(x0, y), y(0) = y0 converges to ¯

y with g(x0, ¯ y) = 0, then there exist α, T > 0 such that the solution of Σ uniformally tends (for t ∈ [α, T]) towards the solution of the reduced system ¯ Σ  ˙ x = f(x, y) = g(x, y) ¯ Σ is the quasi steady-state approximation of Σ.

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 14 / 49

Tikhonov for the vessels

x2 x1 u V1,S1 V2,S2

Input slower than balancing ˙ x1 = u/S1 − 1

εF(x1, x2)/S1

˙ x2 =

1 εF(x1, x2)/S2

x1(0) = x2(0) the exchange satisfies: F(x1, x2) = 0 ⇐ ⇒ x1 = x2 if you assume that x1 = x2, you get: ˙ x1 = u/S1 and ˙ x2 = 0 pb: 1/ε times something small is not zero exact algorithms do not like ≈ 0

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 15 / 49

Tikhonov for the vessels

x2 x1 u V1,S1 V2,S2

Input slower than balancing ˙ x1 = u/S1 − 1

εF(x1, x2)/S1

˙ x2 =

1 εF(x1, x2)/S2

x1(0) = x2(0) The theorem does not apply (not in Tikhonov form). x1 and x2 are both "fast". introducing y = S1x1 + S2x2, the system is under Tikhonov form (y is slow . . . ) indeed, y is the total volume of liquid, so it is not affected by the exchange phenomenon ˙ y = u ˙ x2 =

1 εF((y − S2x2)/S1, x2)/S2

since x1 = (y − S2 x2)/S1 Applying the Tikhonov theorem . . .

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 16 / 49

Tikhonov for the vessels

x2 x1 u V1,S1 V2,S2

Input slower than balancing ˙ x1 = u/S1 − 1

εF(x1, x2)/S1

˙ x2 =

1 εF(x1, x2)/S2

x1(0) = x2(0) The theorem does not apply (not in Tikhonov form). x1 and x2 are both "fast". introducing y = S1x1 + S2x2, the system is under Tikhonov form (y is slow . . . ) indeed, y is the total volume of liquid, so it is not affected by the exchange phenomenon ˙ y = u = F((y − S2x2)/S1, x2)/S2 Since F(x1, x2) = 0 ⇐ ⇒ x1 = x2, and replacing y by its value

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 16 / 49

Tikhonov for the vessels

x2 x1 u V1,S1 V2,S2

Input slower than balancing ˙ x1 = u/S1 − 1

εF(x1, x2)/S1

˙ x2 =

1 εF(x1, x2)/S2

x1(0) = x2(0) The theorem does not apply (not in Tikhonov form). x1 and x2 are both "fast". introducing y = S1x1 + S2x2, the system is under Tikhonov form (y is slow . . . ) indeed, y is the total volume of liquid, so it is not affected by the exchange phenomenon S1 ˙ x1 + S2 ˙ x2 = u x1 = x2 and finally. . .

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 16 / 49

Tikhonov for the vessels

x2 x1 u V1,S1 V2,S2

Input slower than balancing ˙ x1 = u/S1 − 1

εF(x1, x2)/S1

˙ x2 =

1 εF(x1, x2)/S2

x1(0) = x2(0) The theorem does not apply (not in Tikhonov form). x1 and x2 are both "fast". introducing y = S1x1 + S2x2, the system is under Tikhonov form (y is slow . . . ) indeed, y is the total volume of liquid, so it is not affected by the exchange phenomenon ˙ x1 = 1 S1 + S2 u ˙ x2 = 1 S1 + S2 u

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 16 / 49

Elimination can help you

if you do not want to find the change of variables → differential elimination is for you:

x2 x1 u V1,S1 V2,S2

˙ x1 = u/S1 − F1/S1 ˙ x2 = F1/S2 x1 = x2 we do as if we don’t know the law for balancing, but only the law describing the balanced state (x1 = x2) → indeed, the slow system does not depend on the balancing law the elimination of F1 combines the three steps:

finding the change of variables applying the tikhonov theorem in the new variables express the system in the original variables

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 17 / 49

Elimination can help you

if you do not want to find the change of variables → differential elimination is for you:

x2 x1 u V1,S1 V2,S2

˙ x1 = u/S1 − F1/S1 ˙ x2 = F1/S2 x1 = x2 ˙ x1 → u S1 + S2 x2 → x1 F1 → u S2 S1 + S2 ranking: · · · > ˙ F1 > F1 > · · · > ˙ x2 > ˙ x1 > x2 > x1

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 17 / 49

Summary

1

Background Differential elimination Slow/fast dynamics Tikhonov theorem

2

Slow/fast chemical reaction systems reduction Michaelis Menten example General method

3

Application to physics Communicating vessels Diffusion Pendulum Other examples ...

4

MABSys

5

Conclusion

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 18 / 49

Statement of the problem

Input

a set of chemical reactions some are really fast, the others are slow

Consequence

the fast reactions tend to be (almost) at equilibria two timescales: a fast one for the fast reactions, a slow one for the slow ones the fast timescale yields stiff ode systems which are hard to integrate (the integrator tends to oscillate around the correct solution) some variables are not needed

Objective

Eliminate the fast timescale

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 19 / 49

Our algorithm ModelReduce

completely algorithmic (based on differential elimination) available in MABSys (Maple package, Lemaire/Ürgüplü) makes algorithmic ideas from [1, 2, 3].

[1] V. Van Breusegem and G. Bastin. Reduced order dynamical modelling of reaction systems: a singular perturbation approach 30th IEEE Conf. on Decision and Control. pp. 1049-1054, 1991. [2] N. Vora and P . Daoutidis. Nonlinear model reduction of chemical reaction systems. AIChE Journal vol. 47, pp. 2320-2332, 2001. [3] M. Bennet, D. Volfson, L. Tsimring and J. Hasty. Transient Dynamics of Genetic Regulatory Networks Biophysical Journal vol. 92, pp. 3501-3512, 2007 [4] F . Boulier, M. Lefranc, F . Lemaire and P .-E. Morant. Model Reduction of Chemical Reaction Systems using Elimination. MACIS, http://hal.archives-ouvertes.fr/hal-00184558, 2007.

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 20 / 49

A classical enzymatic degradation

The reactions

E + S

k1

− − ⇀ ↽ − −

k−1

C (supposed fast) C

k2

− → E + P (supposed slow) substrate S, product P, complex C, enzyme E

Consequence

E + S − ⇀ ↽ − C is always (almost) at equilibria by “ignoring” E and C, can the system be viewed as S

?

− → P ?

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 21 / 49

A classical enzymatic degradation

The usual reduction

˙ S(t) = − Vm S(t) K + S(t) assuming S ≫ E0 Vm = k2 E0 Briggs-Haldane: K =

k−1 k1

Henri-Michaëlis-Menten: K =

k−1+k2 k1

By Hand

Our reduction

˙ S = − Vm S (K + S) K E0 + (K + S)2 seems valid even if S < E0 Vm = k2 E0 K =

k−1 k1

Automatic

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 22 / 49

A classical enzymatic degradation

The usual reduction

˙ S(t) = − Vm S(t) K + S(t) assuming S ≫ E0 Vm = k2 E0 Briggs-Haldane: K =

k−1 k1

Henri-Michaëlis-Menten: K =

k−1+k2 k1

By Hand

Our reduction

˙ S = − Vm S (K + S) K E0 + (K + S)2 seems valid even if S < E0 Vm = k2 E0 K =

k−1 k1

Automatic

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 22 / 49

From the reactions to a dynamical system

R1 : E + S

k1

− → C R2 : C

k−1

− − → E + S R3 : C

k2

− → E + P Dynamical system: ˙ X = MV Vector of concentrations: X = B B @ C E P S 1 C C A Reaction rate (following mass action law): V = @ k1 E S k−1C k2C 1 A Stoichiometric matrix: M = B B @ 1 −1 −1 −1 1 1 1 −1 1 1 C C A

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 23 / 49

Detail of the reduction

E + S

k1

− − ⇀ ↽ − −

k−1

C (supposed fast) C

k2

− → E + P (supposed slow)

Step 1: build the system

˙ C = F1 − k2C ˙ S = −F1 ˙ E = −F1 + k2C ˙ P = k2C k1ES = k−1C (eq) F1: contribution of the fast reaction (unknown for the moment) (eq) implies that E + S − ⇀ ↽ − C is at equilibria

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 24 / 49

Detail of the reduction

E + S

k1

− − ⇀ ↽ − −

k−1

C (supposed fast) C

k2

− → E + P (supposed slow)

Step 2: eliminate the F1 with Rosenfeld–Gröbner

˙ S = − k2k1ES(k1S + k−1) k−1(k−1 + k1S + k1E) ˙ C = . . . ˙ S = . . . ˙ E = . . . F1 = . . .

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 24 / 49

Detail of the reduction

E + S

k1

− − ⇀ ↽ − −

k−1

C (supposed fast) C

k2

− → E + P (supposed slow)

Step 3: use conservations laws

˙ S = − k1k2E0 S (k1S + k−1) k1k−1 E0 + (k1S + k−1)2 ˙ C = . . . ˙ S = . . . ˙ E = . . . using the conservation laws: C + P + S = C0 + P0 + S0 C + E = C0 + E0 and assuming C0 = 0.

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 24 / 49

Detail of the reduction

E + S

k1

− − ⇀ ↽ − −

k−1

C (supposed fast) C

k2

− → E + P (supposed slow)

Step 4: use symmetries

˙ S = − k1k2E0 S (k1S + k−1) k1k−1 E0 + (k1S + k−1)2 is reformulated to ˙ S = − Vm S (K + S) K E0 + (K + S)2 with Vm = k2 E0 and K =

k−1 k1

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 24 / 49

ModelReduce algorithm

Input

a set of chemical reactions tagged slow and fast some options (conservation laws, ordering, . . . )

Output

a list of dynamical systems (i.e. ˙ X = · · · ) multiple output occurs when the slow variety is split

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 25 / 49

ModelReduce algorithm

Steps

introduce a new symbol Fi ∈ F for each fast reaction build the stoichiometric matrix Mf and the vector rate Vf for the fast reactions do the same for the slow reactions (yields Ms and Vs) build Σ  ˙ X = MsVs + MfF (dynamical system) MfVf = (fast reactions at equilibria) eliminate all the Fi’s and obtain one or several r.d.c. C1, C2, . . . for each Ci, output the dynamical system [. . . , ˙ xj = NF( ˙ xj, Ci), . . .] When the Tikhonov theorem does not apply, the output is not correct (in that case NF( ˙ xj, Ci) might involve 1st-order derivatives, or some Fi)

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 26 / 49

Summary

1

Background Differential elimination Slow/fast dynamics Tikhonov theorem

2

Slow/fast chemical reaction systems reduction Michaelis Menten example General method

3

Application to physics Communicating vessels Diffusion Pendulum Other examples ...

4

MABSys

5

Conclusion

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 27 / 49

Communicating vessels: the system

x2 x1 u V1,S1 V2,S2

Variables

Volumes: V1(t), V2(t) (in m3) Sections: S1, S2 (in m2) Water heights: x1(t), x2(t) (in m) Input in the vessel 1: u(t) (in m3/s)

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 28 / 49

Communicating vessels: the system

x2 x1 u V1,S1 V2,S2

Hypothesis

the filling is much slower than the balancing between the two vessels the dynamics of the balancing law needs to to be known

• ne just assumes the balancing is done when x1 = x2

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 28 / 49

Communicating vessels: the reduction

Equations

˙ V1(t) = u(t) + F1(t) ˙ V2(t) = −F1(t) V1(t) = S1x1(t) V2(t) = S2x2(t) x1(t) = x2(t)

After eliminating F1

˙ x1 = u S1 + S2 For information: F1 = −

S2u S1+S2

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 29 / 49

Fast diffusion: system

2 1 u

Variables

Volumes: V1, V2 (in l) number of moles of a species X: n1(t), n2(t) (in mol) concentration of X: x1(t), x2(t) (in mol/l) addition of X in compartment 1: u(t) (in mol/s)

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 30 / 49

Fast diffusion: system

2 1 u

Hypothesis

the addition of X is much slower than the balancing between the two compartments the diffusion law needs not to be known

• ne just assumes that the diffusion preserves the total number of X and satisfies

x1 = x2

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 30 / 49

Fast diffusion: reduction

Equations

˙ n1(t) = u(t) + F1(t) ˙ n2(t) = −F1(t) n1(t) = V1x1(t) n2(t) = V2x2(t) x1(t) = x2(t)

After eliminating F1

˙ x1 = u V1 + V2 For information: F1 = −

V2u V1+V2

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 31 / 49

The spring pendulum: system

Ty Tx mg x y

Variables

coordinates: x(t), y(t) (in m) mass: m (in kg) spring constant: k (in kg/s2) gravitation: g (in m/s2) free length of the spring : l0

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 32 / 49

The spring pendulum: system

Ty Tx mg x y

Hypothesis: the spring has an infinite spring constant

its length equals the free length (the gravitation and the motion do not stretch it) the spring applies a force (Tx, Ty) (supposed unknown) collinear to the spring direction caution: Tx and Ty are two different functions (no partial derivatives here)

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 32 / 49

Spring pendulum: reduction

Equations

m¨ x(t) = −Tx(t) m¨ y(t) = mg − Ty(t) x(t)2 + y(t)2 = l2 x(t)Ty(t) = y(t)Tx(t)

After eliminating Tx and Ty

y 2 + x2 = l2 ¨ x = ˙ x2 2 „ 1 l0 + x − 1 l0 − x « − yxg l2

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 33 / 49

Singularities

Horizontal position

The following system has a singularity when x(t) = l0 or x(t) = −l0 y 2 + x2 = l2 ¨ x = ˙ x2 2 „ 1 l0 + x − 1 l0 − x « − yxg l2

Vertical position

The following system has a singularity when y(t) = l0 or y(t) = −l0 y 2 + x2 = l2 ¨ y = g − gy 2 l2 − ˙ y 2 2 „ 1 l0 + y − 1 l0 − y «

Numerical considerations

numerical consequence: how to integrate ? By switching model ? linked to our project LEDA (ANR=french research agency), which deals with the treatment of DAE

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 34 / 49

Communicating vessels with varying section

The section is a function of the height x: V1(x1) = R x1

0 S(x)dx.

One needs to solve the problem with compositions of functions

System

˙ V1(x1(t)) = u(t) + F1(t) ˙ V2(x2(t)) = −F1(t) V1(x1(t)) = Z x1(t) S(x)dx V2(x2(t)) = Z x2(t) S(x)dx x1(t) = x2(t) Hard to do because of composition of unknown functions. Easy by hand, not in diffalg

• r DifferentialAlgebra.

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 35 / 49

Summary

1

Background Differential elimination Slow/fast dynamics Tikhonov theorem

2

Slow/fast chemical reaction systems reduction Michaelis Menten example General method

3

Application to physics Communicating vessels Diffusion Pendulum Other examples ...

4

MABSys

5

Conclusion

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 36 / 49

Presentation

Implementation

Authors: François Lemaire and Aslı Ürgüplü (2008-) Coded in Maple Download: www.lifl.fr/~urguplu

Description

handles : models described by chemical reactions provides : routines for reducing/simplifying the model (ODEs) goal : make algorithmic the reductions usually made by hand help the analysis of the model spirit : accessible to non specialists, inspired by collaboration with modelers

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 37 / 49

Typical use of MABSys

System of chemical reactions system of ODEs

user Rate laws or QSSA (App. Reduction)

system of ODEs

Exact Reduction (parameter reduction/semi-rectification) user

Qualitative/quantitative Information

Outside MABSys (User, Extensive computations, Bifurcation Analysis, . . . )

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 38 / 49

Underlying techniques

Techniques

differential elimination for the QSSA

• F. Boulier, M. Lefranc, F. Lemaire, and P

.-E. Morant. Model Reduction of Chemical Reaction Systems using Elimination. MACIS, 2007. http://hal.archives-ouvertes.fr/hal-00184558/fr.

Lie Symmetries (parameter reduction, semi-rectification)

• F. Lemaire and A. Ürgüplü.

A Method for Semi-Rectifying Algebraic and Differential Systems using Scaling type Lie Point Symmetries with Linear Algebra. In Proceedings of ISSAC, 2010. P . J. Olver. Applications of Lie groups to differential equations, second ed., vol. 107 of Graduate Texts in Mathematics. Springer Verlag, 1993.

• A. Sedoglavic.

Reduction of Algebraic Parametric Systems by Rectification of their Affine Expanded Lie Symmetries. In K. Horimoto H. Anai and T. Kutsia, editors, Proceedings of Algebraic Biology 2007, volume 4545 of LNCS, pages 277–291, 2007.

Techniques are hidden to the user ! We want a friendly interface.

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 39 / 49

Model definition

Basic enzymatic degradation

E + S

k1

− − ⇀ ↽ − −

k−1

C (supposed fast) C

k2

− → E + P (supposed slow) substrate S, product P, complex C, enzyme E

> R1 := NewReaction(E+S,C,MassActionLaw(k1),fast=true): > R2 := NewReaction(C,E+S,MassActionLaw(km1),fast=true): > R3 := NewReaction(C,E+P,MassActionLaw(k2)): > RS := [R1,R2,R3]:

Remark: one can use CustomizedLaw for arbitrary rates

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 40 / 49

Basic operations

E + S

k1

− − ⇀ ↽ − −

k−1

C (supposed fast) C

k2

− → E + P (supposed slow)

Vector of rates

> RateVector(RS); [k1 E S] [ ] [km1 C ] [ ] [ k2 C ]

Stoechiometry matrix

> StoichiometricMatrix(RS, [E,S,C,P]); [-1 1 1] [ ] [-1 1 0] [ ] [ 1

• 1
• 1]

[ ] [ 0 1]

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 41 / 49

Basic operations

Conversion to an ODE system

> ReactionSystem2ODEs(RS, [E,S,C,P]); d [-- E(t) = -k1 E(t) S(t) + km1 C(t) + k2 C(t), dt d

• - S(t) = -k1 E(t) S(t) + km1 C(t),

dt d

• - C(t) = k1 E(t) S(t) - km1 C(t) - k2 C(t),

dt d

• - P(t) = k2 C(t)]

dt

Steady points equations

> Equilibria(RS); [k1 E S - km1 C - k2 C, k2 C]

Other basic routines: access reactions information, ODE simulation and plottings,. . .

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 42 / 49

ModelReduce

no conservation laws

> output := ModelReduce(RS, [E,C,P,S]): # only one output > output[1][1]; 2 d km1 C(t) k2 [-- E(t) = ------------------------------, dt 2 km1 S(t) + k1 S(t) + km1 C(t) 2 d km1 C(t) k2 d

• - C(t) = - ------------------------------, -- P(t) = k2 C(t),

dt 2 dt km1 S(t) + k1 S(t) + km1 C(t) d S(t) (km1 + k1 S(t)) k2 C(t)

• - S(t) = - ------------------------------]

dt 2 km1 S(t) + k1 S(t) + km1 C(t) # QSSA Assumptions: > output[1][2]; [k1 E(t) S(t) - km1 C(t)]

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 43 / 49

ModelReduce

with conservation laws

• utput := ModelReduce(RS, [E,C,P,S], useConservationLaws=true):

red_sys := output[1][1]: red_sys := subs(C_0=0, red_sys): red_sys[4]; d E_0 k2 k1 S(t) (km1 + k1 S(t))

• - S(t) = - ---------------------------------------------

dt 2 2 2 k1 S(t) + 2 S(t) km1 k1 + km1 k1 E_0 + km1

• ne has a diff. equation in S(t) only, but still many parameters.

we are far from: ˙ S(t) = − Vm S(t) K + S(t) with Vm = k2 E0 and K =

k−1 k1 (or k−1+k2 k1

)

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 44 / 49

Exact Reductions

Idea

• btain an equivalent system which is "easier" to study

Two routines

InvariantizeByScalings: reduce the number of parameters

reduce the complexity of parameter value exploration helps hand analysis (shown on the enzymatic degradation)

SemiRectifySteadyPoints: reduce the number of parameters on which the steady points depend

some parameters only affect the dynamic

Based on change of coordinates

• ne imposes monomial maps (xi → xα1

1

· · · xαn

n ). Ex: Vm → Vm/k1

ensures equivalence of systems, and positivity of parameters

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 45 / 49

InvariantizeByScalings

InvariantizeByScalings

Input: a dynamical system, list of positive variables, list of remaining variables. Output: a reduced system, a change of variables, list of parameters removed

> red_sys[4]; d E_0 k2 k1 S(t) (k1 S(t) + km1)

• - S(t) = - ---------------------------------------------

dt 2 2 2 k1 S(t) + 2 S(t) km1 k1 + km1 k1 E_0 + km1 > output := InvariantizeByScalings(red_sys, [k1,km1,k2], [], > fixedvars=[E,C,P,S], > scaletime=false): > red_sys2 := output[1]: > red_sys2[4]; d S(t) k2 E_0 (S(t) + km1)

• - S(t) = - -----------------------------------

dt 2 2 S(t) + 2 km1 S(t) + E_0 km1 + km1 > output[2], output[3]; km1 [km1 = ---], [k1] k1

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 46 / 49

InvariantizeByScalings (II)

A last hypothesis

If one assumes S >> E0, one retrieves a simpler formula (classical hypothesis)

> red_sys2[4]; d S(t) k2 E_0 (S(t) + km1)

• - S(t) = - -----------------------------------

dt 2 2 2 km1 S(t) + S(t) + E_0 km1 + km1 > ApproximateSmallQuantity(red_sys2[4], E_0<S(t)); d S(t) k2 E_0

• - S(t) = - -----------

dt S(t) + km1

Better simplification ?

Another step can automatically rename k2 E_0 into k2. We simplify by reducing the number of parameters. Further work: find changes of variables which reduce the size of the equations ( ˙ x(t) = 1 + a + a b x(t) into ˙ x(t) = 1 + a + b x(t))

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 47 / 49

Summary

1

Background Differential elimination Slow/fast dynamics Tikhonov theorem

2

Slow/fast chemical reaction systems reduction Michaelis Menten example General method

3

Application to physics Communicating vessels Diffusion Pendulum Other examples ...

4

MABSys

5

Conclusion

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 48 / 49

Conclusion

Main idea

many systems with constraints are limit case of systems with slow/fast timescales differential elimination can help to remove the fast timescale the dynamics of the fast phenomena is not important, only the equilibria of the fast phenomena is important,

Further work

convince people to use differential elimination for QSSA problems convince physicists/biologists to use packages for reduction (QSSA, parameter reduction, . . . ) write a C version of MABSys with a GPL licence (ModelReduce is done, parameter reduction is not, . . . ) → ANR LEDA

Lemaire (Lille I) Differential Algebra and QSSA DART 2010 49 / 49