Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions

Abstract

The dynamical solution of a well-mixed, nonlinear stochastic chemical kinetic system, described by the Master equation, may be exactly computed using the stochastic simulation algorithm. However, because the computational cost scales with the number of reaction occurrences, systems with one or more "fast" reactions become costly to simulate. This paper describes a hybrid stochastic method that partitions the system into subsets of fast and slow reactions, approximates the fast reactions as a continuous Markov process, using a chemical Langevin equation, and accurately describes the slow dynamics using the integral form of the "Next Reaction" variant of the stochastic simulation algorithm. The key innovation of this method is its mechanism of efficiently monitoring the occurrences of slow, discrete events while simultaneously simulating the dynamics of a continuous, stochastic or deterministic process. In addition, by introducing an approximation in which multiple slow reactions may occur within a time step of the numerical integration of the chemical Langevin equation, the hybrid stochastic method performs much faster with only a marginal decrease in accuracy. Multiple examples, including a biological pulse generator and a large-scale system benchmark, are simulated using the exact and proposed hybrid methods as well as, for comparison, a previous hybrid stochastic method. Probability distributions of the solutions are compared and the weak errors of the first two moments are computed. In general, these hybrid methods may be applied to the simulation of the dynamics of a system described by stochastic differential, ordinary differential, and Master equations.

Tumour suppression by immune system through stochastic oscillations

“…In order to simulate the model, it was necessary to define a proper stochastic simulation algorithm. Indeed, although in the literature we find algorithms for simulating time-dependent propensity functions (Lecca 2006, Anderson 2007 as well as algorithms for simulating hybrid models (Salis andKaznessis 2005, Alfonsi et al 2004), the combination of both is missing. Consequently, for the sake of analysing this model we defined an ad-hoc algorithm based on the ideas of Gillespie's Stochastic Simulation Algorithm (Gillespie 1976, Gillespie 1977 which permits an efficient implementation.…”

“…In order to analyse this model, we need to define an algorithm for computing its time evolution. Although in the literature algorithms for simulating hybrid systems have been presented (Salis andKaznessis 2005, Alfonsi et al 2004), in this model we combine the hybrid approach with time-dependent propensity functions (Lecca 2006, Anderson 2007. Consequently, we define, based on the main ideas of Gillespie's algorithm (Gillespie 1976, Gillespie 1977), a stochastic evolution algorithm for this hybrid model.…”

Hybrid modeling and simulation of stochastic effects on progression through the eukaryotic cell cycle

“…In this method, the fast reaction group is governed by ODEs or CLEs and the slow reaction group is simulated by Gillespie's direct method. A similar strategy was adopted by Salis et al, 29,30 but fast reactions are approximated by CLEs and slow reactions are simulated by Gibson and Bruck's next reaction method. 32 They also developed a more efficient mechanism to monitor the occurrences of slow, discrete events while simultaneously simulating the dynamics of a continuous, stochastic or deterministic process.…”

“…Our work follows the original idea of the hybrid method from Haseltine and Rawlings, 19 and adopts a similar implementation strategy for event handling as in Salis et al 29,30 Suppose the system has N species, and its state vector is denoted by X(t) = (X 1 (t), . .…”

“…19,29,31 Cao et al 31 proposed to partition the system, based simply on the species population numbers. For species whose population numbers are less than a threshold, all related reactions are simulated by the SSA, while other reactions are simulated by the tau-leaping method.…”