Model order reduction of biochemical networks by tropical geometry - - PowerPoint PPT Presentation

Model order reduction of biochemical networks by tropical geometry methods

Ovidiu Radulescu University of Montpellier

Bonn, March 2018, 2 / 27

Models of intracellular biochemistry

Metabolic networks: synthesis, growth Gene networks: Decision making Signalling networks: information transfer

Bonn, March 2018, 3 / 27

Models of intracellular biochemistry: signalling

MAPK signalling, Huang and Ferrell 96 Biochemical reactions: MAPKKK + E1 → MAPKKK : E1 → MAPKKK ∗ + E1, MAPKKK ∗ + E2 → MAPKKK ∗ : E2 → MAPKKK + E2,...

Bonn, March 2018, 4 / 27

Chemical reactions networks (CRN)

◮ x ∈ Rn

+ concentrations vector, xi concentration of the

species i.

Bonn, March 2018, 4 / 27

Chemical reactions networks (CRN)

◮ x ∈ Rn

+ concentrations vector, xi concentration of the

species i.

◮ X ⇄ X +νj reactions νj ∈ Zn, stoichiometric vector

Bonn, March 2018, 4 / 27

Chemical reactions networks (CRN)

◮ x ∈ Rn

+ concentrations vector, xi concentration of the

species i.

◮ X ⇄ X +νj reactions νj ∈ Zn, stoichiometric vector ◮ reaction rates R+/− j

(x) (reactions / unit volume / unit time)

Bonn, March 2018, 4 / 27

Chemical reactions networks (CRN)

◮ x ∈ Rn

+ concentrations vector, xi concentration of the

species i.

◮ X ⇄ X +νj reactions νj ∈ Zn, stoichiometric vector ◮ reaction rates R+/− j

(x) (reactions / unit volume / unit time)

◮ kinetic law, for instance mass action

R+

j (x) = k+ j ∏i∈reactants(j) x

ν+

ji

i

, R−

j (x) = k− j ∏i∈products(j) x

ν−

ji

i

Bonn, March 2018, 4 / 27

Chemical reactions networks (CRN)

◮ x ∈ Rn

+ concentrations vector, xi concentration of the

species i.

◮ X ⇄ X +νj reactions νj ∈ Zn, stoichiometric vector ◮ reaction rates R+/− j

(x) (reactions / unit volume / unit time)

◮ kinetic law, for instance mass action

R+

j (x) = k+ j ∏i∈reactants(j) x

ν+

ji

i

, R−

j (x) = k− j ∏i∈products(j) x

ν−

ji

i ◮ Deterministic dynamics dx dt = ∑r j=1 νj(R+ j (x)− R− j (x)).

For mass action, these are polynomial ODEs

Bonn, March 2018, 5 / 27

Model order reduction produce models with less variables and differential equations; popular method : aggregation of species and reactions

Bonn, March 2018, 6 / 27

Multiple separated scales Widely distributed concentrations, in log scale Widely distributed timescales, in log scale

Produced with the model in Radulescu et al BMC Systems Biol. 2008

Bonn, March 2018, 7 / 27

Slow/fast systems

from Chiavazzo et al Comm.Comp.Phys. 2007 from Hung and Sheperd 12th Ann. Int. Detonation Symp. 2002

In slow/fast systems, fast variables relax quickly, then they are slaved : the system moves on a low dimensional invariant manifold describing the reduced model. The reduced model can change with time.

Bonn, March 2018, 8 / 27

Slow/fast systems : Tikhonov theorem dx dt = 1

ηf(x,y)fast dynamics

dy dt = g(x,y)slow dynamics

η is the fast/slow timescale ratio.

Tikhonov : If for any y the fast dynamics has a hyperbolic point attractor, then after a quick transition the system evolves according to: dy dt = g(x,y) and f(x,y)=0 fast variables are slaved by slow ones

Bonn, March 2018, 9 / 27

Model reduction Several steps:

Bonn, March 2018, 9 / 27

Model reduction Several steps:

◮ Determine the slow/fast decomposition (who are the small

parameter η, the slow and the fast variables?): Jacobian based numerical methods (CSP , ILDM, COPASI implementation); tropical geometry based methods.

Bonn, March 2018, 9 / 27

Model reduction Several steps:

◮ Determine the slow/fast decomposition (who are the small

parameter η, the slow and the fast variables?): Jacobian based numerical methods (CSP , ILDM, COPASI implementation); tropical geometry based methods.

◮ Solve f(x,y) = 0 for x : hard, few symbolic

methods(Grobner bases, sparse polynomial systems (Grigoriev, Weber 2012)); use tropical approximations?

