Efficient simulation of stochastic chemical kinetics with the Stochastic Bulirsch-Stoer extrapolation method

BackgroundBiochemical systems with relatively low numbers of components must be simulated stochastically in order to capture their inherent noise. Although there has recently been considerable work on discrete stochastic solvers, there is still a need for numerical methods that are both fast and accurate. The Bulirsch-Stoer method is an established method for solving ordinary differential equations that possesses both of these qualities.ResultsIn this paper, we present the Stochastic Bulirsch-Stoer method, a new numerical method for simulating discrete chemical reaction systems, inspired by its deterministic counterpart. It is able to achieve an excellent efficiency due to the fact that it is based on an approach with high deterministic order, allowing for larger stepsizes and leading to fast simulations. We compare it to the Euler ?-leap, as well as two more recent ?-leap methods, on a number of example problems, and find that as well as being very accurate, our method is the most robust, in terms of efficiency, of all the methods considered in this paper. The problems it is most suited for are those with increased populations that would be too slow to simulate using Gillespie’s stochastic simulation algorithm. For such problems, it is likely to achieve higher weak order in the moments.ConclusionsThe Stochastic Bulirsch-Stoer method is a novel stochastic solver that can be used for fast and accurate simulations. Crucially, compared to other similar methods, it better retains its high accuracy when the timesteps are increased. Thus the Stochastic Bulirsch-Stoer method is both computationally efficient and robust. These are key properties for any stochastic numerical method, as they must typically run many thousands of simulations

Accurate implementation of leaping in space: The spatial partitioned-leaping algorithm

There is a great need for accurate and efficient computational approaches that can account for both the discrete and stochastic nature of chemical interactions as well as spatial inhomogeneities and diffusion. This is particularly true in biology and nanoscale materials science, where the common assumptions of deterministic dynamics and well-mixed reaction volumes often break down. In this article, we present a spatial version of the partitioned-leaping algorithm (PLA), a multiscale accelerated-stochastic simulation approach built upon the tau-leaping framework of Gillespie. We pay special attention to the details of the implementation, particularly as it pertains to the time step calculation procedure. We point out conceptual errors that have been made in this regard in prior implementations of spatial tau-leaping and illustrate the manifestation of these errors through practical examples. Finally, we discuss the fundamental difficulties associated with incorporating efficient exact-stochastic techniques, such as the next-subvolume method, into a spatial-leaping framework and suggest possible solutions.Comment: 15 pages, 9 figures, 2 table

Computationally-efficient stochastic cluster dynamics method for modeling damage accumulation in irradiated materials

An improved version of a recently developed stochastic cluster dynamics (SCD) method {[}Marian, J. and Bulatov, V. V., {\it J. Nucl. Mater.} \textbf{415} (2014) 84-95{]} is introduced as an alternative to rate theory (RT) methods for solving coupled ordinary differential equation (ODE) systems for irradiation damage simulations. SCD circumvents by design the curse of dimensionality of the variable space that renders traditional ODE-based RT approaches inefficient when handling complex defect population comprised of multiple (more than two) defect species. Several improvements introduced here enable efficient and accurate simulations of irradiated materials up to realistic (high) damage doses characteristic of next-generation nuclear systems. The first improvement is a procedure for efficiently updating the defect reaction-network and event selection in the context of a dynamically expanding reaction-network. Next is a novel implementation of the $\tau$-leaping method that speeds up SCD simulations by advancing the state of the reaction network in large time increments when appropriate. Lastly, a volume rescaling procedure is introduced to control the computational complexity of the expanding reaction-network through occasional reductions of the defect population while maintaining accurate statistics. The enhanced SCD method is then applied to model defect cluster accumulation in iron thin films subjected to triple ion-beam ($\text{Fe}^{3+}$, $\text{He}^{+}$ and $$\text{H\ensuremath{{}^{+}}}$$) irradiations, for which standard RT or spatially-resolved kinetic Monte Carlo simulations are prohibitively expensive

Simulation and inference algorithms for stochastic biochemical reaction networks: from basic concepts to state-of-the-art

