Fully implicit tau-leaping methods for the stochastic simulation of chemical kinetics

The stochastic simulation algorithm (SSA), proposed by Gillespie, is a cardinal simulation method for the chemical kinetics. Because the SSA simulates every reaction event, the amount of the computational time is huge when models have many reaction channels and species. This drawback motivates many attempts to improve the efficiency with the accuracy. The existing "implicit tau-leaping" procedure attempts to accelerate the exact SSA especially for stiff systems. The implicit tau-leaping method uses an implicit discretization for the mean, together with an explicit discretization of the Poisson "noise". It is therefore a partially implicit method. In this paper we propose three fully implicit tau-leaping methods that treat implicitly both the mean part and the variance of the Poisson variables. The three methods considered below are the backward Euler for the mean and backward Euler for the variance of the Poisson variables, trapezoidal for both the mean and the variance of the Poisson variables, and backward Euler for the mean and trapezoidal for the variance of the Poisson variables. These methods are motivated by the theory of weakly convergent discretizations of stochastic differential equations. Numerical results demonstrate the performance of the new fully implicit methods.