The tau-leaping method of simulating the stochastic time evolution of a well-stirred chemically reacting system uses a Poisson approximation to take time steps that leap over many reaction events. Theory implies that tau leaping should be accurate so long as no propensity function changes its value "significantly" during any time step tau. Presented here is an improved procedure for estimating the largest value for tau that is consistent with this condition. This new tau-selection procedure is more accurate, easier to code, and faster to execute than the currently used procedure. The speedup in execution will be especially pronounced in systems that have many reaction channels.

Reactions in real chemical systems often take place on vastly different time scales, with "fast" reaction channels firing very much more frequently than "slow" ones. These firings will be interdependent if, as is usually the case, the fast and slow reactions involve some of the same species. An exact stochastic simulation of such a system will necessarily spend most of its time simulating the more numerous fast reaction events. This is a frustratingly inefficient allocation of computational effort when dynamical stiffness is present, since in that case a fast reaction event will be of much less importance to the system's evolution than will a slow reaction event. For such situations, this paper develops a systematic approximate theory that allows one to stochastically advance the system in time by simulating the firings of only the slow reaction events. Developing an effective strategy to implement this theory poses some challenges, but as is illustrated here for two simple systems, when those challenges can be overcome, very substantial increases in simulation speed can be realized.

We show how stiffness manifests itself in the simulation of chemical reactions at both the continuous-deterministic level and the discrete-stochastic level. Existing discrete stochastic simulation methods, such as the stochastic simulation algorithm and the ͑explicit͒ tau-leaping method, are both exceedingly slow for such systems. We propose an implicit tau-leaping method that can take much larger time steps for many of these problems.

In this paper we examine the different formulations of Gillespie’s stochastic simulation algorithm (SSA) [D. Gillespie, J. Phys. Chem. 81, 2340 (1977)] with respect to computational efficiency, and propose an optimization to improve the efficiency of the direct method. Based on careful timing studies and an analysis of the time-consuming operations, we conclude that for most practical problems the optimized direct method is the most efficient formulation of SSA. This is in contrast to the widely held belief that Gibson and Bruck’s next reaction method [M. Gibson and J. Bruck, J. Phys. Chem. A 104, 1876 (2000)] is the best way to implement the SSA for large systems. Our analysis explains the source of the discrepancy.

An adjoint sensitivity method is presented for parameter-dependent differentialalgebraic equation systems (DAEs). The adjoint system is derived, along with conditions for its consistent initialization, for DAEs of index up to two (Hessenberg). For stable linear DAEs, stability of the adjoint system (for semi-explicit DAEs) or of an augmented adjoint system (for fully implicit DAEs) is shown. In addition, it is shown for these systems that numerical stability is maintained for the adjoint system or for the augmented adjoint system.