 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.

The Gillespie algorithm, also known as the Gillespie stochastic simulation algorithm (SSA), is a numerical method used to simulate the time evolution of systems with probabilistic events, particularly in the context of biochemical reactions. It was developed by Daniel T. Gillespie in 1976 to address the need for modeling the dynamics of chemical systems where the number of molecules is relatively small, and where stochastic effects become significant.