Multilevel Hybrid Split Step Implicit Tau-Leap for Stochastic - PowerPoint PPT Presentation

Multilevel Hybrid Split Step Implicit Tau-Leap for Stochastic Reaction Networks (SRNs) Chiheb Ben Hammouda, Alvaro Moraes and Ra´ ul Tempone DCSE Fall School on ROM and UQ, November 4–8, 2019 1/24

Outline 1 Introduction 2 Multilevel Hybrid Split Step Implicit Tau-Leap (SSI-TL) Estimator 3 Numerical Results 2/24

SRNs applications: epidemic processes [Anderson and Kurtz, 2015] and virus kinetics [Hensel et al., 2009],. . . Biological Models Chemical reactions In-vivo population control: the the expected number of molecules. expected number of proteins. Figure 1.1: DNA transcription and Figure 1.2: Chemical reaction mRNA translation [Briat et al., 2015] network [Briat et al., 2015] 3/24

Statement of the Forward Problem A stochastic reaction network (SRN) is a continuous-time Markov Chain X defined on a probability space (Ω , F , P ) X = ( X 1 , . . . , X d ) : [0 , T ] × Ω → Z d + described by J reactions channels, R j := ( ν j , a j ), where ν j ∈ Z d + : stoichiometric vector (state change vector) a j : R d + → R + : Propensity (jump intensity) functions such that � � � � X ( t ) = x P X ( t + ∆ t ) = x + ν j = a j ( x )∆ t + o (∆ t ) , j = 1 , . . . , J. (1) Goal: Given i) an initial state X 0 = x 0 , ii) a smooth scalar observable of X , g : R d → R , iii) a user-selected tolerance, TOL , and iv) a confidence level 1 − α close to 1, provide accurate estimator ˆ Q of E [ g ( X ( T ))] such that �� � � � E [ g ( X ( T ))] − ˆ � � P Q � < TOL > 1 − α, (2) with near-optimal expected computational work and for a class of systems characterized by having simultaneously fast and slow time scales (Stiff systems). 4/24

Statement of the Forward Problem A stochastic reaction network (SRN) is a continuous-time Markov Chain X defined on a probability space (Ω , F , P ) X = ( X 1 , . . . , X d ) : [0 , T ] × Ω → Z d + described by J reactions channels, R j := ( ν j , a j ), where ν j ∈ Z d + : stoichiometric vector (state change vector) a j : R d + → R + : Propensity (jump intensity) functions such that � � � � X ( t ) = x P X ( t + ∆ t ) = x + ν j = a j ( x )∆ t + o (∆ t ) , j = 1 , . . . , J. (1) Goal: Given i) an initial state X 0 = x 0 , ii) a smooth scalar observable of X , g : R d → R , iii) a user-selected tolerance, TOL , and iv) a confidence level 1 − α close to 1, provide accurate estimator ˆ Q of E [ g ( X ( T ))] such that �� � � � E [ g ( X ( T ))] − ˆ � � P Q � < TOL > 1 − α, (2) with near-optimal expected computational work and for a class of systems characterized by having simultaneously fast and slow time scales (Stiff systems). 4/24

Simulation of SRNs Pathwise-Exact methods: ◮ Stochastic simulation algorithm (SSA) [Gillespie, 1976]. ◮ Modified next reaction algorithm (MNRA) [Anderson, 2007]. Pathwise-approximate methods: ◮ Explicit tau-leap [Gillespie, 2001, J. Aparicio, 2001]. ◮ Drift-implicit tau-leap (explicit-TL) [Rathinam et al., 2003b]. ◮ Split step implicit tau-leap (SSI-TL)[Ben Hammouda et al., 2016]. 5/24

Motivation and Contribution MLMC estimator using explicit TL [Anderson and Higham, 2012]. In systems characterized by having simultaneously fast and slow time scales, exact methods and the explicit-TL method can be very slow and numerically unstable. We propose in [Ben Hammouda et al., 2016] an efficient MLMC method that uses SSI-TL approximation at levels where the explicit-TL method is not applicable due to numerical stability issues. 6/24

The Explicit-TL Method [Gillespie, 2001, J. Aparicio, 2001] From Kurtz’s random time-change representation �� t J � � X ( t ) = x 0 + Y j λ j ( X ( s )) ds ν j , (3) t 0 j =1 where Y j are independent unit-rate Poisson processes, we obtain the explicit-TL method (kind of forward Euler approximation): Given Z exp ( t ) = z ∈ Z d + ,   J � Z exp ( t + τ ) = z +   P j  a j ( z ) τ  ν j , � �� � j =1 λ j P j ( λ j ) are independent Poisson random variables with rate λ j . Caveat: The explicit-TL is not adequate when dealing with stiff problems ( Numerical stability ⇒ τ exp threshold ≪ 1 . ) 7/24

Split Step Implicit Tau-Leap (SSI-TL) Method I The explicit-TL scheme, where z = Z exp ( t ), can be rewritten as follows: J � Z exp ( t + τ ) = z + P j ( a j ( z ) τ ) ν j j =1 J � = z + ( P j ( a j ( z ) τ ) − a j ( z ) τ + a j ( z ) τ ) ν j j =1 J J � � = z + a j ( z ) τ ν j + ( P j ( a j ( z ) τ ) − a j ( z ) τ ) ν j . (4) j =1 j =1 � �� � � �� � zero-mean noise drift 8/24

Split Step Implicit Tau-Leap (SSI-TL) Method II The idea of SSI-TL method is to take only the drift part as implicit while the noise part is left explicit. Let us define z = Z imp ( t ) and define Z imp ( t + τ ) through the following two steps: J � y = z + a j ( y ) τ ν j (Drift-Implicit step) (5) j =1 J � Z imp ( t + τ ) = y + ( P j ( a j ( y ) τ ) − a j ( y ) τ ) ν j j =1 J � = z + P j ( a j ( y ) τ ) ν j (Tau-leap step) j =1 9/24

