A higher-order numerical framework for stochastic simulation of chemical reaction systems

Background: In this paper, we present a framework for improving the accuracy of fixed-step methods for Monte Carlo simulation of discrete stochastic chemical kinetics. Stochasticity is ubiquitous in many areas of cell biology, for example in gene regulation, biochemical cascades and cell-cell interaction. However most discrete stochastic simulation techniques are slow. We apply Richardson extrapolation to the moments of three fixed-step methods, the Euler, midpoint and theta-trapezoidal tau-leap methods, to demonstrate the power of stochastic extrapolation. The extrapolation framework can increase the order of convergence of any fixed-step discrete stochastic solver and is very easy to implement; the only condition for its use is knowledge of the appropriate terms of the global error expansion of the solver in terms of its stepsize. In practical terms, a higher-order method with a larger stepsize can achieve the same level of accuracy as a lower-order method with a smaller one, potentially reducing the computational time of the system. Results: By obtaining a global error expansion for a general weak first-order method, we prove that extrapolation can increase the weak order of convergence for the moments of the Euler and the midpoint tau-leap methods, from one to two. This is supported by numerical simulations of several chemical systems of biological importance using the Euler, midpoint and theta-trapezoidal tau-leap methods. In almost all cases, extrapolation results in an improvement of accuracy. As in the case of ordinary and stochastic differential equations, extrapolation can be repeated to obtain even higher-order approximations. Conclusions: Extrapolation is a general framework for increasing the order of accuracy of any fixed-step stochastic solver. This enables the simulation of complicated systems in less time, allowing for more realistic biochemical problems to be solved.

Propagation Models for Biochemical Reaction Networks

Maria-Emanuela-Canini Mateescu

In this thesis we investigate different ways of approximating the solution of the chemical master equation (CME). The CME is a system of differential equations that models the stochastic transient behaviour of biochemical reaction networks. It does so by describing the time evolution of probability distribution over the states of a Markov chain that represents a biological network, and thus its stochasticity is only implicit. The transient solution of a CME is the vector of probabilities over the states of the corresponding Markov chain at a certain time point t, and it has traditionally been obtained by applying methods that are general to continuous-time Markov chains: uniformization, Krylov subspace methods, and general ordinary differential equation (ODE) solvers such as the fourth order Runge-Kutta method. Even though biochemical reaction networks are the main application of our work, some of our results are presented in the more general framework of propagation models (PM), a computational formalism that we introduce in the first part of this thesis. Each propagation model N has two associated propagation processes, one in discrete-time and a second one in continuous-time. These propagation processes propagate a generic mass through a discrete state space. For example, in order to model a CME, N propagates probability mass. In the discrete-time case the propagation is done step-wise, while in the continuous-time case it is done in a continuous flow defined by a differential equation. Again, in the case of the chemical master equation, this differential equation is the equivalent of the chemical master equation itself where probability mass is propagated through a discrete state space. Discrete-time propagation processes can encode methods such as the uniformization method and the fourth order Runge-Kutta integration method that we have mentioned above, and thus by optimizing propagation algorithms we optimize both of these methods simultaneously. In the second part of our thesis, we define stochastic hybrid models that approximate the stochastic behaviour of biochemical reaction networks by treating some variables of the system deterministically. This deterministic approximation is done for species with large populations, for which stochasticity does not play an important role. We propose three such hybrid models, which we introduce from the coarsest to the most refined one: (i) the first one replaces some variables of the system with their overall expectations, (ii) the second one replaces some variables of the system with their expectations conditioned on the values of the stochastic variables, (iii) and finally, the third one, splits each variable into a stochastic part (for low valuations) and a deterministic part (for high valuations), while tracking the conditional expectation of the deterministic part. For each of these algorithms we give the corresponding propagation models that propagate not only probabilities but also the respective continuous approximations for the deterministic variables.

EPFL2011

Scalable Simulation Methodologies for Many-Core Heterogeneous Systems

Shivani Raghav

With increasing complexity and performance demands of emerging compute-intensive data-parallel workloads, many-core computing systems are becoming a popular trend in computer design. Fast and scalable simulation methods are needed to make meaningful predictions of design alternatives to prepare early software development and to assess the performance of a system before the real hardware is available. However, simulation technology for large-scale, many-core systems is lagging behind. Most of the existing simulation techniques are slow and complex or have poor scalability and high cost of development, which leads to an unacceptable performance when simulating a large number of cores. New challenges brought by emerging many-core systems demand new methods for simulating these platforms. This thesis investigates the techniques for improving the performance of parallel simulators to model many-core heterogeneous architectures. With the recent advances in programmability and processing power of graphical hardware, General Purpose Graphical Processing Units (GPGPUs) are becoming a popular host platform for solving broad range of computationally expensive problems. The easy availability of low cost, highly parallel GPGPU platforms presents an opportunity for accelerating architectural simulations. In this thesis, I propose and investigate a novel idea of using GPGPUs as a host platform for accelerating simulation of many-core systems. GPGPUs are most suitable for certain specific highly parallel workloads with fine-grained multi-threading and high computational intensity. Architectural simulation, on the other hand, is an extremely complex task with numerous dependencies and rules. To this end, I first determine the feasibility of accelerating various target many-core architectures at different levels of abstraction details, considering the particularities of prevalent GPU architectures and their programming models. Based on this study, I target a heterogeneous architecture of a multi-core CPU connected to many-core coprocessors for GPGPU acceleration. I propose ways to exploit ample parallelism in many-core simulation and present a comprehensive framework to develop an architecture-independent, fast, scalable, easy-to-use, full-system simulator. In particular, I outline some of the challenges faced and insights gained using our proposed approach. I present specific methods and approaches that maximize the throughput of GPU power and memory bandwidth allowing optimized parallelization with minimal synchronization. Our results show high scalability and high speed up in performance for simulating up to eight thousand cores.

EPFL2014

Explicit Stabilized Multirate Method For Stiff Stochastic Differential Equations

Assyr Abdulle, Giacomo Rosilho De Souza

Stabilized explicit methods are particularly efficient, for large systems of stiff stochastic differential equations (SDEs) due to their extended stability domain. However, they lose their efficiency when a severe stiffness is induced by very few "fast" degrees of freedom, as the stiff and nonstiff terms are evaluated concurrently. Therefore, inspired by [A. Abdulle, M. J. Grote, and G. Rosilho de Souza, Explicit stabilized multirate method for stiff differential equations, Math. Comp., in press, 2022], we introduce a stochastic modified equation whose stiffness depends solely on vi the "slow" terms. By integrating this modified equation with a stabilized explicit scheme, we devise a multirate method which overcomes the bottleneck caused by a few severely still terms and recovers the efficiency of stabilized schemes for large systems of nonlinear SDEs. The scheme is not based on any scale separation assumption of the SDE. Therefore, it is employable for problems stemming from the spatial discretization of stochastic parabolic partial differential equations on locally refined grids. The multirate scheme has strong Order 1/2, weak order 1, and its stability is proved on a model problem. Numerical experiments confirm the efficiency and accuracy of the scheme.