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 =( X 1 , . . . , X d ) : [0 , T ] → Z d + described by J reaction channels , R j := ( ν j , a j ), where • ν j ∈ Z d is called stoichiometric vector + → R + is called propensity function • a j : Z d such that � � � � X ( t ) = x X ( t + dt ) = x + ν j = a j ( x ) dt + o ( dt ) . P Goal: Given an observable g : Z d + → 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. 1000 • M − − − → M + P , mRNA is then being translated into proteins. 0 . 001 • P + P − − − → D , finally the proteins produce stable Dimers. 0 . 1 1 • M − − → ∅ , P → ∅ degradation of mRNA and proteins, resp. − Estimate the expected number of Dimers at time T =1 up to certain tolerance TOL , with high probability.   T   1 0 0 25 10 3 M 0 1 0         ν =( ν j ) J 10 − 3 P ( P − 1) j =1 := 0 − 2 1 and a ( X ):= ,         − 1 0 0 0 . 1 M     0 − 1 0 P (1) where X ( t ) = ( M ( t ) , P ( t ) , D ( t )). 6 / 39

Example: Gene transcription and translation 6000 4 10 mRNA 5000 Proteins 3 10 Dimers 4000 Mean M Species P Species 3000 2 D 10 2000 1 10 1000 0 0 10 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Time Time 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 Z d + . 1 In state x at time t , compute ( a j ( x )) J j =1 and their sum a 0 ( x ) = � 1 ≤ j ≤ J a j ( x ). 2 Simulate the time to the next reaction, τ , as an exponential random variable with rate a 0 ( x ). Observe: E [ τ ] = ( a 0 ( x )) − 1 could be very small!. 3 Simulate independently the next reaction, ν j , according to the probability mass function values ( a j ( x ) / a 0 ( 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], �� t J � � X ( t ) = X (0) + Y j a j ( X ( s )) ds ν j , 0 j =1 where Y j are independent unit-rate Poisson processes, we obtain the Tau-leap method (kind of forward Euler approximation): Given ¯ X ( t )= x ∈ Z d + ,   J � ¯   X ( t + τ ) = x + P j  a j ( x ) τ  ν j , � �� � j =1 λ j where P j ( λ 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. � ¯ � � � ¯ ∈ Z d X ( t + τ ) / X ( t ) = x ≤ ChBnd ( x , τ ) < δ. P + Idea: Use the moment generating function of a linear combination of independent Poisson random variables to produce a Chernoff bound ChBnd ( x , τ ). For i = 1 , . . . , d , solve numerically J � � � a j ( x )( e − s ν ji − 1) sup s > 0 exp inf − sx i + τ i = δ/ d . τ i > 0 j =1 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 � e sQ � � e sQ ≥ e sn � ≤ E P ( Q ≥ n ) = P , e sn � � s > 0 {− sn + λ ( e s − 1) } P ( Q ≥ n ) ≤ exp inf . When λ ∈ (0 , n ), inf s > 0 {− sn + λ ( e s − 1) } is achieved at ˜ s = log( n /λ ), and its value is n (1 − log( n /λ ) − λ ). 11 / 39

Exact pre-leaping: The Chernoff bound 1 0.1 0.01 Klar � 1D � 0.001 Chernoff � this work � Poisson � exact � 10 � 4 Gaussian � approximation � 10 � 5 10 � 6 10 Λ 4 6 8 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 ) = ( a 0 ( 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 Z d + (non physical results). Pre-leap: adjust the time step to control the one-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. T 0 is the next grid point. 1: Compute τ SSA . 2: if K 1 τ SSA > T 0 − t then Use SSA 3: 4: else Compute τ Ch 5: if τ Ch ≥ K 2 τ SSA then 6: Use Tau-Leap 7: else 8: Use SSA 9: end if 10: 11: end if Note: K 1 is the cost of computing τ Ch ( x ) divided by the cost of taking an SSA step. K 2 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 Chernoff tau−leap and SSA regions Chernoff tau−leap and SSA regions Chernoff tau−leap and SSA regions 40 40 40 Tau−leap Tau−leap Tau−leap SSA SSA SSA 35 35 35 30 30 30 25 25 25 Proteins Proteins Proteins 20 20 20 15 15 15 10 10 10 5 5 5 0 0 0 0 10 20 30 40 0 10 20 30 40 0 10 20 30 40 mRNA mRNA mRNA 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 ( t 0 )= x 0 , simulates a hybrid path. 1: while t < T do method ← One-step rule 2: if method = SSA then 3: Apply SSA and advance one step 4: else 5: ¯ t ← ¯ X t + P ( a ( ¯ X ′ X t ) τ Ch ) ν 6: if ¯ X ′ t �∈ Z d + then 7: return 8: end if 9: X t ← ¯ ¯ X ′ 10: t t ← t + τ Ch 11: N TL ← N TL + 1 12: end if 13: 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 Z d + . We show in Ref 1 that P ( A c ) ≤ δ E [ N TL ] − δ 2 � � N 2 − E [ N TL ]) + o ( δ 2 ) . 2 ( E TL In practice, we use δ E [ N TL ] as an upper bound of P ( A c ). We also prove that E [ N TL ] is bounded for polynomial propensity functions and tends to zero when δ → 0. Remark: The role of δ is to turn A c into a rare event. Direct sampling of hybrid paths to estimate P ( A c ) is non feasible, while the estimate of E [ N TL ] is straightforward. 17 / 39

Global Error Decomposition Define the Monte Carlo estimator M as M M := 1 � � � g ( ¯ X ( T )) 1 A (¯ ω m ) . M m =1 Define the computational global error, E , as E := E [ g ( X ( T ))] − M . g := g ( ¯ Global Error decomposition (notation: ¯ X ( T ))) E = E [ g ( X ( T ))( 1 A + 1 A c )] ± E [¯ g 1 A ] − M = E [ g ( X ( T )) 1 A c ] + E [( g ( X ( T )) − ¯ g ) 1 A ] + E [¯ g 1 A ] −M . � �� � � �� � � �� � =: E E (exit) =: E I (discretization) =: E S (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