SYSMETAB Non Stationary Metabolic flux analysis in isotope labeling - - PowerPoint PPT Presentation

▶
sysmetab
SMART_READER_LITE
LIVE PREVIEW

SYSMETAB Non Stationary Metabolic flux analysis in isotope labeling - - PowerPoint PPT Presentation

Aug 11, 2023 471 likes •792 views

SYSMETAB Non Stationary Metabolic flux analysis in isotope labeling experiments using the adjoint approach S. Mottelet 1 , G. Gaullier 2 Universit e de Technologie de Compi` egne 1 Laboratoire TIMR EA 4297 2 Laboratoire LMAC EA 2222 S.

slide-1
SLIDE 1

SYSMETAB

Non Stationary Metabolic flux analysis in isotope labeling experiments using the adjoint approach

  • S. Mottelet1, G. Gaullier2

Universit´ e de Technologie de Compi` egne

1Laboratoire TIMR EA 4297 2Laboratoire LMAC EA 2222

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 1 / 16

slide-2
SLIDE 2

The project

2013-2016

Post-doctoral fellows G. Sadaka (2013-2014), G. Gaullier (2015-2016) This work was performed, in partnership with the SAS PIVERT, within the frame of the French Institute for the Energy Transition (Institut pour la Transition Energ´ etique - ITE) P .I.V.E.R.T. (www.institut-pivert.com) selected as an Investment for the Future (”Investissements d’Avenir”). This work was supported, as part of the Investments for the Future, by the French Government under the reference ANR-001-01” Work Package 4 ”Metabolism of lipids from the plant to microorganisms” (sub-task 4.1.22 ”Metabolic Modelling” joint work with INRA-BIMLip, UMR 6022 GEC - UTC).

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 2 / 16

slide-3
SLIDE 3

Outline

A model problem and some important notions Non stationary MFA : motivation, software Sysmetab, the adjoint approach and other innovations Trends and conclusion

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 3 / 16

slide-4
SLIDE 4

A model problem

and some important notions

Continuous enzymatic reaction

Q

  • ! S

v

  • ! P

Q

  • !

S0 = v(S) + Q, t > 0 P0 = S0, S(0) = P(0) = 0.

Monod/Michaelis kinetics

v(S) = VmaxS K + S , Q < Vmax Metabolic stationary state lim

t!1 S(t) = S⇤

How to identify K et Vmax when the metabolic stationary state is reached and maintained ?

  • ! feed the reactor with a mix of two isotopomers of S = S1 + S2 :

Q1

  • ! S1

v1

  • ! P1

Q1

  • !

Q2

  • ! S2

v2

  • ! P2

Q2

  • !

Q1 = αQ, Q2 = (1 α)Q, S = S1 + S2, P = P1 + P2.

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 4 / 16

slide-5
SLIDE 5

A model problem

and some important notions

Continuous enzymatic reaction

Q

  • ! S

v

  • ! P

Q

  • !

S0 = v(S) + Q, t > 0 P0 = S0, S(0) = P(0) = 0.

Monod/Michaelis kinetics

v(S) = VmaxS K + S , Q < Vmax Metabolic stationary state lim

t!1 S(t) = S⇤

How to identify K et Vmax when the metabolic stationary state is reached and maintained ?

  • ! feed the reactor with a mix of two isotopomers of S = S1 + S2 :

Q1

  • ! S1

v1

  • ! P1

Q1

  • !

Q2

  • ! S2

v2

  • ! P2

Q2

  • !

Q1 = αQ, Q2 = (1 α)Q, S = S1 + S2, P = P1 + P2.

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 4 / 16

slide-6
SLIDE 6

A model problem

and some important notions

Continuous enzymatic reaction

Q

  • ! S

v

  • ! P

Q

  • !

S0 = v(S) + Q, t > 0 P0 = S0, S(0) = P(0) = 0.

Monod/Michaelis kinetics

v(S) = VmaxS K + S , Q < Vmax Metabolic stationary state lim

t!1 S(t) = S⇤

How to identify K et Vmax when the metabolic stationary state is reached and maintained ?

  • ! feed the reactor with a mix of two isotopomers of S = S1 + S2 :

