Center for Uncertainty Quantification Thuwal, KSA, Jan. 6, 2015 1 - - PowerPoint PPT Presentation

Hybrid Multilevel Monte Carlo Simulation of Stochastic Reaction Networks

Alvaro Moraes joint work with Ra´ ul Tempone and Pedro Vilanova http://sri-uq.kaust.edu.sa

Center for Uncertainty Quantification

Thuwal, KSA, Jan. 6, 2015

1 / 39

Our Contribution

This presentation is based on: Ref 1: A. Moraes, R. Tempone and P. Vilanova,“Hybrid Chernoff Tau-leap”, SIAM Multiscale Modeling and Simulation, 12(2):581-615, 2014. [Moraes et al., 2014a] Ref 2: A. Moraes, R. Tempone and P. Vilanova,“Multilevel Hybrid Chernoff Tau-leap”, preprint arXiv: 1403.2943, 2014 (submitted). [Moraes et al., 2014c] Ref 3: A. Moraes, R. Tempone and P. Vilanova,“Multilevel adaptive reaction-splitting simulation method for stochastic reaction networks”, preprint arXiv:1406.1989, 2014 (submitted). [Moraes et al., 2014b]

2 / 39

Part I

1 State the problem 2 Exact and approximate simulation of Stochastic Reaction

Networks

3 Chernoff Tau-leap 4 Hybrid paths 5 Global error control 6 Work optimization strategies 7 Some numerical results

3 / 39

Statement of the problem

Let X be a Stochastic Reaction Network X=(X1, . . . , Xd) : [0, T] → Zd

+

described by J reaction channels, Rj := (νj, aj), where

• νj ∈ Zd is called stoichiometric vector
• aj : Zd

+ → R+ is called propensity function

such that P

• X(t + dt) = x+νj
• X(t) = x
• = aj(x)dt + o (dt) .

Goal: Given an observable g : Zd

+ → R, a tolerance TOL>0, and a

small number α > 0, find an estimator M of E [g(X(T))] such that P (|E [g(X(T))] − M| > TOL) < α with near-optimal computational work.

4 / 39

Total error |E [g(X(T))] − M| vs TOL

Figure : P (|E [g(X(T))] − M| > TOL) < α = 5%

5 / 39

Example: Gene transcription and translation

• ∅ 25

− → M, a gene is being transcribed into mRNA.

• M

1000

− − − → M+P, mRNA is then being translated into proteins.

• P+P

0.001

− − − → D, finally the proteins produce stable Dimers.

• M

0.1

− − → ∅, P

1

− → ∅ degradation of mRNA and proteins, resp. Estimate the expected number of Dimers at time T=1 up to certain tolerance TOL, with high probability. ν=(νj)J

j=1:=

      1 1 −2 1 −1 −1      

T

and a(X):=       25 103M 10−3P(P−1) 0.1 M P       , (1) where X(t) = (M(t), P(t), D(t)).

6 / 39

Example: Gene transcription and translation

0.2 0.4 0.6 0.8 1 1000 2000 3000 4000 5000 6000 Time Species mRNA Proteins Dimers Mean

0.2 0.4 0.6 0.8 1 10 10

1

10

2

10

3

10

4

Time Species M P D

Figure : Exact paths of the time evolution of the species in the gene transcription and translation (GTT) example.

7 / 39

The Stochastic Simulation Algorithm (SSA)[Gillespie, 1976]

It produces statistically correct path-simulations of X. Remember that X is a continuous-time Markov chain living in Zd

+. 1 In state x at time t, compute (aj(x))J j=1 and their sum

a0(x) =

1≤j≤J aj(x). 2 Simulate the time to the next reaction, τ, as an exponential

random variable with rate a0(x). Observe: E [τ] = (a0(x))−1 could be very small!.

3 Simulate independently the next reaction, νj, according to the

probability mass function values (aj(x)/a0(x))J

j=1. 4 Update: t ← t + τ and x ← x+νj. 5 Record (t, x). Return to step one if t < T, otherwise end the

simulation.

8 / 39

The Tau-leap method(TL) [Aparicio and Solari, 2001, Gillespie, 2001]

From Kurtz’s random time change representation [Kurtz, 1978], X(t) = X(0) +

J

• j=1

Yj t aj(X(s))ds

• νj,

