Adaptive hybrid simulations for multiscale stochastic reaction networks

Abstract

The probability distribution describing the state of a Stochastic Reaction Network evolves according to the Chemical Master Equation (CME). It is common to estimated its solution using Monte Carlo methods such as the Stochastic Simulation Algorithm (SSA). In many cases these simulations can take an impractical amount of computational time. Therefore many methods have been developed that approximate the Stochastic Process underlying the Chemical Master Equation. Prominent strategies are Hybrid Models that regard the firing of some reaction channels as being continuous and applying the quasi-stationary assumption to approximate the dynamics of fast subnetworks. However as the dynamics of a Stochastic Reaction Network changes with time these approximations might have to be adapted during the simulation. We develop a method that approximates the solution of a CME by automatically partitioning the reaction dynamics into discrete/continuous components and applying the quasi-stationary assumption on identifiable fast subnetworks.Our method does not require user intervention and it adapts to exploit the changing timescale separation between reactions and/or changing magnitudes of copy numbers of constituent species. We demonstrate the efficiency of the proposed method by considering examples from Systems Biology and showing that very good approximations to the exact probability distributions can be achieved in significantly less computational time.Keywords: adaptive hybrid models, stochastic reaction networks, quasi-stationary, adaptive piecewise deterministic markov processes 1

Sensitivity Analysis for Multiscale Stochastic Reaction Networks Using Hybrid Approximations

“…To accommodate these considerations a scaling parameter N 0 is chosen that corresponds to the system volume or the magnitude of the population-size of an abundant species. Thereafter each species S i is assigned an abundance factor α i ≥ 0 and each reaction k is assigned a scaling constant β k ∈ R. These parameters serve as normalizing constants, in the sense that if X i (t) denotes the copy-number of species S i at time t, then N −α i 0 X i (N γ 0 t) is roughly of order 1 or O(1) on the timescale of interest γ ∈ R. Similarly the scaled rate constant κ k = κ k N −β k 0 is O(1) for each reaction k. There exists computational approaches to automatically select these normalizing parameters α i -s and β k -s [35].…”

“…In this section we briefly discuss how the PDMP (Z(t)) t≥0 defined by (2.11) can be simulated. Many strategies exist for efficient simulation of PDMPs [25,36,35]. In this paper we shall simulate PDMPs by adapting Algorithm 2 in [36], which is a generalization of the next reaction method (NRM) [17] and is based on representation (2.11) using time-changed Poisson processes.…”

Stochastic hybrid models of gene regulatory networks – A PDE approach

“…Therefore, hybrid models have become popular in which for highly-abundant species only average counts are tracked while discrete random variables are used to represent species with low copy numbers. These hybrid approaches allow for faster and yet accurate Monte Carlo sampling that stochastically selects counts of species with low copy numbers and numerically integrates average counts of all other species [HGK15,CDR09,HH12,PK04].…”

Geometric Ergodicity for Piecewise Contracting Processes with Applications for Tropical Stochastic Lattice Models

“…In principle, piecewise contraction holds for systems that are otherwise contracting if there is no stochastic turbulence. This includes a wide range of models in engineering, biochemistry, economics, and natural science [10,19,36,38]. On the other hand, it is relatively easy to verify as it depends solely on .…”