Q1

  • ! S1

v1

  • ! P1

Q1

  • !

Q2

  • ! S2

v2

  • ! P2

Q2

  • !

Q1 = αQ, Q2 = (1 α)Q, S = S1 + S2, P = P1 + P2.

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 4 / 16

slide-7
SLIDE 7

A model problem

and some important notions

Continuous enzymatic reaction

Q

  • ! S

v

  • ! P

Q

  • !

S0 = v(S) + Q, t > 0 P0 = S0, S(0) = P(0) = 0.

Monod/Michaelis kinetics

v(S) = VmaxS K + S , Q < Vmax Metabolic stationary state lim

t!1 S(t) = S⇤

How to identify K et Vmax when the metabolic stationary state is reached and maintained ?

  • ! feed the reactor with a mix of two isotopomers of S = S1 + S2 :

Q1

  • ! S1

v1

  • ! P1

Q1

  • !

Q2

  • ! S2

v2

  • ! P2

Q2

  • !

Q1 = αQ, Q2 = (1 α)Q, S = S1 + S2, P = P1 + P2.

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 4 / 16

slide-8
SLIDE 8

A model problem

and some important notions

Isotopic labeling S1 = 13C6H12O6, S2 = 12C6H12O6 Under the hypothesis justifying the Monod/Michaelis kinetics S0 = Vmax S K + S + Q, S0

1 = Vmax

S1 K + S + αQ, S0

2 = Vmax

S2 K + S + (1 α)Q. Isotopic labeling experiment

I t < Tm, α = 0 until the metabolic stationary state is reached I t ≥ Tm, α = 1 until the isotopic stationary state is reached

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 5 / 16

slide-9
SLIDE 9

A model problem

and some important notions

Isotopic labeling S1 = 13C6H12O6, S2 = 12C6H12O6 Under the hypothesis justifying the Monod/Michaelis kinetics S0 = Vmax S K + S + Q, S0

1 = Vmax

S1 K + S + αQ, S0

2 = Vmax

S2 K + S + (1 α)Q. Isotopic labeling experiment

I t < Tm, α = 0 until the metabolic stationary state is reached I t ≥ Tm, α = 1 until the isotopic stationary state is reached

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 5 / 16

slide-10
SLIDE 10

A model problem

and some important notions

Isotopic labeling S1 = 13C6H12O6, S2 = 12C6H12O6 Under the hypothesis justifying the Monod/Michaelis kinetics S0 = Vmax S K + S + Q, S0

1 = Vmax

S1 K + S + αQ, S0

2 = Vmax

S2 K + S + (1 α)Q. Isotopic labeling experiment

I t < Tm, α = 0 until the metabolic stationary state is reached I t ≥ Tm, α = 1 until the isotopic stationary state is reached

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 5 / 16

slide-11
SLIDE 11

A model problem

and some important notions

For t 50 a first order kinetics is observed for S1 and S2 S0

1 = µS1 + αQ,

S0

2 = µS2 + (1 α)Q,

µ = Vmax K + S⇤, S⇤ = KQ Vmax Q µ can be estimated from S1 and/or S2 for t 2 [50, 60]

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 6 / 16

slide-12
SLIDE 12

A model problem

and some important notions

For t 50 a first order kinetics is observed for S1 and S2 S0

1 = µS1 + αQ,

S0

2 = µS2 + (1 α)Q,

µ = Vmax K + S⇤, S⇤ = KQ Vmax Q µ can be estimated from S1 and/or S2 for t 2 [50, 60]

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 6 / 16

slide-13
SLIDE 13

A model problem

Actual state variables in MFA

Isotopomer fractions : si = Si S Complete kinetics (metabolic/isotopic) S0 = v(S) + Q, (Ss1)0 = v(S)s1 + αQ, (Ss2)0 = v(S)s2 + (1 α)Q. the reaction rate is shared proportionally to isotopomer fractions Sole isotopic kinetics S⇤s0

1 = v(S⇤)s1 + αQ,

S⇤s0