where Yj are independent unit-rate Poisson processes, we obtain the Tau-leap method (kind of forward Euler approximation): Given ¯ X(t)=x ∈ Zd

+,

¯ X(t+τ) = x +

J

• j=1

Pj   aj(x)τ

λj

   νj, where Pj(λj) are independent Poisson random variables with rate λj. Danger: ¯ X(t+τ) may have some negative components!

9 / 39

A Chernoff bound for the Tau-leap method I

Problem: Given δ>0, find the largest τ = τ(x, δ) s.t. P ¯ X(t+τ) / ∈ Zd

+

• ¯

X(t) = x

• ≤ ChBnd(x, τ) < δ.

Idea: Use the moment generating function of a linear combination

• f independent Poisson random variables to produce a Chernoff

bound ChBnd(x, τ). For i = 1, . . . , d, solve numerically sup

τi>0

inf

s>0 exp

• −sxi + τi

J

• j=1

aj(x)(e−sνji−1)

• = δ/d.

We define τCh := min{τ1, . . . , τd}.

10 / 39

What is a Chernoff bound?

Single-reaction case Let Q ∼ Poisson(λ), λ > 0. Given n ≥ 0 integer, the Chernoff bound is given by P (Q ≥ n) ≤ exp

• n(1 − log(n/λ) − λ)
• ,

valid for λ < n (otherwise the bound is trivial). Proof: First note that for s > 0 the Markov inequality gives P (Q ≥ n) = P

• esQ ≥ esn

≤ E

• esQ

esn , P (Q ≥ n) ≤ exp

• inf

s>0{−sn + λ(es − 1)}

• .

When λ ∈ (0, n), infs>0{−sn + λ(es − 1)} is achieved at ˜ s = log(n/λ), and its value is n(1 − log(n/λ) − λ).

11 / 39

Exact pre-leaping: The Chernoff bound

4 6 8 10 Λ 106 105 104 0.001 0.01 0.1 1

Klar 1D Chernoff this work Poisson exact Gaussian approximation

Figure : Let n = 10 and λ ∈ (2, 10). Semi-logarithmic plot of P (Q(λ) ≥ n) ≤ ChBnd(n, λ) = exp

• n(1 − log(n/λ) − λ)
• . See Klar’s

bound in [Klar, 2000]. See also [Cao et al., 2005a]

12 / 39

A Hybrid Chernoff Tau-Leap Algorithm

Main issues:

• Exact algorithms like SSA and MNRM may be not

computationally feasible, since the expected inter-arrival time between transitions, τSSA(x) = (a0(x))−1, where x = X(t) could be very small.

• Approximate algorithms that evolve with fixed time steps,

like the Tau-leap, may be faster [Aparicio and Solari, 2001, Gillespie, 2001], but introduces discretization error and can jump outside Zd

+ (non physical

results). Pre-leap: adjust the time step to control the

• ne-step exit probability, possibly enforcing too small steps.

Main idea:

• We propose a hybrid algorithm that, at each time step,

adaptively switches between SSA and Tau-leap.

13 / 39

One-step switching rule

Algorithm 1 From ¯ X(t)=x take one step. T0 is the next grid point.

1: Compute τSSA. 2: if K1τSSA > T0−t then 3:

Use SSA

4: else 5:

Compute τCh

6:

if τCh ≥ K2τSSA then

7:

Use Tau-Leap

8:

else

9:

Use SSA

10:

end if

11: end if

Note: K1 is the cost of computing τCh(x) divided by the cost of taking an SSA step. K2 is the cost of taking a Chernoff Tau-leap step divided by the cost of taking an SSA step plus the cost of computing τCh(x).

14 / 39

One-step switching rule in the GTT model

10 20 30 40 5 10 15 20 25 30 35 40

mRNA Proteins Chernoff tau−leap and SSA regions

Tau−leap SSA

10 20 30 40 5 10 15 20 25 30 35 40

mRNA Proteins Chernoff tau−leap and SSA regions

Tau−leap SSA

10 20 30 40 5 10 15 20 25 30 35 40

mRNA Proteins Chernoff tau−leap and SSA regions

Tau−leap SSA

Figure : Regions of the one-step switching rule in the GTT model. The blue and red dots show the Chernoff tau-leap and the SSA regions,

• respectively. From left to right, δ = 10−2, 10−4, 10−6, respectively.