Multilevel Monte Carlo [Giles, 2008] Goal : Reduce the variance of the standard Monte Carlo estimator. A hierarchy of nested meshes of the time interval [0 , T ], indexed by ℓ = 0 , 1 , . . . , L . h ℓ = M − ℓ h 0 : The size of the subsequent time steps for levels ℓ ≥ 1 , where M> 1 is a given integer constant and h 0 the step size used at level ℓ = 0. Z ℓ : The approximate process generated using a step size of h ℓ . Consider now the following telescoping decomposition of E [ g ( Z L ( T ))]: L � E [ g ( Z L ( T ))] = E [ g ( Z 0 ( T ))] + E [ g ( Z ℓ ( T )) − g ( Z ℓ − 1 ( T ))] . (6) ℓ =1 By defining  N 0 �  ˆ 1 Q 0 := g ( Z 0 , [ n 0 ] ( T ))   N 0  n 0 =1 (7) N ℓ � � �  ˆ 1  Q ℓ := g ( Z ℓ, [ n ℓ ] ( T )) − g ( Z ℓ − 1 , [ n ℓ ] ( T )) ,   N ℓ n ℓ =1 we arrive at the unbiased MLMC estimator, ˆ Q , of E [ g ( Z L ( T ))]: L � ˆ ˆ Q := Q ℓ . (8) ℓ =0 10/24

Multilevel Hybrid SSI-TL I Our multilevel hybrid SSI-TL estimator is defined as: L int − 1 L � � Q := ˆ ˆ Q ℓ + ˆ ˆ ˆ Q L imp + Q L int + Q ℓ , (9) c ℓ = L imp ℓ = L int +1 +1 c where  N i,L imp  c  � ˆ 1 g ( Z imp  Q L imp := , [ n ] ( T ))   N i,L imp L imp  c c  n =1 c   N ii,ℓ  � �  � + 1 ≤ ℓ ≤ L int − 1 ˆ g ( Z imp ℓ, [ n ℓ ] ( T )) − g ( Z imp , L imp 1  Q ℓ := ℓ − 1 , [ n ℓ ] ( T ))  c  N ii,ℓ n ℓ =1 N ie,L int � �  �  ˆ g ( Z exp L int , [ n ] ( T )) − g ( Z imp 1  Q L int := L int − 1 , [ n ] ( T ))   N ie,L int  n =1    N ee,ℓ � �  , L int + 1 ≤ ℓ ≤ L.  � ˆ 1 g ( Z exp ℓ, [ n ℓ ] ( T )) − g ( Z exp  Q ℓ := ℓ − 1 , [ n ℓ ] ( T ))   N ee,ℓ n ℓ =1 (10) 11/24

Coupling (Idea) [Kurtz, 1982, Anderson and Higham, 2012] To couple two Poisson random variables, P 1 ( λ 1 ), P 2 ( λ 2 ), with rates λ 1 and λ 2 , respectively, we define λ ⋆ := min { λ 1 , λ 2 } and we consider the decomposition � P 1 ( λ 1 ) := Q ( λ ⋆ ) + Q 1 ( λ 1 − λ ⋆ ) P 2 ( λ 2 ) := Q ( λ ⋆ ) + Q 2 ( λ 2 − λ ⋆ ) where Q ( λ ⋆ ), Q 1 ( λ 1 − λ ⋆ ) and Q 2 ( λ 2 − λ ⋆ ) are three independent Poisson random variables. We have small variance between the coupled processes �� �� Var [ P 1 ( λ 1 ) − P 2 ( λ 2 )] = Var Q 1 ( λ 1 − λ ⋆ ) − Q 2 ( λ 2 − λ ⋆ ) = | λ 1 − λ 2 | . Observe: If P 1 ( λ 1 ) and P 2 ( λ 2 ) are independent, then, we have a larger variance Var [ P 1 ( λ 1 ) − P 2 ( λ 2 )] = λ 1 + λ 2 . 12/24

Multilevel Hybrid SSI-TL II L int − 1 L � � Q := ˆ ˆ Q ℓ + ˆ ˆ ˆ Q L imp + Q L int + Q ℓ , (11) c ℓ = L imp ℓ = L int +1 +1 c 1 L imp , the coarsest discretization level. c 2 L int , the interface level. 3 L , the finest discretization level. , { N ii,ℓ } L int − 1 4 N := { N i,L imp +1 , N ie,L int , { N ee,ℓ } L ℓ = L int +1 } , the number of ℓ = L imp c c samples per level. 13/24

Estimation procedure I Coarsest discretization level , L imp , is determined by the numerical c stability constraint of our MLMC estimator, two conditions must be satisfied: The stability of a single path: determined by a linearized stability analysis of the backward Euler method applied to the deterministic ODE model corresponding to our system. The stability of the variance of the coupled paths of our MLMC estimator: expressed by � � � � Var g ( Z L imp +1 ) − g ( Z L imp ) ≪ Var g ( Z L imp ) . c c c 14/24

Estimation procedure II Total number of levels , L , and the set of the number of samples per level , N , are selected to satisfy the accuracy constraint given by (2), with near-optimal expected computational work. As a result, the MLMC algorithm should bound the bias and the statistical error as follows [Collier et al., 2014]: � �� � g ( X ( T )) − ˆ � � � E Q � ≤ (1 − θ ) TOL, (12) � θ TOL � 2 � � ˆ Var Q ≤ (13) C α for some given confidence parameter, C α , such that Φ( C α ) = 1 − α/ 2 ; here, Φ is the cumulative distribution function of a standard normal random variable. 15/24