2 = v(S⇤)s2 + (1 α)Q.

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 7 / 16

slide-14
SLIDE 14

A model problem

Actual state variables in MFA

Isotopomer fractions : si = Si S Complete kinetics (metabolic/isotopic) S0 = v(S) + Q, (Ss1)0 = v(S)s1 + αQ, (Ss2)0 = v(S)s2 + (1 α)Q. the reaction rate is shared proportionally to isotopomer fractions Sole isotopic kinetics S⇤s0

1 = v(S⇤)s1 + αQ,

S⇤s0

2 = v(S⇤)s2 + (1 α)Q.

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 7 / 16

slide-15
SLIDE 15

A model problem

Actual state variables in MFA

Isotopomer fractions : si = Si S Complete kinetics (metabolic/isotopic) S0 = v(S) + Q, (Ss1)0 = v(S)s1 + αQ, (Ss2)0 = v(S)s2 + (1 α)Q. the reaction rate is shared proportionally to isotopomer fractions Sole isotopic kinetics S⇤s0

1 = v(S⇤)s1 + αQ,

S⇤s0

2 = v(S⇤)s2 + (1 α)Q.

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 7 / 16

slide-16
SLIDE 16

Non stationary MFA

Motivation

Why make some experiments once the metabolic state is reached ?

I The particular form of the kinetics (and associated parameters) disappears

and is replaced by asymptotic reaction rates, called fluxes

I Remaining unknowns : fluxes and asymptotic concentration of metabolites

called pool sizes

I If the reaction network is exhaustive, the model is reliable

Why make some experiments before the isotopic stationary state is reached ?

I Cost of the labeled substrate I Time constants can be very large for some organism I Qualitative and quantitative gain of information (narrower confidence

intervals for fluxes)

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 8 / 16

slide-17
SLIDE 17

Non stationary MFA

Software

INCA/MFA suite (Young et al., Vanderbilt Univ., Nashville, 10/2013)

I BFGS inspired, Matlab implementation, only MS and H1-RMN

measurements, closed source

Sysmetab (Mottelet et al., UTC, 1/2016)

I BFGS inspired / feasible directions (FSQP), Scilab implementation, adjoint

approach, open-source

influx si (Sokol et al., LISBP Toulouse, 4/2016)

I Gauss-Newton inspired / feasible directions (NLSIC), R/C++ implementation,

  • pen-source
  • S. Mottelet et al., Metabolic Flux Analysis in Isotope Labeling Experiments using the Adjoint

Approach, IEEE/ACM Transactions on Computational Biology and Bioinformatic, 2016 (to appear, DOI : 10.1109/TCBB.2016.2244299)

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 9 / 16

slide-18
SLIDE 18

Non stationary MFA

The identification problem in continuous time

State equation for t 2 [0, T] D(s⇤)x0(t) = f(x(t), v), x(t) 2 RN is the vector of cumomers x1(t), . . . , xn(t) Observation at times τi, i = 1 . . . m y(τi) = C0 + Cx(t) = C0 +

n

X

k=1

Ckxk(τi), i = 1 . . . m Cost function (where zi is the measurement at time τi) J(v, s⇤) =

m

X

i=1

ky(τi) zik2 Identification problem : find v, s⇤ (p parameters) minimizing J(v, s⇤)

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 10 / 16

slide-19
SLIDE 19

Non stationary MFA

The identification problem in continuous time

State equation for t 2 [0, T] D(s⇤)x0(t) = f(x(t), v), x(t) 2 RN is the vector of cumomers x1(t), . . . , xn(t) Observation at times τi, i = 1 . . . m y(τi) = C0 + Cx(t) = C0 +

n

X

k=1

Ckxk(τi), i = 1 . . . m Cost function (where zi is the measurement at time τi) J(v, s⇤) =

m

X

i=1

ky(τi) zik2 Identification problem : find v, s⇤ (p parameters) minimizing J(v, s⇤)

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 10 / 16

slide-20
SLIDE 20

Non stationary MFA

The identification problem in continuous time