15 / 39

Hybrid Tau-leap algorithm

Algorithm 2 Given ¯ X(t0)=x0, simulates a hybrid path.

1: while t < T do 2:

method ← One-step rule

3:

if method = SSA then

4:

Apply SSA and advance one step

5:

else

6:

¯ X ′

t ← ¯

Xt + P(a( ¯ Xt)τCh)ν

7:

if ¯ X ′

t ∈ Zd + then

8:

return

9:

end if

10:

¯ Xt ← ¯ X ′

t

11:

t ← t + τCh

12:

NTL ← NTL + 1

13:

end if

14: end while

16 / 39

An exiting path is a rare event: the role of δ

Let A be the event that a hybrid trajectory arrived to final time T without exiting Zd

+. We show in Ref 1 that

P (Ac) ≤ δE [NTL] − δ2 2 (E

• N2

TL

• − E [NTL]) + o(δ2).

In practice, we use δE [NTL] as an upper bound of P (Ac). We also prove that E [NTL] is bounded for polynomial propensity functions and tends to zero when δ → 0. Remark: The role of δ is to turn Ac into a rare event. Direct sampling of hybrid paths to estimate P(Ac) is non feasible, while the estimate of E [NTL] is straightforward.

17 / 39

Global Error Decomposition

Define the Monte Carlo estimator M as M := 1 M

M

• m=1
• g( ¯

X(T))1A

• (¯

ωm). Define the computational global error, E, as E := E [g(X(T))] − M. Global Error decomposition (notation: ¯ g:=g( ¯ X(T))) E = E [g(X(T))(1A + 1Ac)] ± E [¯ g1A] − M = E [g(X(T))1Ac]

• =:EE (exit)

+ E [(g(X(T))−¯ g) 1A]

• =:EI (discretization)

+ E [¯ g1A] −M

• =:ES (statistical)

. Problem: Given TOL> 0 and α> 0, find the parameters for computing M such that |E|<TOL with confidence level 1−α, and with nearly optimal computational work.

18 / 39

Estimation procedure

We have a 3-phase estimation procedure: Phase 1 (off-line) Calibration of virtual machine-dependent quantities, for examples Ci, i = 1, 2, 3. Phase 2 Given a prescribed TOL and a confidence level α, find:

• δ (upper bound for the one-step exit probability

for controlling the global exit error).

• a time mesh of size h (for controlling the time

discretization error).

• M (the number of Monte Carlo realizations for

controlling the statistical error). such that the estimator M can be computed with near-optimal computational work. Phase 3 Estimation of E [g(X(T))].

19 / 39

About Phase 2, I

Define the expected cost of a hybrid trajectory, Ψ(h, δ) := C1E [NSSA,K1(h, δ)] + C2E [NSSA,K2(h, δ)] + C3E [NTL(h, δ)] +

J

• j=1

E

• [0,T]

CP(aj( ¯ X(s))τCh( ¯ X(s), δ))1TL( ¯ X(s))ds

• ,

we want to find an approximate solution of    minM,h,δ MΨ(h, δ) s.t. EI + EE + ES ≤ TOL . where C1, C2, C3 are machine-dependent work quantities (Alg. 1) and CP is the computational work of the Gamma simulation method for Poisson variates [Ahrens and Dieter, 1974],

20 / 39

About Phase 2, II

500 1000 1500 2 3 4 5 6 x10

−4

λ CP(λ) Poisson random variates computational work model Actual simulation runtimes Least squares fit 5 10 15 1 2 x10

−4

λ CP(λ) Poisson random variates computational work model Actual simulation runtimes Least squares fit

Figure : Computational work of generating Poisson deviates (computer and MATLAB dependent).

21 / 39

About Phase 2, III

For a given TOL > 0, we refine δ until EE ≤ TOL2,    minM,h MΨ(h, δ0) s.t. EI(h, δ0) + ES ≤ TOL−TOL2 . We now refine the mesh until EI(h0, δ0) + ES ≤ TOL−TOL2, using an intermediate value for Maux=Maux(h, δ, Ψ), to finally get M(h0, δ0) =

• CA
• S2 (g(X(T)); Ms)

TOL − TOL2 − EI(h0, δ0) 2 . Finally, if MΨ < MSSAE [NSSA∗] use Hybrid, otherwise use SSA.

