An efficient method for computing steady state solutions with Gillespie's direct method

Gillespie's direct method is a stochastic simulation algorithm that may be used to calculate the steady state solution of a chemically reacting system. Recently the all possible states method was introduced as a way of accelerating the convergence of the simulations. We demonstrate that while the all possible states (APS) method does reduce the number of required trajectories, it is actually much slower than the original algorithm for most problems. We introduce the elapsed time method, which reformulates the process of recording the species populations. The resulting algorithm yields the same results as the original method, but is more efficient, particularly for large models. In implementing the elapsed time method, we present robust methods for recording statistics and empirical probability distributions. We demonstrate how to use the histogram distance to estimate the error in steady state solutions.