State equation for t 2 [0, T] D(s⇤)x0(t) = f(x(t), v), x(t) 2 RN is the vector of cumomers x1(t), . . . , xn(t) Observation at times τi, i = 1 . . . m y(τi) = C0 + Cx(t) = C0 +

n

X

k=1

Ckxk(τi), i = 1 . . . m Cost function (where zi is the measurement at time τi) J(v, s⇤) =

m

X

i=1

ky(τi) zik2 Identification problem : find v, s⇤ (p parameters) minimizing J(v, s⇤)

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 10 / 16

slide-21
SLIDE 21

Non stationary MFA

The identification problem in continuous time

State equation for t 2 [0, T] D(s⇤)x0(t) = f(x(t), v), x(t) 2 RN is the vector of cumomers x1(t), . . . , xn(t) Observation at times τi, i = 1 . . . m y(τi) = C0 + Cx(t) = C0 +

n

X

k=1

Ckxk(τi), i = 1 . . . m Cost function (where zi is the measurement at time τi) J(v, s⇤) =

m

X

i=1

ky(τi) zik2 Identification problem : find v, s⇤ (p parameters) minimizing J(v, s⇤)

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 10 / 16

slide-22
SLIDE 22

Sysmetab

Innovations

Computation of the gradient of J(v, s⇤) (for BFGS inspired methods)

I Sysmetab uses the adjoint method : compute the adjoint state λ(t)

x(t) ∈ RN, λ(t) ∈ RN : ODE of size N + N

I others use the direct method : compute the sensitivity matrix ∂vx(t)

x(t) ∈ RN, ∂vx(t) ∈ MN,p(R) : ODE of size N + Np

Practical implementation in Sysmetab

I Time discretization + SDIRK ODE schemes (Singly Diagonally Implicit

Runge-Kutta) of order 1, 2 and 4 :

D(s⇤)(xk+1 xk) = Φ(xk, v), xk ⇡ x(tk) How to choose order and stepsize tk+1 tk = h = T

nT of the scheme ?

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 11 / 16

slide-23
SLIDE 23

Sysmetab

Innovations

Computation of the gradient of J(v, s⇤) (for BFGS inspired methods)

I Sysmetab uses the adjoint method : compute the adjoint state λ(t)

x(t) ∈ RN, λ(t) ∈ RN : ODE of size N + N

I others use the direct method : compute the sensitivity matrix ∂vx(t)

x(t) ∈ RN, ∂vx(t) ∈ MN,p(R) : ODE of size N + Np

Practical implementation in Sysmetab

I Time discretization + SDIRK ODE schemes (Singly Diagonally Implicit

Runge-Kutta) of order 1, 2 and 4 :

D(s⇤)(xk+1 xk) = Φ(xk, v), xk ⇡ x(tk) How to choose order and stepsize tk+1 tk = h = T

nT of the scheme ?

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 11 / 16

slide-24
SLIDE 24

Sysmetab

Innovations

Real error (identified parameters - true parameters, central metabolism of E. Coli, 10 measurement points)

10 20 30 40 50 60 10

  • 6

10

  • 5

10

  • 4

10

  • 3

10

  • 2

10

  • 1

10 10 Error (est. param. vs true values) sdirk4b (slope=-4.29) sdirk2a (slope=-1.91) Euler (slope=-0.68) Comparison of schemes nT (size of time grid)

The error on identified parameters inherits the ODE scheme order : this allows a smart error estimation method !

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 12 / 16

slide-25
SLIDE 25

Sysmetab

Innovations

Error estimation by the Richardson method ˆ Θ(h) is the vector of estimated parameters with ODE scheme of order q and step h (α = 1.5, for example) ˆ Θ(h) = Θ + Chq + O(hq+1), ˆ Θ(αh) = Θ + C(αh)q + O(hq+1), We obtain the estimate C ⇡ ˆ Θ(h) ˆ Θ(αh) (1 αq)hq , which gives ˆ Θ(h) Θ ⇡ ˆ Θ(h) ˆ Θ(αh) (1 αq) This estimation has to be considered relatively to Θ estimated standard deviation

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 13 / 16

slide-26
SLIDE 26

Sysmetab

