D-leaping: Accelerating stochastic simulation algorithms for reactions with delays Basil Bayati, Philippe Chatelain, Petros Koumoutsakos * Chair of Computational Science, ETH Zurich, CH-8092, Switzerland article info Article history: Received 20 December 2008 Received in revised form 26 April 2009 Accepted 4 May 2009 Available online 12 May 2009 Keywords: Stochastic simulation Delayed reactions Tau-leaping R-leaping abstract We propose a novel, accelerated algorithm for the approximate stochastic simulation of biochemical systems with delays. The present work extends existing accelerated algo- rithms by distributing, in a time adaptive fashion, the delayed reactions so as to minimize the computational effort while preserving their accuracy. The accuracy of the present algo- rithm is assessed by comparing its results to those of the corresponding delay differential equations for a representative biochemical system. In addition, the fluctuations produced from the present algorithm are comparable to those from an exact stochastic simulation with delays. The algorithm is used to simulate biochemical systems that model oscillatory gene expression. The results indicate that the present algorithm is competitive with exist- ing works for several benchmark problems while it is orders of magnitude faster for certain systems of biochemical reactions. Ó 2009 Elsevier Inc. All rights reserved. 1. Introduction Several reactions in eukaryotic cells are not instantaneous, but rather their reactants are subject to transcription, process- ing, and synthesis before they can react with other chemical species [20,19]. These biochemical processes can be modeled via a system of reactions with delays. Delayed reactions can also be used as models of spatially dependent stochastic processes [9] when proteins may need to diffuse to a distant compartment in the cell in order to react with other proteins. These reactions can be formulated as a continuous-time Markov process with a discrete set of states that can be expressed by the so-called Master Equation (ME) [15,11]. Exact realizations of the ME can be obtained via the Stochastic Simulation Algo- rithm [5,10,12] (SSA). The connections between SSA and Molecular Dynamics, as well as the classical Langevin, Fokker– Planck, and reaction-rate equations, have been recently reviewed in [14]. The SSA is exact since it independently simulates all reaction events but it can be computationally expensive for large systems. In order to accelerate the SSA, several approximate algorithms have been proposed. These algorithms accelerate the SSA by either prescribing a larger time-step [13,22,7,9,8] or the number of reactions per time-step [2]. Recently, there has been interest in extending the SSA to incorporate delays. Delays in the stochastic process render it, by definition, non-Markovian, and suitable modifications to the SSA are necessary in order to produce the correct dynamics [6,3,1]. Cai [6] and Anderson [1] have proposed exact, delayed SSAs and, additionally, Leier et al. [16] have developed a delayed, accelerated, approximate, SSA (DAA-SSA). In this paper, a time-adaptive generalization of DAA-SSA is proposed (D-leaping). D-leaping is shown to converge to the Delay Differential Equation (DDE) and preserve the correct statistical fluctuations. Furthermore, the algorithm is adaptive in time, and is shown to be, for certain chemical reactions, orders of magnitude faster than the algorithm presented in [16]. The D-leaping algorithm can be combined with R-Leaping [2] and s-Leaping [13] as described in this work. 0021-9991/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2009.05.004 * Corresponding author. Tel.: +41 1 632 5258; fax: +41 1 632 1703. E-mail address: petros@ethz.ch (P. Koumoutsakos). Journal of Computational Physics 228 (2009) 5908–5916 Contents lists available at ScienceDirect Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp