Multilevel Monte Carlo for the continuous time Markov chain models arising in biology David F. Anderson ∗ ∗ anderson@math.wisc.edu Department of Mathematics University of Wisconsin - Madison Work is joint with Des Higham MCQMC 2012 February 17th, 2012
Overview ◮ Most common stochastic models of biochemical reaction systems are continuous time Markov chains. ◮ Common examples include: 1. Gene regulatory networks. 2. Models of viral infection. 3. General population models (epidemic, predator-prey, etc.) ◮ These models are becoming dominant in molecular biology. Problem: extend multi-level Monte Carlo to this setting.
Why? 1. Why would the be useful? ◮ Short answer: the common simulation methods can be very computationally expensive. ◮ Many reactions happening over time-frame of interest. 2. Why is this non-trivial? ◮ In SDE case, we have a useful representation � t � t X ( t ) = X ( 0 ) + b ( X ( s )) ds + σ ( X ( s )) dW ( s ) 0 0 � t � t Z ℓ ( t ) = Z ℓ + b ( Z ℓ ◦ η ( s )) ds + σ ( Z ℓ ◦ η ( s )) dW ( s ) . 0 0 ◮ Two processes are coupled via Brownian motion W . ◮ Clear how to simulate two processes with different time steps. ◮ We need a pathwise representation to get a good coupling that: 2.1 Couples closely. 2.2 Easy to analyze.
The Poisson process A Poisson process, Y ( · ) : (a) Let { ξ i } be i.i.d. exponential random variables with parameter one. (b) Now, put points down on a line with spacing equal to the ξ i : x x x x x x x x ↔ ↔ ← → · · · ξ 1 ξ 2 ξ 3 t ◮ Let Y ( t ) denote the number of points hit by time t . ◮ In the figure above, Y ( t ) = 6.
The Poisson process Let Y be a unit rate Poisson process and λ ≥ 0. Then Y ( λ t ) is a Poisson process with parameter λ . x x x x x x x x ↔ ↔ ← → · · · λ t ξ 1 ξ 2 ξ 3 There is no reason λ needs to be constant in time, in which case �� t � Y λ ( t ) ≡ Y λ ( s ) ds 0 is an non-homogeneous Poisson process with intensity λ ( t ) .
Build up model: Random time change representation of Tom Kurtz Consider the simple system A + B → C where one molecule each of A and B is being converted to one of C . Simple book-keeping: if X ( t ) = ( X A ( t ) , X B ( t ) , X C ( t )) T gives the state at time t , − 1 , − 1 X ( t ) = X ( 0 ) + R ( t ) 1 where ◮ R ( t ) is the # of times the reaction has occurred by time t , and ◮ X ( 0 ) is the initial condition.
Build up model: Random time change representation of Tom Kurtz Assuming intensity or propensity of reaction is λ ( X ( s )) = κ X A ( s ) X B ( s ) , we can model �� t � R ( t ) = Y κ X A ( s ) X B ( s ) ds 0 where Y is a unit-rate Poisson point process. Hence �� t � X A ( t ) − 1 ≡ X ( t ) = X ( 0 ) + Y X B ( t ) − 1 κ X A ( s ) X B ( s ) ds . 0 X C ( t ) 1
Build up model: Random time change representation of Tom Kurtz • Now consider a network of reactions involving d chemical species, S 1 , . . . , S d : ν ′ k 1 S 1 + ν ′ k 2 S 2 + · · · + ν ′ ν k 1 S 1 + ν k 2 S 2 + · · · + ν kd S d − → kd S d Denote reaction vector as ζ k = ν ′ k − ν k , • The intensity (or propensity) of k th reaction is λ k : Z d ≥ 0 → R . • By analogy with before �� t � � X ( t ) = X ( 0 ) + Y k λ k ( X ( s )) ds ζ k , 0 k Y k are independent, unit-rate Poisson processes.
Example Consider a model of gene transcription and translation: 25 G → G + M , ( Transcription ) 1000 → M + P , ( Translation ) M 0 . 001 P + P → D , ( Dimerization ) 0 . 1 M → ∅ , ( Degradation of mRNA ) 1 P → ∅ ( Degradation of Protein ) . Then, if X = [ X M , X P , X D ] T , � � � t 1 0 + Y 2 X ( t ) = X ( 0 ) + Y 1 ( 25 t ) 0 1000 X M ( s ) ds 1 0 0 0 � � � t 0 + Y 3 0 . 001 X P ( s )( X P ( s ) − 1 ) ds − 2 0 1 � � � � t � � t − 1 0 + Y 5 + Y 4 0 . 1 X M ( s ) ds 0 1 X P ( s ) ds − 1 0 0 0 0
Computing Expectations Problem: Approximate E f ( X ( T )) to some desired tolerance, ǫ > 0. Easy! ◮ Simulate the CTMC exactly, ◮ generate independent paths, X [ i ] ( t ) , use the unbiased estimator n � µ n = 1 � f ( X [ i ] ( t )) . n i = 1 ◮ stop when desired confidence interval is ± ǫ . Total computational complexity = (cost per path) × (# paths) = O ( N × ǫ − 2 ) .
Tau-leaping: Euler’s method Recall: �� t � � X ( t ) = X ( 0 ) + Y k λ k ( X ( s )) ds ζ k . 0 k � t λ k ( X ( s )) ds : Tau-leaping is essentially an Euler approximation of 0 �� h � � Z ( h ) = Z ( 0 ) + Y k λ k ( Z ( s )) ds ζ k 0 k � � � ≈ Z ( 0 ) + Y k λ k ( Z ( 0 )) h ζ k k � � � d = Z ( 0 ) + Poisson λ k ( Z ( 0 )) h ζ k . k
Euler’s method Path-wise representation for Z ( t ) generated by Euler’s method is �� t � � Z ( t ) = X ( 0 ) + Y k λ k ( Z ◦ η ( s )) ds ζ k , 0 k where t n ≤ s < t n + 1 = t n + h η ( s ) = t n if is a step function giving left endpoints of time discretization.
Return to approximating E f ( X ( T )) Let n � µ n = 1 f ( Z L , [ i ] ( t )) . � n i = 1 We have µ n ) = (Bias) 2 MSE ( � + (variance) , If have an order one method Total computational complexity = (cost per path) × (# paths) = O ( ǫ − 1 × ǫ − 2 ) = O ( ǫ − 3 ) .
Benefits/drawbacks Benefits: 1. Can drastically lower the computational complexity of a problem Drawbacks: 1. Convergence results usually give order of convergence. Can’t give a precise h L . 2. Tau-leaping has problems: what happens if you go negative? 3. Gone away from an unbiased estimator.
Multi-level Monte Carlo Usual MLMC: L − 1 � E f ( X ( t )) ≈ E [ f ( Z L ( t )) − f ( Z L − 1 ( t ))]+ E [ f ( Z ℓ ( t )) − f ( Z ℓ − 1 ( t ))]+ E f ( Z ℓ 0 ( t )) . ℓ = ℓ 0 + 1 In our setting: L � E f ( X ( t )) = E [ f ( X ( t )) − f ( Z L ( t ))] + E [ f ( Z ℓ ( t )) − f ( Z ℓ − 1 ( t ))] + E f ( Z ℓ 0 ( t )) . ℓ = ℓ 0 + 1 gives an unbiased estimator.
Multi-level Monte Carlo: an unbiased estimator L � E f ( X ( t )) = E [ f ( X ( t )) − f ( Z L ( t ))] + E [ f ( Z ℓ ( t )) − f ( Z ℓ − 1 ( t ))] + E f ( Z ℓ 0 ( t )) . ℓ = ℓ 0 + 1 Estimators: n E � 1 � def = ( f ( X [ i ] ( T ) − f ( Z L , [ i ] ( T ))) , Q E n E i = 1 n ℓ � 1 � def Q ℓ = ( f ( Z ℓ, [ i ] ( T )) − f ( Z ℓ − 1 , [ i ] ( T ))) , for ℓ ∈ { ℓ 0 + 1 , . . . , L } n ℓ i = 1 n 0 � 1 � def Q 0 = f ( Z ℓ 0 , [ i ] ( T )) . n 0 i = 1 Question: what is the coupling and the computational cost?
How do we couple Poisson processes Suppose I want to generate: ◮ A Poisson process with intensity 13.1. ◮ A Poisson process with intensity 13. ◮ We could let Y 1 and Y 2 be independent, unit-rate Poisson processes, and set Z 13 . 1 ( t ) = Y 1 ( 13 . 1 t ) , Z 13 ( t ) = Y 2 ( 13 t ) , Using this representation, these processes are independent and, hence, not coupled. The variance of difference is large: Var ( Z 13 . 1 ( t ) − Z 13 ( t )) = Var ( Y 1 ( 13 . 1 t )) + Var ( Y 2 ( 13 t )) = 26 . 1 t .
How do we generate processes simultaneously Suppose I want to generate: ◮ A Poisson process with intensity 13.1. ◮ A Poisson process with intensity 13. ◮ We could let Y 1 and Y 2 be independent unit-rate Poisson processes, and set Z 13 . 1 ( t ) = Y 1 ( 13 t ) + Y 2 ( 0 . 1 t ) Z 13 ( t ) = Y 1 ( 13 t ) , The variance of difference is much smaller: Var ( Z 13 . 1 ( t ) − Z 13 ( t )) = Var ( Y 2 ( 0 . 1 t )) = 0 . 1 t .
How do we generate processes simultaneously More generally, suppose we want 1. non-homogeneous Poisson process with intensity f ( t ) and 2. non-homogeneous Poisson process with intensity g ( t ) . Let Y 1 , Y 2 , and Y 3 be independent, unit-rate Poisson processes and define �� t � �� t � Z f ( t ) = Y 1 f ( s ) ∧ g ( s ) ds + Y 2 f ( s ) − ( f ( s ) ∧ g ( s )) ds , 0 0 �� t � �� t � Z g ( t ) = Y 1 f ( s ) ∧ g ( s ) ds + Y 3 g ( s ) − ( f ( s ) ∧ g ( s )) ds , 0 0
Back to our processes �� t � � X ( t ) = X ( 0 ) + Y k λ k ( X ( s )) ds ζ k , 0 k �� t � � Z ( t ) = X ( 0 ) + λ k ( Z ◦ η ( s )) ds Y k ζ k . 0 k Now couple �� t � � λ k ( X ( s )) ∧ λ k ( Z ℓ ◦ η ℓ ( s )) ds X ( t ) = X ( 0 ) + Y k , 1 ζ k 0 k �� t � � + Y k , 2 λ k ( X ( s )) − λ k ( X ( s )) ∧ λ k ( Z ℓ ◦ η ℓ ( s )) ds ζ k 0 k �� t � � λ k ( X ( s )) ∧ λ k ( Z ℓ ◦ η ℓ ( s )) ds Z ℓ ( t ) = Z ℓ ( 0 ) + Y k , 1 ζ k 0 k �� t � � + Y k , 3 λ k ( Z ℓ ◦ η ℓ ( s )) − λ k ( X ( s )) ∧ λ k ( Z ℓ ◦ η ℓ ( s )) ds ζ k 0 k Algorithm for simulation is equivalent to simulating CTMC.
Recommend
More recommend