22 / 39

Numerical results: global error vs predicted and actual work

10

1

10

2

10

3

10

4

10

−2

10

−1

Work (runtime, seconds) Error bound Predicted/Actual work vs. Error bound, Genes model SSA predicted Hybrid predicted SSA Hybrid

Figure : As the TOL goes to zero, the error bound goes to zero and the expected computational work of the Hybrid method tends to the expected work of the SSA. Observe the agreement between predicted and actual work. (1)

23 / 39

Part II

1 The Multilevel Monte Carlo idea 2 A result about complexity of our MLMC algorithm 3 Some numerical results

24 / 39

Introducing some notation

Consider a hierarchy of nested time meshes of the interval [0, T], indexed by ℓ = 0, 1, . . . , L, such that

1 ∆t0 is the size of the coarsest time mesh. 2 ∆tℓ = 2−ℓ∆t0, ℓ = 1, . . . , L.

Let

• ¯

Xℓ(·):= ¯ X(·; ∆tℓ, δ) be a hybrid Chernoff tau-leap path generated using a time mesh of size ∆tℓ and one-step exit probability bound δ.

• Aℓ:={¯

ω ∈ ¯ Ω : ¯ Xℓ(t) ∈ Zd

+, ∀t ∈ [0, T]}.

• gℓ := g( ¯

Xℓ(T)).

25 / 39

The MLMC estimator and Global Error Decomposition I

Consider the following telescoping decomposition (see Ref 2): E [gL1AL] = E [g01A0] +

L

• ℓ=1

E

• gℓ1Aℓ − gℓ−11Aℓ−1
• ,

which motivates our MLMC estimator of E [g(X(T))], ML := 1 M0

M0

• m0=1

g01A0(ωm0) +

L

• ℓ=1

1 Mℓ

Mℓ

• mℓ=1

[gℓ1Aℓ − gℓ−11Aℓ−1](ωmℓ). We define the computational global error, EL, as EL := E [g(X(T))] − ML.

26 / 39

The MLMC estimator and Global Error Decomposition II

Now, consider the following decomposition of EL EL = E

• g(X(T))(1AL + 1Ac

L)

• ± E [gL1AL] − ML

= E

• g(X(T))1Ac

L

• =:EE,L (exit)

+ E [(g(X(T))−gL) 1AL]

• =:EI,L (discretization)

+ E [gL1AL] −ML

• =:ES,L (statistical)

. Problem: Given TOL> 0, find the parameters for computing ML such that |EL|<TOL with high probability, and with nearly optimal computational work. Issues:

1 Simulate coupled pairs (gℓ, gℓ−1) for ℓ = 1, . . . , L. 2 Estimate accurately the global error components.

27 / 39

Computational Complexity I

Theorem: Computational complexity of the Multilevel Hybrid Chernoff Tau-Leap is wL(TOL) =

• CA

θ

L

• ℓ=0
• Vℓψℓ

2 TOL−2 = O

• TOL−2

. where V0 := Var [g01A0], Vℓ := Var

• gℓ1Aℓ − gℓ−11Aℓ−1
• , ℓ ≥ 1

and ψℓ is the expected computational work per level of simulating coupled paths.

28 / 39

Numerical results: global error vs actual work

10

1

10

2

10

3

10

4

10

−2

10

−1

Actual work (runtime, seconds) Error bound Actual work vs. Error bound, Genes model SSA Hybrid Hybrid ML slope 1/2 Figure : Actual work for each one of the one hundred calibrations.

29 / 39

Part III

1 This is a new extension of our hybrid Multilevel Monte Carlo

method.

2 The Reaction-splitting strategy is intended to deal with stiff

problems.

3 A novel Control Variate based on a deterministic-time change

approximation.

4 We preserve the computational complexity O

• TOL−2

.

5 Some numerical results.

30 / 39

Example: Virus kinetics

• E

1

− → E+G, the viral template (E) forms a viral genome (G).

• G

0.025

− − − → E, the genome generates a new template.

• E

1000

− − − → E+S, a viral structural protein (S) is generated.

• G+S

7.5×10−6

− − − − − − → V , the virus (V) is produced.

• E

0.25

− − → ∅, S

2

− → ∅ degradation reactions. ν:=         1 −1 1 1 −1 −1 1 −1 −1        

