This paper builds a convergence analysis of explicit tau-leaping schemes for simulat-ing chemical reactions from the viewpoint of stochastic differential equations. Mathematically, thechemical reaction process is a pure jump process on a lattice with state-dependent intensity. Thestochastic differential equation form of the chemical master equation can be given via Poisson ran-dom measures. Based on this form, different types of tau-leaping schemes can be proposed. In orderto make the problem well-posed, a modified explicit tau-leaping scheme is considered. It is shownthat the mean square strong convergence is of order 1/2 and the weak convergence is of order 1 forthis modified scheme. The novelty of the analysis is to handle the non-Lipschitz property of thecoefficients and jumps on the integer lattice.

In this paper, we present a class of explicit numerical methods for stiff Itˆo stochasticdifferential equations (SDEs). These methods are as simple to program and to use as the well-knownEuler-Maruyama method, but much more efficient for stiff SDEs.For such problems, it is wellknown that standard explicit methods face step-size reduction. While semi-implicit methods canavoid these problems at the cost of solving (possibly large)nonlinear systems, we show that the step-size reduction phenomena can be reduced significantly for explicit methods by using stabilizationtechniques. Stabilized explicit numerical methods calledS-ROCK (for stochastic orthogonal Runge-Kutta Chebyshev) have been introduced in [C. R. Acad. Sci. Paris, vol. 345, no. 10, 2007] asan alternative to (semi-) implicit methods for the solutionof stiff stochastic systems. In this paperwe discuss a genuine Itˆo version of the S-ROCK methods whichavoid the use of transformationformulas from Stratonovich to Itˆo calculus. This is important for many applications. We present twofamilies of methods for one-dimensional and multi-dimensional Wiener processes. We show that forstiff problems, significant improvement over classical explicit methods can be obtained. Convergenceand stability properties of the methods are discussed and numerical examples as well as applicationsto the simulation of stiff chemical Langevin equations are presented.

We aim to construct higher order tau-leaping methods for numerically simulating stochastic chemical kinetic systems in this paper. By adding a random correction to the primitive tau-leaping scheme in each time step, we greatly improve the accuracy of the tau-leaping approximations. This gain in accuracy actually comes from the reduction in the local truncation error of the scheme in the order of τ, the marching time step size. While the local truncation error of the primitive tau-leaping method is

Stabilized or Chebyshev explicit methods have been widely used in the past to solve stiffordinary differential equations. Making use of special properties of Chebyshev-like poly-nomials, these methods have favorable stability properties compared to standard explicitmethods while remaining explicit. A new class of such methods, called ROCK, introducedin [Numer. Math., 90, 1-18, 2001] has recently been extended to stiff stochastic differentialequations under the name S-ROCK [C. R. Acad. Sci. Paris, 345(10), 2007 and Commun.Math. Sci, 6(4), 2008]. In this paper we discuss the extension of the S-ROCK methodsto systems with discrete noise and propose a new class of methods for such problems, theτ-ROCK methods. One motivation for such methods is the simulation of multi-scale orstiff chemical kinetic systems and such systems are the focus of this paper, but our newmethods could potentially be interesting for other stiff systems with discrete noise. Twoversions of theτ-ROCK methods are discussed and their stability behavior is analyzed ona test problem. Compared to theτ-leaping method, a significant speed-up can be achievedfor some stiff kinetic systems. The behavior of the proposed methods are tested on severalnumerical experiments.

There are two scalings for the convergence analysis of tau-leaping methods in the literature. This paper attempts to resolve this debate in the paper. We point out the shortcomings of both scalings. We systematically develop the weak Ito-Taylor expansion based on the infinitesimal generator of the chemical kinetic system and generalize the rooted tree theory for ODEs and SDEs driven by Brownian motion to rooted directed graph theory for the jump processes. We formulate the local truncation error analysis based on the large volume scaling. We find that even in this framework the midpoint tau-leaping does not improve the weak local order for the covariance compared with the explicit tau-leaping. We propose a procedure to explain the numerical order behavior by abandoning the dependence on the volume constant V from the leading error term. The numerical examples validate our arguments. We also give a general global weak convergence analysis for the explicit tau-leaping type methods in the large volume scaling.

