Accelerating Exact Stochastic Simulation of Biochemical Systems

The ability to accurately and efficiently simulate computer models of biochemical systems is of growing importance to the molecular biology and pharmaceutical research communities. Exact stochastic simulation is a popular approach for simulating such systems because it properly represents genetic noise and it accurately represents systems with small populations of chemical species. Unfortunately, the computational demands of exact stochastic simulation often limit its applicability. To enable next-generation whole-cell and multi-cell stochastic modeling, advanced tools and techniques must be developed to increase simulation efficiency. This work assesses the applicability of a variety of hardware and software acceleration approaches for exact stochastic simulation including serial algorithm improvements, parallel computing, reconfigurable computing, and cluster computing. Through this analysis, improved simulation techniques for biological systems are explored and evaluated

BioSimulator.jl: Stochastic simulation in Julia

Biological systems with intertwined feedback loops pose a challenge to mathematical modeling efforts. Moreover, rare events, such as mutation and extinction, complicate system dynamics. Stochastic simulation algorithms are useful in generating time-evolution trajectories for these systems because they can adequately capture the influence of random fluctuations and quantify rare events. We present a simple and flexible package, BioSimulator.jl, for implementing the Gillespie algorithm, $\tau$-leaping, and related stochastic simulation algorithms. The objective of this work is to provide scientists across domains with fast, user-friendly simulation tools. We used the high-performance programming language Julia because of its emphasis on scientific computing. Our software package implements a suite of stochastic simulation algorithms based on Markov chain theory. We provide the ability to (a) diagram Petri Nets describing interactions, (b) plot average trajectories and attached standard deviations of each participating species over time, and (c) generate frequency distributions of each species at a specified time. BioSimulator.jl's interface allows users to build models programmatically within Julia. A model is then passed to the simulate routine to generate simulation data. The built-in tools allow one to visualize results and compute summary statistics. Our examples highlight the broad applicability of our software to systems of varying complexity from ecology, systems biology, chemistry, and genetics. The user-friendly nature of BioSimulator.jl encourages the use of stochastic simulation, minimizes tedious programming efforts, and reduces errors during model specification.Comment: 27 pages, 5 figures, 3 table

Efficient Finite Difference Method for Computing Sensitivities of Biochemical Reactions

Sensitivity analysis of biochemical reactions aims at quantifying the dependence of the reaction dynamics on the reaction rates. The computation of the parameter sensitivities, however, poses many computational challenges when taking stochastic noise into account. This paper proposes a new finite difference method for efficiently computing sensitivities of biochemical reactions. We employ propensity bounds of reactions to couple the simulation of the nominal and perturbed processes. The exactness of the simulation is reserved by applying the rejection-based mechanism. For each simulation step, the nominal and perturbed processes under our coupling strategy are synchronized and often jump together, increasing their positive correlation and hence reducing the variance of the estimator. The distinctive feature of our approach in comparison with existing coupling approaches is that it only needs to maintain a single data structure storing propensity bounds of reactions during the simulation of the nominal and perturbed processes. Our approach allows to computing sensitivities of many reaction rates simultaneously. Moreover, the data structure does not require to be updated frequently, hence improving the computational cost. This feature is especially useful when applied to large reaction networks. We benchmark our method on biological reaction models to prove its applicability and efficiency.Comment: 29 pages with 6 figures, 2 table

An FPGA Based Implementation of the Exact Stochastic Simulation Algorithm

