Gillespie algorithm

The Gillespie algorithm samples statistically exact trajectories of a well-mixed chemical reaction network modeled by a continuous-time Markov chain. At copy-number state $x$, each reaction propensity function $a_r(x)$ is a nonnegative event rate and each stoichiometric vector ${\nu }_r$ specifies the update. Assume the rates have no explicit time dependence. Put $a_0(x)={\sum }_r{a}_r(x)$. If $a_0(x)=0$, no further chemical reaction can occur: the state is an absorbing state.

For $a_0(x)>0$, sample a waiting time $\tau$ from the exponential distribution with rate $a_0(x)$, independently choose channel $r$ with probability $a_r(x)/a_0(x)$, and update $(t,x)\leftarrow (t+\tau ,x+{\nu }_r)$. Recompute the reaction propensity functions and repeat. The no-event probability is $S(\tau )=e^{-a_0\tau }$ because $S^{\prime }=-a_{0}S$; multiplying by $a_r$ gives the joint waiting-time/channel probability density function $a_r{e}^{-a_0\tau }$. This factorization explains the two draws.^[ref]

For $A\to B$ with reaction propensity function $a(x)=k{x}_A$, the waiting time has mean $1/(k{x}_A)$, and $x_A=0$ is an absorbing state. Exactness concerns the specified chemical master equation; spatial mixing errors and explicitly time-varying rates require a different model or sampling rule.

Tau-leaping accelerates the Gillespie algorithm by approximating reaction rates as constant over an interval $\tau$, then drawing all reaction counts for that interval at once. At state $x$, with reaction propensity functions $a_r(x)$ and stoichiometric vectors ${\nu }_r$, the explicit step isThe leap condition requires the propensities to change little during the interval. Unlike event-by-event exact simulation, this introduces an approximation error controlled through the choice of $\tau$.^[ref]

For $A\to B$ with rate $k{x}_A$, the approximate number converted is Poisson with mean $k{x}_A\tau$. The exact number converted over the interval is instead binomial with parameters $x_A$ and $1-e^{-k\tau }$. The two agree to first order for small $k\tau$, but the Poisson draw can exceed $x_A$ and produce a negative count. Smaller or adaptive leaps and modified sampling rules address this limitation; freezing the rate is not automatically safe at low copy number.