Recently Anderson and Mattingly [Comm. Math. Sci. 9, 301 (2011)] proposed a method which can solve chemical Langevin equations with weak second order accuracy. We extend their work to the discrete chemical jump processes. With slight modification, the method can also solve discrete chemical kinetic systems with weak second order accuracy in the large volume scaling. Especially, this method achieves higher order accuracy than both the Euler τ-leaping and mid-point τ-leaping methods in the sense that the local truncation error for the covariance is of order τ3V−1 when τ = V−β (0 < β < 1) and the system size V → ∞. We present the convergence analysis, numerical stability analysis, and numerical examples. Overall, in the authors’ opinion, the new method is easy to be implemented and good in performance, which is a good candidate among the highly accurate τ-leaping type schemes for discrete chemical reaction systems.

Motivation: Cancer is well known to be the end result of somatic mutations that disrupt normal cell division. The number of such mutations that have to be accumulated in a cell before cancer develops depends on the type of cancer. The waiting time Tm until the appearance of m mutations in a cell is thus an important quantity in population genetics models of carcinogenesis. Such models are often difficult to analyze theoretically because of the complex interactions of mutation, drift and selection. They are also computationally expensive to simulate because of the large number of cells and the low mutation rate. Results: We develop an efficient algorithm for simulating the waiting time Tm until m mutations under a population genetics model of cancer development. We use an exact algorithm to simulate evolution of small cell populations and coarse-grained τ-leaping approximation to handle large populations. We compared our hybrid simulation algorithm with the exact algorithm in small populations and with available asymptotic results for large populations. The comparison suggested that our algorithm is accurate and computationally efficient. We used the algorithm to study the waiting time for up to 20 mutations under a Moran model with variable population sizes. Our new algorithm may be useful for studying realistic models of carcinogenesis, which incorporates variable mutation rates and fitness effects.

We consider the nucleation of one-dimensional stochastic Cahn-Hilliard dynamics with the standard double well potential. We design the string method for computing the most probable transition path in the zero temperature limit based on large deviation theory. We derive the nucleation rate formula for the stochastic Cahn-Hilliard dynamics through finite dimensional discretization. We also discuss the algorithmic issues for calculating the nucleation rate, especially the high dimensional sampling for computing the determinant ratios.

We focus on the nucleation rate calculation for diblock copolymers by studying the two-dimensional stochastic Cahn--Hilliard dynamics with a Landau--Brazovskii energy functional. To do this, we devise the string method to compute the minimal energy path of nucleation events and the gentlest ascent dynamics to locate the saddle point on the path in Fourier space. Both methods are combined with the semi-implicit spectral method and hence are very effective. We derive the nucleation rate formula in the infinite-dimensional case and prove the convergence under numerical discretizations. The computation of the determinant ratio is also discussed for obtaining the rate. The algorithm is successfully applied to investigate the nucleation from the lamellar phase to the cylinder phase in the mean field theory for diblock copolymer melts. The comparison with projected stochastic Allen--Cahn dynamics is also discussed.

Genetic switching driven by noise is a fundamental cellular process in genetic regulatory networks. Quantitatively characterizing this switching and its fluctuation properties is a key problem in computational biology. With an autoregulatory dimer model as a specific example, we design a general methodology to quantitatively understand the metastability of gene regulatory system perturbed by intrinsic noise. Based on the large deviation theory, we develop new analytical techniques to describe and calculate the optimal transition paths between the on and off states. We also construct the global quasi-potential energy landscape for the dimer model. From the obtained quasi-potential, we can extract quantitative results such as the stationary distributions of mRNA, protein and dimer, the noise strength of the expression state, and the mean switching time starting from either stable state. In the final stage, we apply this procedure to a transcriptional cascades model. Our results suggest that the quasi-potential energy landscape and the proposed methodology are general to understand the metastability in other biological systems with intrinsic noise.