Bonn, March 2018, 9 / 27

Model reduction Several steps:

◮ Determine the slow/fast decomposition (who are the small

parameter η, the slow and the fast variables?): Jacobian based numerical methods (CSP , ILDM, COPASI implementation); tropical geometry based methods.

◮ Solve f(x,y) = 0 for x : hard, few symbolic

methods(Grobner bases, sparse polynomial systems (Grigoriev, Weber 2012)); use tropical approximations?

◮ Pool reactions, prune species : use elementary modes

(Radulescu et al 2008).

Bonn, March 2018, 9 / 27

Model reduction Several steps:

◮ Determine the slow/fast decomposition (who are the small

parameter η, the slow and the fast variables?): Jacobian based numerical methods (CSP , ILDM, COPASI implementation); tropical geometry based methods.

◮ Solve f(x,y) = 0 for x : hard, few symbolic

methods(Grobner bases, sparse polynomial systems (Grigoriev, Weber 2012)); use tropical approximations?

◮ Pool reactions, prune species : use elementary modes

(Radulescu et al 2008).

◮ Pool species, prune reactions : use conservation laws

(Gorban, Radulescu, Zinovyev 2010).

Bonn, March 2018, 10 / 27

Scaling and tropical geometry Polynomial kinetics (mass action law) dxi dt = Pi(x) =

Mi

∑

j=1

kijxαij, xαij = x

α1

ij

1 ...x

αn

ij

n , 1 ≤ i ≤ n.

αij ∈ Nn are multi-indices, kij > 0 are kinetic parameters.

xi are positive species concentrations. D = {xi ≥ 0,1 ≤ i ≤ n} is a positive invariant set of the ODEs.

Bonn, March 2018, 10 / 27

Scaling and tropical geometry Polynomial kinetics (mass action law) dxi dt = Pi(x) =

Mi

∑

j=1

kijxαij, xαij = x

α1

ij

1 ...x

αn

ij

n , 1 ≤ i ≤ n.

αij ∈ Nn are multi-indices, kij > 0 are kinetic parameters.

xi are positive species concentrations. D = {xi ≥ 0,1 ≤ i ≤ n} is a positive invariant set of the ODEs. Consider ε > 0 a small scaling parameter.

Bonn, March 2018, 10 / 27

Scaling and tropical geometry Polynomial kinetics (mass action law) dxi dt = Pi(x) =

Mi

∑

j=1

kijxαij, xαij = x

α1

ij

1 ...x

αn

ij

n , 1 ≤ i ≤ n.

αij ∈ Nn are multi-indices, kij > 0 are kinetic parameters.

xi are positive species concentrations. D = {xi ≥ 0,1 ≤ i ≤ n} is a positive invariant set of the ODEs. Consider ε > 0 a small scaling parameter. Rescale kij = ¯ kijεγij with ¯ kij ∼ 1, for instance

γij = round(d log(kij)/log(ε))/d, d ∈ N∗.

Bonn, March 2018, 10 / 27

Scaling and tropical geometry Polynomial kinetics (mass action law) dxi dt = Pi(x) =

Mi

∑

j=1

kijxαij, xαij = x

α1

ij

1 ...x

αn

ij

n , 1 ≤ i ≤ n.

αij ∈ Nn are multi-indices, kij > 0 are kinetic parameters.

xi are positive species concentrations. D = {xi ≥ 0,1 ≤ i ≤ n} is a positive invariant set of the ODEs. Consider ε > 0 a small scaling parameter. Rescale kij = ¯ kijεγij with ¯ kij ∼ 1, for instance

γij = round(d log(kij)/log(ε))/d, d ∈ N∗.

Suppose the same is possible for concentrations: xi = εai ¯

• xi. Not

possible if xi are zero or get exponentially close to zero. Some persistence condition is needed:

∃ai,Ci > 0, xi(t) > Ciεai,∀t ≥ 0.

Bonn, March 2018, 11 / 27

Scaling and tropical geometry Rescaled equations

εai d¯

xi dt =

Mi

∑

j=1

εµij(a)¯

kij¯ xαij, 1 ≤ i ≤ n. where

µij(a) = γij +

n

∑

l=1

αij

l al

Bonn, March 2018, 12 / 27

Crazy-quilt: the big picture of nonlinear, multiscale dynamics.

Abstract representation of metastability as itinerant trajectory in a patchy phase space landscape.

Bonn, March 2018, 13 / 27

Tropical equilibration problem

◮ Find a such that there is (j,j′), j = j′ such that

i) µij(a) = µij′(a), ii) µij(a) ≤ µil(a) for all l = j,j′, iii) kijkij′ < 0. at least one negative and one positive dominant monomials have the same orders