Results for the central metabolism of E. Coli PGA FruBP Suc

M0 M1 M2 M4 M5 M6 M3

Glc6P Fru6P ICit PEP Time (s) Time (s) Time (s) Time (s)

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 14 / 16

slide-27
SLIDE 27

Sysmetab

Results for the central metabolism of E. Coli

nT = 10, SDIRK4 scheme (flux and error estimation take 40 seconds) Data Nonstationary Stationary

Free fluxes

Value

  • st. dev.
  • rel. err. (%)

Value

  • std. dev.

Glucupt 1.n 0.808 0.0019 5.2 0.808 0.0037 pyk.n 1.537 0.1817 0.1 1.548 0.1336 zwf.n 0.157 0.0212 2.5 0.129 0.0323 gnd.n 0.154 0.0105 4.2 0.129 0.0278

  • ut Ac.n

0.213 1e-5 1e-7 0.213 1e-5 ald.x 0.670 0.0427 2.8 0.842 0.1010 eno.x 999 1 999 1 ta.x 0.548 0.0203 0.7 0.553 0.0387 tk1.x 0.189 0.0216 4.7 0.164 0.0502 tk2.x 0.0120 0.0278 fum a.x 0.106 0.1846 0.05 0.184 ppc.x 0.224 0.0948 3.6 0.169 0.0798

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 15 / 16

slide-28
SLIDE 28

Sysmetab

Results for the central metabolism of E. Coli

nT = 10, SDIRK4 scheme (flux and error estimation take 40 seconds) Data Nonstationary Stationary

Free fluxes

Value

  • st. dev.
  • rel. err. (%)

Value

  • std. dev.

Glucupt 1.n 0.808 0.0019 5.2 0.808 0.0037 pyk.n 1.537 0.1817 0.1 1.548 0.1336 zwf.n 0.157 0.0212 2.5 0.129 0.0323 gnd.n 0.154 0.0105 4.2 0.129 0.0278

  • ut Ac.n

0.213 1e-5 1e-7 0.213 1e-5 ald.x 0.670 0.0427 2.8 0.842 0.1010 eno.x 999 1 999 1 ta.x 0.548 0.0203 0.7 0.553 0.0387 tk1.x 0.189 0.0216 4.7 0.164 0.0502 tk2.x 0.0120 0.0278 fum a.x 0.106 0.1846 0.05 0.184 ppc.x 0.224 0.0948 3.6 0.169 0.0798

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 15 / 16

slide-29
SLIDE 29

Sysmetab

Results for the central metabolism of E. Coli

nT = 10, SDIRK4 scheme (flux and error estimation take 40 seconds) Data Nonstationary Stationary

Free fluxes

Value

  • st. dev.
  • rel. err. (%)

Value

  • std. dev.

Glucupt 1.n 0.808 0.0019 5.2 0.808 0.0037 pyk.n 1.537 0.1817 0.1 1.548 0.1336 zwf.n 0.157 0.0212 2.5 0.129 0.0323 gnd.n 0.154 0.0105 4.2 0.129 0.0278

  • ut Ac.n

0.213 1e-5 1e-7 0.213 1e-5 ald.x 0.670 0.0427 2.8 0.842 0.1010 eno.x 999 1 999 1 ta.x 0.548 0.0203 0.7 0.553 0.0387 tk1.x 0.189 0.0216 4.7 0.164 0.0502 tk2.x 0.0120 0.0278 fum a.x 0.106 0.1846 0.05 0.184 ppc.x 0.224 0.0948 3.6 0.169 0.0798

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 15 / 16

slide-30
SLIDE 30

Trends and conclusion

Sysmetab uses up to date numerical methods and can compete with other software (stationary and non stationary MFA)

I Adjoint computation of the gradient −

→ fast identification

I High order implicit ODE scheme + error estimation I Handling of compartments and dilution effects

Ongoing developpements

I Speed improvements (switch from Scilab to Julia) I Documentation (...) I CellDesigner plug-in

http://forge.scilab.org/p/sysmetab

  • S. Mottelet, G. Gaullier (UTC)

SYSMETAB 16 / 16