tr

and a(X):=         E 0.025 G 1000 E 7.5×10−6G S 0.25 E 2 S         , where X(t) = (G(t), S(t), E(t), V (t)).

31 / 39

Example: Virus kinetics

Estimate the expected number of produced viruses (V) at time T=20 (in days) up to a given TOL with high confidence.

5 10 15 20 1000 2000 3000 4000 5000 6000

20 exact paths Time Number of particles

G S E V

5 10 15 20 10 10

1

10

2

10

3

10

4

20 exact paths Time Number of particles (log scale)

G S E V

Figure : Exact paths of the time evolution of the species in the virus kinetics (VRK) example.

32 / 39

Numerical results: global error vs predicted work

10

4

10

6

10

8

10

−3

10

−2

10

−1

Predicted work Error bound Predicted work vs. Error bound, Virus model SSA Mixed ML slope 1/2 Asymptotic Figure : The quotient ˆ WMix/ ˆ WHyb starts at 2% and reaches its asymptotic value of 6%.

33 / 39

Final comments

• A novel mixed algorithm that allows us to approach problems

in which there is a natural partition on the reaction channels.

• An 3-phase estimation procedure that provides the elements

for the multilevel Monte Carlo setting (one-step exit probabilities, time meshes and number of mixed paths), which

• ptimizes the computational work.
• The complexity of the method is of order O
• TOL−2

, like in the SSA, but with a better constant!

• A control variate for the variance of the QoI when using

single-level paths (applies for the hybrid and mixed methods).

• We achieved speedups of 102 − 104 with the control variate

and the mixed MLMC for large systems (200 reactions & species).

34 / 39

References I

Ahrens, J. and Dieter, U. (1974). Computer methods for sampling from gamma, beta, Poisson and bionomial distributions. Computing, 12:223–246. Anderson, D. and Higham, D. (2012). Multilevel Monte Carlo for continuous Markov chains, with applications in biochemical kinetics. SIAM Multiscal Model. Simul., 10(1). Anderson, D. F. (2007). A modified next reaction method for simulating chemical systems with time dependent propensities and delays. The Journal of Chemical Physics, 127(21):214107.

35 / 39

References II

Aparicio, J. P. and Solari, H. (2001). Population dynamics: Poisson approximation and its relation to the langevin processs. Physical Reveiw Letters, 86(18):4183–4186. Cao, Y., Gillespie, D., and Petzold, L. (2005a). Avoiding negative populations in explicit Poisson tau-leaping. The Journal of Chemical Physics, 123:054104. Cao, Y., Gillespie, D. T., and Petzold, L. R. (2005b). The slow-scale stochastic simulation algorithm. The Journal of Chemical Physics, 122(1):–. Giles, M. (2008). Multi-level Monte Carlo path simulation. Operations Research, 53(3):607–617.

36 / 39

References III

Gillespie, D. T. (1976). A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of Computational Physics, 22:403–434. Gillespie, D. T. (2001). Approximate accelerated stochastic simulation of chemically reacting systems. Journal of Chemical Physics, 115:1716–1733. Karlsson, J., Katsoulakis, M., Szepessy, A., and Tempone, R. (2010). Automatic Weak Global Error Control for the Tau-Leap Method. preprint arXiv:1004.2948v3, pages 1–22.

37 / 39

References IV

Klar, B. (2000). Bounds on tail probabilities of discrete distributions.

• Probab. Eng. Inf. Sci., 14:161–171.

Kurtz, T. G. (1978). Strong approximation theorems for density dependent markov chains. Stochastic Processes and their Applications, 6(3):223 – 240. Kurtz, T. G. (1982). Representation and approximation of counting processes.

• Adv. Filtering OptimalStochastic Control 42.

Moraes, A., Tempone, R., and Vilanova, P. (2014a). Hybrid chernoff tau-leap. Multiscale Modeling & Simulation, 12(2):581–615.

38 / 39

References V

Moraes, A., Tempone, R., and Vilanova, P. (2014b). A multilevel adaptive reaction-splitting simulation method for stochastic reaction networks. submitted to SIAM Journal of Scientific Computing, preprint arXiv:1406.1989. Moraes, A., Tempone, R., and Vilanova, P. (2014c). Multilevel hybrid Chernoff tau-leap. arXiv:1403.2943.

39 / 39