Stochastic projective methods for simulating stiff chemical reacting systems

Abstract

In this paper, stochastic projective methods are proposed to improve the stability and efficiency in simulating stiff chemical reacting systems. The efficiency of existing explicit tau-leaping methods can often severely be limited by the stiffness in the system, forcing the use of small time steps to maintain stability. The methods presented in this paper, namely stochastic projective (SP) and telescopic stochastic projective (TSP) method, can be considered as more general stochastic versions of the recently developed stable projective numerical integration methods for deterministic ordinary differential equations. SP and TSP method are developed by fully re-interpreting and extending the key projective integration steps in the deterministic regime under a stochastic context. These new stochastic methods not only automatically reduce to the original deterministic stable methods when applied to simulating ordinary differential equations, but also carry the enhanced stability property over to the stochastic regime. In some sense, the proposed methods are stochastic generalizations to their deterministic counterparts. As such, SP and TSP method can adopt a much larger effective time step than is allowed for explicit tau-leaping, leading to noticeable runtime speedup. The explicit nature of the proposed stochastic simulation methods relaxes the need for solving any coupled nonlinear systems of equations at each leaping step, making them more efficient than the implicit tau-leaping method with similar stability characteristics. The efficiency benefits of SP and TSP method over the implicit tau-leaping is expected to grow even more significantly for large complex stiff chemical systems involving hundreds of active species and beyond.

Introduction

In microscopic systems where several active chemical reactions exist, populations of some reactant molecules can be small enough to make their inherent random fluctuations recognizable and significant. As a result, the system is dominated by stochastic molecular behaviors, which classical deterministic approaches to simulation of chemical reacting systems are not able to capture. The stochastic simulation algorithm (SSA) [1], [2] proposed by Gillespie numerically simulates the time evolution of coupled chemical reactions in an exact manner based on the microphysical premise that underlies the chemical master equation. It requires that the numbers of molecules be updated upon every reaction. While the random characteristic of chemical reactions is traced accurately, high computational time is required especially when some of the species exist in large numbers.

In light of the high computational complexity of SSA, explicit tau-leaping method was proposed by Gillespie [3] as a tradeoff for efficiency. The general idea is to fire as many reactions at each time interval as allowed by the condition that tau, or the time step size, is small enough so that no propensity function undergoes appreciable changes in its value. The state of species is updated according to the number of firings for each reaction, which follows a Poisson distribution. While the accuracy of the tau-leaping method can be controlled with the error tolerance specified by the user, there is no guarantee for the global stability, especially for stiff chemical reacting systems.

To tackle the stiffness of stochastic simulations, several approaches have been developed. The implicit leaping method including the implicit tau-leaping method and the trapezoidal tau-leaping method have been proposed by Rathinam et al. in [4] and Cao et al. in [5], respectively. Due to the inherent nature of variance suppression, these methods allow for much larger step sizes without causing instability, speeding up the simulation dramatically. Like the implicit Euler method, the implicit tau-leaping and trapezoidal tau-leaping methods need to solve coupled nonlinear equations using an iterative method such as Newton–Ransom method at each time step. For a small chemical reactive system involving only a few species, the overhead of solving nonlinear equations is more than compensated by using larger time steps in the simulation. However, it is expected that this benefit would diminish or even become negative for systems involving many species where solving nonlinear equations is a significant cost. There also exist approaches that partition a chemical reacting system into subsets of fast and slow reactions, for example the slow-scale SSA proposed by Cao et al. [6] and the nested SSA proposed by Weinan E et al. [7]. With the stochastic partial equilibrium assumption for fast reactions, in these methods only the slow reactions are simulated either by exact SSA or explicit tau-leaping method. Their propensity functions involving fast species are estimated either analytically or numerically. Since accurate simulations for fast reactions are avoided, great improvement of efficiency can be obtained. Note that one needs an efficient and automatic procedure to detect the reactions in equilibrium. These multi-scale stochastic simulation algorithms work rather efficiently when there are distinctive gaps between the time scales involved in the chemical reacting system. In another line of work, the R-leaping method [8] has been proposed by Auger et al. as an accelerated stochastic simulation algorithm. This algorithm proceeds the simulations by a predefined number of reaction firings at each step instead of a predefined time interval, with the number of reaction firings sampled from correlated binomial distributions. Systems with disparate rates are efficiently handled by sorting the propensity functions and sampling first with the large propensity functions. As a result, the average sampling times for binomial distributions are significantly smaller than the number of existing reactions. Mjolsness et al. [9] proposed ER-leap as an accelerated exact stochastic simulation method based on R-leap method. It was demonstrated to be exact and have an asymptotically sublinear simulation runtime with respect to molecule number. Bayati et al. [10] proposed the falvorized SSA method which is defined by the composition of two SSA steps. Fast variables are turned off in the second step so as to speed up simulations.

In this paper, we address the stochastic simulation challenges imposed by stability and stiffness issues along the line of tau-leaping methods. This is in the sense that the proposed methods are not explicitly constructed as a multi-scale method, but a stochastic integration formula with enhanced stability. Our starting point is the recently developed projective integration methods targeting stiff ODEs in the deterministic regime [11], [12]. The proposed stochastic projective (SP) and telescopic stochastic projective (TSP) method are explicit in nature just as the explicit tau-leaping method, but with improved stability via more sophisticated multi-step multi-level time marching. The effective step size of SP and TSP is no longer limited by the fastest time scale in the system, allowing for much larger step sizes without jeopardizing the stability. In a different angle, SP and TSP method can be considered as more general, stochastic versions of their deterministic counterparts; they are constructed by re-interpreting and extending the key projective integration steps in the deterministic regime under a stochastic context. When being applied to simulating deterministic ODEs, these methods automatically reduce to the deterministic projective and telescopic projective integration methods and hence maintain the good stability established thereof. Intuitively, such characteristic is a very desirable property for any stable stochastic method: it shall be at least stable for deterministic ODEs, which are special cases of the chemical master equation (CME).

As will be shown theoretically and experimentally in the paper, the improved stability of the proposed methods is also carried over the stochastic regime, making them appealing improvements over the explicit tau-leaping method. For a reversible isomerization test problem, we theoretically show that if certain stability condition is satisfied, the mean of the species population given by the proposed stochastic projective method converges to the exact value, regardless of the leaping step size. Under the same condition, the histogram variance of the state distributions given by the method converges to the exact value as the leaping step size approaches zero. Furthermore, the stability region of the stochastic projective method is identical to that of its deterministic counterpart. The results suggest that the proposed methods well generalize the deterministic projective methods to the stochastic regime. On the computational side, unlike the implicit tau-leaping method or trapezoidal tau-leaping method, the presented SP and TSP method are explicit and hence no coupled nonlinear equations need to be solved during the course of simulation. This represents a significant runtime efficiency benefit, especially for large reaction networks with many species.

The rest of the paper is organized as follows. In Section 2, we first provide a background for the SSA, tau-leaping methods, and deterministic projective and telescopic integration methods. The stochastic projective method and its telescopic version are then derived in Section 3 and Section 4, respectively. We then present the theoretical analysis and the experimental results of two numerical examples to verify the new approach in Sections 5 Absolute stability and convergence for the stochastic projective method, 6 Numerical examples, followed by concluding remarks in Section 7.