Variants of Stochastic Simulation Algorithm

1. Variants of Stochastic Simulation Algorithm Henry Ato Ogoe Department of Computer Science Åbo Akademi University

2. The Stochastic Framework • Assume N molecular species {s1,...,SN} • State Vector X (t) = (X1(t),…,XN (t)) where Xi (t) is the number of molecules of species Si at time t. • M reaction channels R1,…,RN • Assume system is well-stirred and in thermal equilibrium

3. The Stochastic Framework • Dynamics of reaction channel Ri is characterized by its • propensity functionaj, and • state change vectorvj=(v1j,…,vNj), where vij gives change in population of Si induced by Rj, such that • aj(x)dt is the probability that, given X(t) = x, one reaction will occur in the next infinitesimal time interval [t ,t + dt] • R(x) is a jump Markov process

4. The Stochastic Framework • The time evolution of the probabilities of each state is defined by the Chemical Master Equation (CME); • where P(x,t|x0,t0) is the probability that X(t)=x given X(t0) = x0 • CME is impractical to solve especially for large systems • Alternative approaches???

5. Alternative Approaches to the CME • Exact Simulations • Inexact Simulations/Approximations

6. Exact Stochastic Simulation • Startingfrom the initial states, X(t0) the SSA simulated the trajectory by repeatedly updating the states after estimating • τ, the time the next reaction will fire, and • μ,the index of the firing reaction • Both τ and μ can be estimated probabilistically from the probability density function P(μ,τ) that the next reaction is μ and it occurs at τ.

7. Exact Stochastic simulation • Let • It can be shown that • Integrating P(μ,τ) over all τ from 0 to ∞ P(μ = j) = aj/a0 • Summing P(μ,τ) over all μ The two distributions above leads to Gillespie‘s SSA and other mathematically equivalent

8. variants with different computational efficiency

9. First Reaction Method (FRM)-Gillespie, 1977 • Generate a putative time τk for each reaction channel Rk according to • where k = 1,…,M; r1,…,rM are M statistically independent random samplings of U(0,1) • τ = min{τ1,…,τM} • μ = index ofmin{τ1,…,τM} • Update X X + Vμ

10. Flaws ???? • Uses M random numbers per time step • Uses O (M) to update the ak’s • Uses O (M) to identify smallest τμ

11. Direct Method (DM)-Gillespie, 1977 • Draw two independent samples r1 and r2 from U(0,1) • The index of the firing reaction is the smallest integer satisfying

12. Flaws???? • Unnecessary recalculation of all propensities • Slow, search depth (the no. of steps taken to identify ) ≈O (M)

13. Next Reaction Method (NRM) –Gibson & Bruck (2000)