Bonn, March 2018, 13 / 27

Tropical equilibration problem

◮ Find a such that there is (j,j′), j = j′ such that

i) µij(a) = µij′(a), ii) µij(a) ≤ µil(a) for all l = j,j′, iii) kijkij′ < 0. at least one negative and one positive dominant monomials have the same orders

◮ Tropically truncated system (describes the dominant, fast

dynamics) d¯ xi dt = ενi(a)(|¯ kij|¯ xαij −|¯ kij′|¯ xαij′

+ terms of same order) νi(a) = µi(a)− ai gives the timescale of xi : small νi

means fast xi.

Bonn, March 2018, 14 / 27

Tropical equilibration problem Equivalent formulation min

kij>0(µij(a)) = min kij<0(µij(a)), 1 ≤ i ≤ n

µij(a) = γij +

n

∑

l=1

αij

l al

system of polynomial equations in min-plus algebra a⊕ b = min(a,b) a⊗ b = a+ b

Bonn, March 2018, 15 / 27

Tropical equilibrations : geometrical interpretation

◮ polynomial

−k1x1 + k2x1x2 + k3x2

◮ orders of magnitudes

(valuations) x1 = εa1, x2 = εa2, k1 = εγ1, k2 = εγ2, k3 = εγ3.

◮ k1x1 = εγ1+a1,

k2x1x2 = εγ2+a1+a2, k3x2 = εγ3+a2.

◮ tropical curve : Min(γ1 +

a1,γ2 + a1 + a2,γ3 + a2) is attained at least twice

◮ tropical equilibration :

Min(γ2 +a1 +a2,γ3 +a2) =

γ1 + a1

logx1 logx2 k2x1x2 dominant k1x1 dominant k3x2 dominant k1x1 = k2x1x2 k1x1 = k3x2 k2x1x2 = k3x2

Bonn, March 2018, 16 / 27

Tropical equilibrations : quasi-equilibrium

dx1 dt

= −k1x1 + k2x1x2 + k3x2

dx2 dt

=

k1x1 − k2x1x2 −(k3 + k′

3)x2

with k′

3 << k3, in other words

γ′

3 > γ3 logx1 logx2 k2x1x2 dominant k1x1 dominant k3x2 dominant k1x1 = k2x1x2 k1x1 = k3x2 k2x1x2 = k3x2

We show (Samal et al BMB 2015) that x1 and x2 are fast and x1 + x2 is slow

◮ k1x1 = k2x2 =

⇒ x2 = k1x1/k3 (linear regime)

◮ k1x1 = k2x1x2 =

⇒ x2 = k1/k2 (saturation)

◮ −k1x1 + k2x1x2 + k3x2 = 0 =

⇒ x2 = k1x1/(k3 + k2x1) (two branches

reduction)

Bonn, March 2018, 17 / 27

Partial tropical equilibrations : quasi-steady state

dx1 dt

= −k1x1 + k2x1x2 + k3x2

dx2 dt

=

k1x1 − k2x1x2 −(k3 + k′

3)x2

with k′

3 >> k3, in other words

γ′

3 < γ3 logx1 logx2 k1x1 = k2x1x2 k1x1 = k3x2 k1x1 = k′

3x 2

k2x1x2 = k3x2 k2x1x2 = k′

3x2

We show (Samal et al BMB 2015) that if a1 < min(γ3 −γ2,γ1 −γ2) or

γ3 −γ2 < min(a1,γ1 −γ2), then x2 is fast and x1 is slow

◮ k1x1 = k′

3x2 =

⇒ x2 = k1x1/k′

3 (linear regime)

◮ k1x1 = k2x1x2 =

⇒ x2 = k1/k2 (saturation)

◮ −k1x1 + k2x1x2 + k′

3x2 = 0 =

⇒ x2 = k1x1/(k′

3 + k2x1) (two branches

reduction)

Bonn, March 2018, 18 / 27

Model reduction and tropical geometry What can we do with tropical methods?

◮ Compute timescales, identify fast and slow variables. ◮ Obtain simplified, truncated algebraic equations for the

elimination of fast variables.

◮ One branch of tropical solutions → one reduced model,

but several branches can be combined to extend the validity of the approximation.

◮ Are the tropical equilibration + separated timescales

conditions sufficient for accurate reduction?

Bonn, March 2018, 19 / 27

Rigorous asymptotic results

Suppose we have the formal timescales of variables, including identification

• f fast cycles (combinations of fast variables conserved by the fast dynamics).

Bonn, March 2018, 19 / 27

Rigorous asymptotic results

Suppose we have the formal timescales of variables, including identification