Stochasticity is a key characteristic of intracellular processes such as gene regulation and chemical signalling. Therefore, characterising stochastic effects in biochemical systems is essential to understand the complex dynamics of living things. Mathematical idealisations of biochemically reacting systems must be able to capture stochastic phenomena. While robust theory exists to describe such stochastic models, the computational challenges in exploring these models can be a significant burden in practice since realistic models are analytically intractable. Determining the expected behaviour and variability of a stochastic biochemical reaction network requires many probabilistic simulations of its evolution. Using a biochemical reaction network model to assist in the interpretation of time course data from a biological experiment is an even greater challenge due to the intractability of the likelihood function for determining observation probabilities. These computational challenges have been subjects of active research for over four decades. In this review, we present an accessible discussion of the major historical developments and state-of-the-art computational techniques relevant to simulation and inference problems for stochastic biochemical reaction network models. Detailed algorithms for particularly important methods are described and complemented with MATLAB implementations. As a result, this review provides a practical and accessible introduction to computational methods for stochastic models within the life sciences community

A hybrid multiscale Monte Carlo algorithm (HyMSMC) to cope with disparity in time scales and species populations in intracellular networks

<p>Abstract</p> <p>Background</p> <p>The fundamental role that intrinsic stochasticity plays in cellular functions has been shown via numerous computational and experimental studies. In the face of such evidence, it is important that intracellular networks are simulated with stochastic algorithms that can capture molecular fluctuations. However, separation of time scales and disparity in species population, two common features of intracellular networks, make stochastic simulation of such networks computationally prohibitive. While recent work has addressed each of these challenges separately, a generic algorithm that can <it>simultaneously </it>tackle disparity in time scales <it>and </it>population scales in stochastic systems is currently lacking. In this paper, we propose the hybrid, multiscale Monte Carlo (HyMSMC) method that fills in this void.</p> <p>Results</p> <p>The proposed HyMSMC method blends stochastic singular perturbation concepts, to deal with potential stiffness, with a hybrid of exact and coarse-grained stochastic algorithms, to cope with separation in population sizes. In addition, we introduce the computational singular perturbation (CSP) method as a means of systematically partitioning fast and slow networks and computing relaxation times for convergence. We also propose a new criteria of convergence of fast networks to stochastic low-dimensional manifolds, which further accelerates the algorithm.</p> <p>Conclusion</p> <p>We use several prototype and biological examples, including a gene expression model displaying bistability, to demonstrate the efficiency, accuracy and applicability of the HyMSMC method. Bistable models serve as stringent tests for the success of multiscale MC methods and illustrate limitations of some literature methods.</p

Modelling non-Markovian dynamics in biochemical reactions

Biochemical reactions are often modelled as discrete-state continuous-time stochastic processes evolving as memoryless Markov processes. However, in some cases, biochemical systems exhibit non-Markovian dynamics. We propose here a methodology for building stochastic simulation algorithms which model more precisely non-Markovian processes in some specific situations. Our methodology is based on Constraint Programming and is implemented by using Gecode, a state-of-the-art framework for constraint solving

Stochastic modelling of eukaryotic cell cycle

Stochastic models are developed to capture the inherent stochasticity of the biochemical networks associated to many biological processes. The objective of the present thesis is to present a detailed picture of stochastic approach for the mathematical modeling of eukaryotic cell cycle, to demonstrate an important application of such model in chemotherapy and to present a methodology for selecting the model parameters. The stochastic cell cycle model, developed using stochastic chemical kinetics approach, leads to the formation of an infinite dimensional differential equation in probabilities of system being in a specific state. Using Monte Carlo simulations of this model, dynamics of populations of eukaryotic cells such as yeasts or mammalian cells are obtained. Simulations are stochastic in nature and therefore exhibit variability among cells that is similar to the variability observed in natural populations. The model’s capability to predict heterogeneities in cell populations is used as a basis to implement it in a chemotherapic modeling framework to demonstrate how the model can be used to assist in the drug development stage by investigating drug administration strategies that can have different killing effect on cancer cells and healthy cells. Finally, basic cell cycle model is refined in a systematic way to make it more suitable for describing the population characteristics of budding yeast. Selection of model parameters using an evolutionary optimization strategy referred to as insilico evolution is described. The benefits of this approach lie in the fact that it generates an initial guess of reasonable set of parameters which in turn can be used in the least squares fitting of model to the steady state distributions obtained from flow cytometry measurements. The Insilco evolution algorithm serves as a tool for sensitivity analysis of the model parameters and leads to a synergistic approach of model and experiments guiding each other. To conclude, the stochastic model based on single cell kinetics will be useful for predicting the population distribution on whole organism level. Such models find applications in wide areas of biological and biomedical applications. Evolutionary optimization strategies can be used in parameter estimation methods based on steady state distributions