Princeton

University

Equation-Free Uncertainty Quantification: An Application to Yeast Glycolytic Oscillations

Katherine A. Bold, Yu Zou, Ioannis G. Kevrekidis Department of Chemical Engineering and PACM Princeton University Michael A. Henson Department of Chemical Engineering University of Massachusetts, Amherst WCCM VII, LA July 16-22, 2006

Department of Chemical Engineering and PACM

Princeton

University Outline 1.

2.

3.

4.

5.

6.

Background for Uncertainty Quantification Fundamentals of Polynomial Chaos Stochastic Galerkin Method Equation-Free Uncertainty Quantification Application to Yeast Glycolytic Oscillations Remarks

Department of Chemical Engineering and PACM

Princeton

University Background for Uncertainty Quantification Uncertain Phenomena in science and engineering * Inherent Uncertainty: Uncertainty Principle of quantum mechanics, Kinetic theory of gas, … * Uncertainty due to lack of knowledge: randomness of BC, IC and parameters in a mathematical model, measurement errors associated with an inaccurate instrument, … Scopes of application * Estimate and predict propagation of probabilities for model variables: chemical reactants, biological oscillators, stock and bond values, structural random vibration,… * Design and decision making in risk management: optimal selection of parameters in a manufacturing process, assessment of an investment to achieve maximum profit,...

* Evaluate and update model predictions via experimental data: validate accuracy of a stochastic model based on experiment, data assimilation, … Modeling Techniques * Sampling methods (non-intrusive): Monte Carlo sampling, Quasi Monte Carlo, Latin Hypercube Sampling, Quadrature/Cubature rules * Non-sampling methods (intrusive) : perturbation methods, higher-order moment analysis, stochastic Galerkin method

Department of Chemical Engineering and PACM

Princeton

University Fundamentals of Polynomial Chaos The functional of independent random variables f (

ξ

),

ξ

 (  1 (  ),  2 (  ),  , 

n

(  )) can be used to represent a random variable, a random field or process.

Spectral expansion (Ghanem and Spanos, 1991)  (

ξ

) f (

ξ

is the probability measure of

ξ

)  .

j

   0

a j



j

if i≠j. The inner product <·, ·> is defined as (

ξ

 )

a

j ’s are PC coefficients, Ψ

j

’s are orthogonal polynomial functions with <Ψ

i

,Ψ

j

>=0

f

(

ξ

),

g

(

ξ

)   

f

(

ξ

)

g

(

ξ

)

d

 (

ξ

) , Notes Selection of Ψ j is dependent on the probability measure or distribution of e.g., (Xiu and Karniadarkis, 2002) if  is a Gaussian measure, then Ψ j if  (

ξ

(

ξ

) ) is a Lesbeque measure, then Ψ j are Hermite polynomials; are Legendre polynomials.

ξ

,

Department of Chemical Engineering and PACM

Princeton

University Stochastic Galerkin (PC expansion) Method (Ghanem and Spanos, 1991) Preliminary Formulation Input: random IC, BC, parameters Model Response: Solution * Model: e.g., ODE 

(

x

;

K

)

 .

x

 g

(

x

;

K

)



0

* Represent the input in terms of expansion of independent r.v.’s (KL, SVD, POD): e.g., time-dependent parameter

K n



m

 1  

m m



m

* Represent the response in terms of the truncated PC expansion

x

(

ξ

,

t

) 

j P

  0 

j

(

t

) 

j

(

ξ

) * The solution process involves solving for the PC coefficients

α

j (t), j=1,2,…,P

Department of Chemical Engineering and PACM

Princeton

University Stochastic Galerkin Method Solution technique: Galerkin projection  

(

j P

  0 

j

(

t

)



j

(

ξ

) ;

ξ

),



i

(

ξ

)



0

resulting in coupled ODE’s for

α j

(t), where .

A (

t

)  G ( A (

t

))  0 

i

A

(

t

)



(

 0

(

t

),

 1

(

t

),...,



P

(

t

))

T

Advantages and weakness * PC expansion has exponential convergence rate * Model reduction * Free of moment closure problems ? The coupled ODE’s of PC coefficients may not be obtained explicitly

Department of Chemical Engineering and PACM

Princeton

University Equation-free Uncertainty Quantification Coarse time-stepper (Kevrekidis

et al

., 2003, 2004) * Lifting (MC, quadrature/cubature): * Microsimulation:

x

(

ξ

i

,

t

0 ) 

j P

  0 

j

(

t

0 ) 

j

(

ξ

i

) ,

i

 1 , 2 ,  ,

N e

* Restriction: .

x

(

ξ

i

 g

x i t K i



0 ,

i



1, 2, ,

N e



j

(

t

)



x

(

ξ

,

t

),



j

(

ξ

)



/

 

j

(

ξ

),



j

(

ξ

)



,



x

(

ξ

,

t

), 

j

(

ξ

) 

N e i

   1

i x

(

ξ

i

,

t

) 

j

(

ξ

i

) For Monte Carlo sampling, 

i



1 /

N e

 with each sampling point.

j



0 , 1 ,



,

P

Department of Chemical Engineering and PACM

Princeton

University Projective Integration (Kevrekidis

et al

., 2003, 2004) Fixed-point Computation (Kevrekidis

et al

., 2003, 2004)   Lifting Restriction

Department of Chemical Engineering and PACM

Princeton

University Yeast Glycolytic Oscillations (Wolf and Heinrich, Biochem. J. (2000) 345, p321-334) Reaction scheme for a single cell glucose J 0 glucose ATP v 1 ADP glyceraldehyde-3-P/ dihydroxyacetone-P NAD + v 2 NADH NADH 1,3-bisphospho-glycerate v 3 ADP ATP pyruvate/acetaldehyde J pyruvate/acetaldehyde ex cytosol v 5 NAD + NADH v 7 v 6 glycerol Notation: A 2 A 3 N 1 N 2 S 1 S 2 - ADP - ATP, A - NAD + 2 +A - NADH, N - glucose 1 3 = A(const) +N 2 = N(const) - glyceraldehyde-3-P/ dihydroxyacetone-P S 3 - 1,3-bisphospho -glycerate S 4 S 4 ex - pyruvate/acetaldehyde - pyruvate/acetaldehyde ex J 0 - influx of glucose J - outflux of pyruvate/ acetaldehyde Reaction rates: NAD v 4 + ethanol external environment v v v v v v v 1 2 3 4 5 6 7 = k 1 S 1 A 3 [1+(A 3 /K I ) q ] -1 = k 2 S 2 N 1 = k 3 S 3 A 2 = k 4 S 4 N 2 = k 5 A 3 = k 6 S 2 N 2 = kS 4 ex

Department of Chemical Engineering and PACM

Princeton

University

Yeast Glycolytic Oscillations

Coupled ODEs for multicellular species concentrations

dS 1 , dt i



J 0 , i



v 1 , i



J 0 , i



k 1 S 1 , i A 3 , i

   

1



dS 2 , i dt



2 v 1 , i



v 2 , i



v 6 , i



2 k 1 S 1 , i A 3 , i

   

1



A 3 , i K I A 3 , i K I q

   

1 q

   

1



k 2 S 2 , i ( N



N 2 , i )



k 6 S 2 , i N 2 , i dS 3 , i



v 2 , i



v 3 , i



k 2 S 2 , i ( N



N 2 , i )



k 3 S 3 , i ( A



A 3 , i ) dt dS 4 , i



v 3 , i



v 4 , i



J i



k 3 S 3 , i ( A



A 3 , i )



k 4 S 4 , i N 2 , i



J i dt dN 2 , i



v 2 , i



v 4 , i



v 6 , i



k 2 S 2 , i ( N



N 2 , i )



k 4 S 4 , i N 2 , i



k 6 S 2 , i N 2 , i dt dA 3 , i dt dS 4 , ex dt

 

2 v 1 , i



2 v 3 , i



v 5 , i

 

2 k 1 S 1 , i A 3 , i

   

1

  

M i M

 

1 J i



v 7

 

M i M

 

1



( S 4 , i



A 3 , i K I S 4 , ex )

 

S 4 , ex q

   

1



2 k 3 S 3 , i ( A



A 3 , i )



k 5 A 3 , i

Department of Chemical Engineering and PACM

Princeton

University

Yeast Glycolytic Oscillations

Heterogeneity of the coupled model:

J 0



J 0

 

J

 Polynomial Chaos expansion of the solution:

x

(



, t )



( S 1 (



, t ), S 2 (



, t ), S 1 (



, t ), S 2 (



, t ), N 2 (



, t ), A 3 (



, t )) T

x



t



j

3   0

α

j

( ) 

j

Coarse variables:

α

j and S 4 ex (25 variables totally) Lifting:

x



i t



j

3   0

α

j

 Fine variables:

S

1 ( 

i

),

j



i S

2 ( 

i

),

i

 1, 2,...,

M S

3 ( 

i

),

S

4 ( 

i

),

N

2 ( 

i

),

A

3 ( 

i

),

S

4

ex

M – number of cells; M = 1000 (6M+1 variables) Restriction: Minimizing

|| x

(  ,

t

) 

j

3   0

α

j

(

t

) 

j

(  ) ||

L

2 to obtain

α

j

Department of Chemical Engineering and PACM

Princeton

University

Yeast Glycolytic Oscillations

Full ensemble simulation

J 0



2 .

3 ,



J



0 .

001

Department of Chemical Engineering and PACM

Princeton

University

Yeast Glycolytic Oscillations

Projective integration of PC coefficients S 1 S 2 S 3 S 4 N 2 A 3 t t Time histories of zeroth-order PC coef’s restricted from the full-ensemble simulation A 3 N 2 A phase map of zeroth-order PC coef’s through projective integration

Department of Chemical Engineering and PACM

Princeton

University

Yeast Glycolytic Oscillations

Limit-cycle computation

Poincaré section

Poincaré section: PC coef. of N 2 cell is constant

limit cycle

In the space of coarse variables, zeroth-order is constant

fixed point

A 3 In the space of fine variables, N 2 of a single ___ limit cycle in the space of PC coefficients xxx restricted PC coefficients of a limit cycle of the full-ensemble simulation N 2 Phase maps of zeroth-order PC coef’s through limit-cycle computation

Department of Chemical Engineering and PACM

Princeton

University

Yeast Glycolytic Oscillations

Flow map 

T

(

x

0 ) 

x

(

T

;

x

0 ), Stability of limit cycles

T

– period of the limit cycle x - PC coefficients o - full-ensemble simulation real Eigenvalues of Jacobians of the flow maps in the coarse and fine variable spaces

Department of Chemical Engineering and PACM

Princeton

University

Yeast Glycolytic Oscillations

Free oscillator

J

0  2.1, 

J

 0.08

Zeroth-order PC coefficient of N 2 (CPI)

Department of Chemical Engineering and PACM

Princeton

University Remarks • EF UQ is applied to the biological oscillations.

• The case of only one random parameter is studied. The work can be possibly extended to situations with multiple random parameters or random processes. More advantageous sampling techniques, such as cubature rules and Quasi Monte Carlo, may be used. Reference Bold, K.A., Zou, Y., Kevrekidis, I.G., and Henson, M.A., Efficient simulation of coupled biological oscillators through Equation-Free Uncertainty Quantification, in preparation, available at http://arnold.princeton.edu/~yzou

Department of Chemical Engineering and PACM