• f fast cycles (combinations of fast variables conserved by the fast dynamics).

Introduce the slow variables yk = ∑l cklxl, where k ∈ [1,r].

Bonn, March 2018, 19 / 27

Rigorous asymptotic results

Suppose we have the formal timescales of variables, including identification

• f fast cycles (combinations of fast variables conserved by the fast dynamics).

Introduce the slow variables yk = ∑l cklxl, where k ∈ [1,r]. We obtain the truncated system d¯ xi dt = ενi Fi(Xr, ¯

y,Xs,ε)

i ∈ [1,f − r],fast species dynamics d¯ xi dt = ενi Si(Xr, ¯

y,Xs,ε),

i ∈ [f + 1,n],slow species dynamics d¯ yk dt = ερk Yk(Xr, ¯

y,Xs,ε),

k ∈ [1,r],slow dynamics of cycle masses where ν1 ≤ ν2 ≤ ··· ≤ νn

Bonn, March 2018, 20 / 27

Rigorous asymptotic results: conditions

Assume the following conditions hold

◮ Gap conditions νf+1 > νf−r,ρk > νf−r. ◮ Uniform hyperbolicity condition For all ¯

y and Xs the fast truncated

system has a stable hyperbolic steady state

¯

xi = φi(¯

y,Xs)+ higher order terms,

i ∈ [1,f − r], such that d > C0εκ where d is the distance from the spectrum of linearisation of the fast truncated system to the imaginary axis; κ ≥ 0 is small enough and C0 is independent on ε.

Bonn, March 2018, 21 / 27

Rigorous asymptotic results: main theorem

Then for sufficiently small ε > 0 the full system has a locally attracting and locally invariant normally hyperbolic Cp (p > 1) smooth manifold defined by

¯

xi = φi(¯

y,Xs,ε)+ higher order terms,

i ∈ [1,f − r], (1) and the dynamics of the slow variables ¯

y,¯

xf+1,¯ xf+2,...,¯ xn for large times takes the form d¯ xi dt = ενi Si(φ, ¯

y,Xs,ε)+ higher order terms,

i ∈ [f + 1,n] d¯ yk dt = ερk Yk(φ, ¯

y,Xs,ε)+ higher order terms,

k ∈ [1,r] (2) (Radulescu, Vakulenko, Grigoriev MMNP 2015)

Bonn, March 2018, 22 / 27

A medium size example

TGFβ signaling model (important in cancer research).

Bonn, March 2018, 23 / 27

A medium size example

B1 B2 B3 B4 B5 B6 B7 B8 B9

Minimal tropical branches and geometric/algebraic connections.

Bonn, March 2018, 24 / 27

A medium size example Numerically computed trajectories and distances to minimal branches.

Bonn, March 2018, 25 / 27

Using tropical equilibrations for model reduction of the TGFβ model

Full model trajectories (blue) from 21 ODEs. Reduced model trajectories (red) from 8 ODEs. The reduction is based on tropical branches B2, B1.

Bonn, March 2018, 26 / 27

Conclusion

◮ Tropical methods solve crucial steps of chemical reaction

networks order reduction: slow/fast decomposition, elimination of fast variables.

Bonn, March 2018, 26 / 27

Conclusion

◮ Tropical methods solve crucial steps of chemical reaction

networks order reduction: slow/fast decomposition, elimination of fast variables.

◮ We have implemented algorithms for the computation of

tropical equilibrations and of timescales (talk by Chris Lueders).

Bonn, March 2018, 26 / 27

Conclusion

◮ Tropical methods solve crucial steps of chemical reaction

networks order reduction: slow/fast decomposition, elimination of fast variables.

◮ We have implemented algorithms for the computation of

tropical equilibrations and of timescales (talk by Chris Lueders).

◮ The slow/fast decomposition has to be fully formalized and

implemented; lacks an algorithm for detecting fast cycles; lacks a methodology for partial equilibrations;

Bonn, March 2018, 26 / 27

Conclusion

◮ Tropical methods solve crucial steps of chemical reaction

networks order reduction: slow/fast decomposition, elimination of fast variables.

◮ We have implemented algorithms for the computation of

tropical equilibrations and of timescales (talk by Chris Lueders).

◮ The slow/fast decomposition has to be fully formalized and

implemented; lacks an algorithm for detecting fast cycles; lacks a methodology for partial equilibrations;

◮ We need algorithms for the elimination of fast variables

(solving the truncated system of algebraic equations). One branch or several branches?

Bonn, March 2018, 27 / 27

Partial tropical equilibrations, extra conditions

Tropical curves and dominant vector field for Tyson 91 cell cycle model