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 function aj, and
• state change vector vj (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,tx0,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

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

7
Exact Stochastic simulation

• Let
• It can be shown that
• Integrating P(µ,t) over all t from 0 to 8
• P(µ j) aj/a0
• Summing P(µ,t) over all µ
• The two distributions above leads to Gillespies
  SSA and other mathematically equivalent

8

• variants with different computational efficiency

9
First Reaction Method (FRM)-Gillespie, 1977

• Generate a putative time tk for each reaction
  channel Rk according to
• where k 1,,M r1,,rM are M statistically
  independent random samplings of U(0,1)
• t mint1,,tM
• µ index of mint1,,tM
• Update X X Vµ

10
Flaws ????

• Uses M random numbers per time step
• Uses O (M) to update the aks
• Uses O (M) to identify smallest tµ

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)