Mathematical and statistical modeling of biological systems is a desired goal for many years. Many biochemical models are often evaluated using a deterministic approach, which uses differential equations to describe the chemical interactions. However, such an approach is inaccurate for small species populations as it neglects the discrete representation of population values, presents the possibility of negative populations, and does not represent the stochastic nature of biochemical systems. The Stochastic Simulation Algorithm (SSA) developed by Gillespie is able to properly account for these inherent noise fluctuations. Due to the stochastic nature of the Monte Carlo simulations, large numbers of simulations must be run in order to get accurate statistics for the species populations and reactions. However, the algorithm tends to be computationally heavy and leads to long simulation runtimes for large systems. Therefore, this thesis explores implementing the SSA on a Field Programmable Gate Array (FPGA) to improve performance. Employing the Field programmable Gate Arrays exploits the parallelism present in the SSA, providing speedup over the software implementations that execute sequentially. In contrast to prior work that requires re-construction and re-synthesis of the design to simulate a new biochemical system, this work explores the use of reconfigurable hardware in implementing a generic biochemical simulator

Reconfigurable Hardware Acceleration of Exact Stochastic Simulation

This thesis explores the use of reconfigurable hardware in modeling chemical species reacting in a spatially homogeneous environment. The time evolution of biochemical models is often evaluated using a deterministic approach that uses differential equations to describe the chemical interactions of the model. However, such an approach treats species as continuous valued concentrations, is inaccurate for small species populations, and neglects the stochastic nature of biochemical systems. The Stochastic Simulation Algorithm (SSA) developed by Gillespie is able to properly account for these inherent noise fluctuations. This allows the SSA to accurately project the time evolution of a biochemical model. Unfortunately, the SSA can be computationally intensive and require a substantial amount of time to complete. Therefore, it has been proposed that the SSA be implemented on a Field Programmable Gate Array (FPGA) to improve performance. Employing an FPGA allows parallelism to be exploited within the SSA providing a speedup over software implementations executing instructions sequentially. Recent work in this area has focused on implementing the SSA on an FPGA to simulate specific biochemical models. However, this requires re-constructing and re-synthesizing the design in order to simulate a new biochemical system. This work examines the use of a reconfigurable computing platform to allow an implementation of the SSA on an FPGA to simulate a variety of models. The designs presented herein demonstrate a speedup of roughly 1.5X

Accelerating the Gillespie τ-Leaping Method Using Graphics Processing Units

The Gillespie τ-Leaping Method is an approximate algorithm that is faster than the exact Direct Method (DM) due to the progression of the simulation with larger time steps. However, the procedure to compute the time leap τ is quite expensive. In this paper, we explore the acceleration of the τ-Leaping Method using Graphics Processing Unit (GPUs) for ultra-large networks ( reaction channels). We have developed data structures and algorithms that take advantage of the unique hardware architecture and available libraries. Our results show that we obtain a performance gain of over 60x when compared with the best conventional implementations

Recycling random numbers in the stochastic simulation algorithm

The stochastic simulation algorithm (SSA) was introduced by Gillespie and in a different form by Kurtz. Since its original formulation there have been several attempts at improving the efficiency and hence the speed of the algorithm. We briefly discuss some of these methods before outlining our own simple improvement, the recycling direct method (RDM), and demonstrating that it is capable of increasing the speed of most stochastic simulations. The RDM involves the statistically acceptable recycling of random numbers in order to reduce the computational cost associated with their generation and is compatible with several of the pre-existing improvements on the original SSA. Our improvement is also sufficiently simple (one additional line of code) that we hope will be adopted by both trained mathematical modelers and experimentalists wishing to simulate their model systems

Uniformisation techniques for stochastic simulation of chemical reaction networks

This work considers the method of uniformisation for continuous-time Markov chains in the context of chemical reaction networks. Previous work in the literature has shown that uniformisation can be beneficial in the context of time-inhomogeneous models, such as chemical reaction networks incorporating extrinsic noise. This paper lays focus on the understanding of uniformisation from the viewpoint of sample paths of chemical reaction networks. In particular, an efficient pathwise stochastic simulation algorithm for time-homogeneous models is presented which is complexity-wise equal to Gillespie's direct method. This new approach therefore enlarges the class of problems for which the uniformisation approach forms a computationally attractive choice. Furthermore, as a new application of the uniformisation method, we provide a novel variance reduction method for (raw) moment estimators of chemical reaction networks based upon the combination of stratification and uniformisation