Scaling methods for accelerating kinetic Monte Carlo simulations of chemical reaction networks

TL;DR

This paper introduces a novel adaptive and heterogeneous scaling algorithm called partial scaling to accelerate kinetic Monte Carlo simulations of chemical reaction networks, reducing computational cost while maintaining accuracy.

Contribution

The paper presents a new partial scaling method that adaptively adjusts reaction rates without prior classification of species, implemented in BioNetGen, enhancing simulation efficiency.

Findings

Significantly accelerates biological system models

Maintains accuracy with reduced reaction events

Applicable to diverse biological models

Abstract

Various kinetic Monte Carlo algorithms become inefficient when some of the population sizes in a system are large, which gives rise to a large number of reaction events per unit time. Here, we present a new acceleration algorithm based on adaptive and heterogeneous scaling of reaction rates and stoichiometric coefficients. The algorithm is conceptually related to the commonly used idea of accelerating a stochastic simulation by considering a sub-volume $\lambda \Omega$ ($0<\lambda<1$) within a system of interest, which reduces the number of reaction events per unit time occurring in a simulation by a factor $1/\lambda$ at the cost of greater error in unbiased estimates of first moments and biased overestimates of second moments. Our new approach offers two unique benefits. First, scaling is adaptive and heterogeneous, which eliminates the pitfall of overaggressive scaling. Second, there is no need for an \emph{a priori} classification of populations as discrete or continuous (as in a hybrid method), which is problematic when discreteness of a chemical species changes during a simulation. The method requires specification of only a single algorithmic parameter, $N_c$, a global critical population size above which populations are effectively scaled down to increase simulation efficiency. The method, which we term partial scaling, is implemented in the open-source BioNetGen software package. We demonstrate that partial scaling can significantly accelerate simulations without significant loss of accuracy for several published models of biological systems. These models characterize activation of the mitogen-activated protein kinase ERK, prion protein aggregation, and T-cell receptor signaling.

Figures

Tables

Table 1. Table 1: Efficiency and accuracy of standard scaling and partial scaling methods for stochastic simulation based on models for ERK activation Kochańczyk  et al. ( 2017 ) , prion protein aggregation Rubenstein  et al. ( 2007 ) , and TCR signaling Lipniacki  et al. ( 2008 ) . CPU time was measured using the C++ clock() function. Summary statistics for first-moment errors, ⟨ ε 1 ⟩ delimited-⟨⟩ subscript 𝜀 1 \left\langle\varepsilon_{1}\right\rangle , and Van Valen coefficients of variation, ⟨ C ​ V VV ⟩ delimited-⟨⟩ 𝐶 subscript 𝑉 VV \left\langle CV_{\text{VV}}\right\rangle , were calculated as described in the Appendix.

Model	Scaling method	CPU time	$⟨{\epsilon }_1⟩$	$⟨C{V}_{\text{VV}}⟩$
ERK	None (Gillespie SSA)	$5.30\times {10}^2\pm 1.17\times {10}^1$	$0.00\times {10}^0$	$2.98\times {10}^{-2}$
ERK	Standard ($\lambda =0.3$)	$1.57\times {10}^2\pm 6.10\times {10}^0$	$4.36\times {10}^{-1}$	$5.39\times {10}^{-2}$
ERK	Standard ($\lambda =0.1$)	$5.21\times {10}^1\pm 2.42\times {10}^0$	$1.31\times {10}^0$	$9.44\times {10}^{-2}$
ERK	Standard ($\lambda =0.03$)	$1.52\times {10}^1\pm 1.10\times {10}^0$	$3.99\times {10}^{+2}$	$1.77\times {10}^{-1}$
ERK	Standard ($\lambda =0.01$)	$4.63\times {10}^0\pm 5.12\times {10}^{-1}$	$1.20\times {10}^{+3}$	$2.79\times {10}^{-1}$
ERK	Partial ($N_c=3.0\times {10}^4$)	$4.83\times {10}^2\pm 1.01\times {10}^1$	$3.44\times {10}^{-3}$	$3.03\times {10}^{-2}$
ERK	Partial ($N_c=1.0\times {10}^4$)	$4.14\times {10}^2\pm 9.25\times {10}^0$	$1.75\times {10}^{-3}$	$3.17\times {10}^{-2}$
ERK	Partial ($N_c=3.0\times {10}^3$)	$2.82\times {10}^2\pm 8.42\times {10}^0$	$1.47\times {10}^{-2}$	$3.49\times {10}^{-2}$
ERK	Partial ($N_c=1.0\times {10}^3$)	$9.95\times {10}^1\pm 4.42\times {10}^0$	$1.67\times {10}^{-3}$	$4.50\times {10}^{-2}$
ERK	Partial ($N_c=3.0\times {10}^2$)	$4.47\times {10}^1\pm 2.48\times {10}^0$	$4.80\times {10}^{-1}$	$6.94\times {10}^{-2}$
Prion	None (Gillespie SSA)	$4.39\times {10}^2\pm 5.24\times {10}^1$	$0.00\times {10}^0$	$2.03\times 10-2$
Prion	Standard ($\lambda =0.3$)	$1.25\times {10}^2\pm 1.38\times {10}^1$	$1.22\times {10}^{-1}$	$3.00\times 10-1$
Prion	Standard ($\lambda =0.1$)	$2.61\times {10}^1\pm 7.61\times {10}^{-1}$	$2.87\times {10}^{-1}$	$3.92\times 10-2$
Prion	Standard ($\lambda =0.03$)	$8.19\times {10}^0\pm 4.90\times {10}^{-1}$	$2.35\times {10}^{-1}$	$5.73\times 10-2$
Prion	Standard ($\lambda =0.01$)	$3.08\times {10}^0\pm 2.68\times {10}^{-1}$	$2.39\times {10}^{-1}$	$5.96\times 10-2$
Prion	Partial ($N_c=3.0\times {10}^3$)	$6.29\times {10}^2\pm 6.51\times {10}^0$	$1.34\times {10}^{-4}$	$2.18\times 10-2$
Prion	Partial ($N_c=1.0\times {10}^3$)	$4.50\times {10}^2\pm 5.67\times {10}^0$	$1.33\times {10}^{-4}$	$2.79\times 10-2$
Prion	Partial ($N_c=3.0\times {10}^2$)	$4.18\times {10}^2\pm 5.56\times {10}^0$	$1.32\times {10}^{-4}$	$3.95\times 10-2$
Prion	Partial ($N_c=1.0\times {10}^2$)	$3.35\times {10}^2\pm 1.58\times {10}^1$	$1.93\times {10}^{-4}$	$7.84\times 10-2$
Prion	Partial ($N_c=3.0\times {10}^1$)	$1.69\times {10}^2\pm 3.34\times {10}^0$	$1.81\times {10}^{-4}$	$1.08\times 10-1$
Prion	Partial ($N_c=1.0\times {10}^1$)	$4.76\times {10}^1\pm 3.57\times {10}^0$	$5.44\times {10}^{-4}$	$1.94\times 10-1$
TCR	None (Gillespie SSA)	$1.04\times {10}^2\pm 2.27\times {10}^1$	$0.00\times {10}^0$	$3.06\times {10}^{-1}$
TCR	Standard ($\lambda =0.3$)	$2.06\times {10}^1\pm 3.82\times {10}^0$	$2.82\times {10}^{-1}$	$2.75\times {10}^{-1}$
TCR	Standard ($\lambda =0.1$)	$4.86\times {10}^0\pm 9.44\times {10}^{-1}$	$4.40\times {10}^{-1}$	$2.29\times {10}^{-1}$
TCR	Partial ($N_c=1.0\times {10}^3$)	$3.62\times {10}^1\pm 5.02\times {10}^0$	$1.19\times {10}^{-3}$	$3.16\times {10}^{-1}$
TCR	Partial ($N_c=3.0\times {10}^2$)	$1.63\times {10}^1\pm 1.58\times {10}^0$	$1.20\times {10}^{-3}$	$2.99\times {10}^{-1}$
TCR	Partial ($N_c=1.0\times {10}^2$)	$8.48\times {10}^0\pm 8.25\times {10}^{-1}$	$9.46\times {10}^{-4}$	$3.10\times {10}^{-1}$
TCR	Partial ($N_c=3.0\times {10}^1$)	$3.23\times {10}^0\pm 5.42\times {10}^{-1}$	$1.63\times {10}^{-3}$	$3.05\times {10}^{-1}$

Equations

$$\varnothing \stackrel{{\kappa }_1}{\to }$$

$$X_1+X_1\stackrel{{\kappa }_2}{\to }$$

$$X_2\stackrel{{\kappa }_3}{\to }$$

$$\frac{\text{d}{n}_1}{\text{d}t}=$$

$$\frac{\text{d}{n}_2}{\text{d}t}=$$

$$(N_1,N_2)\stackrel{h_1(\mathbf{N})\times {\bar{\kappa }}_1}{\to }$$

$$(N_1,N_2)\stackrel{h_2(\mathbf{N})\times {\bar{\kappa }}_2}{\to }$$

$$(N_1,N_2)\stackrel{h_3(\mathbf{N})\times {\bar{\kappa }}_3}{\to }$$

$$\frac{\text{d}}{\text{d}t}{\mathbb{P}}_{N_1,N_2}=$$

$$-\frac{{\kappa }_2}{\Omega }\frac{N_1\left(N_1-1\right)}{2}{\mathbb{P}}_{N_1,N_2}$$

$$+\frac{{\kappa }_2}{\Omega }\frac{\left(N_1+2\right)\left(N_1+1\right)}{2}{\mathbb{P}}_{N_1+2,N_2-1}$$

$$-{\kappa }_3\left[N_2{\mathbb{P}}_{N_1,N_2}-\left(N_2+1\right){\mathbb{P}}_{N_1,N_2+1}\right].$$

$${\partial }_t\rho \left(n_1,n_2,t\right)=-\Omega {\kappa }_1\left[\rho \left(n_1,n_2,t\right)-\rho \left(n_1-\text{d}n,n_2,t\right)\right]$$

$$-{\kappa }_2\Omega \frac{n_1\left(n_1-\text{d}n\right)}{2}\rho \left(n_1,n_2,t\right)$$

$$+{\kappa }_2\Omega \frac{\left(n_1+2\text{d}n\right)\left(n_1+\text{d}n\right)}{2}\rho \left(n_1+2\text{d}n,n_2-\text{d}n,t\right)$$

$$-{\kappa }_3\Omega \left[n_2\rho \left(n_1,n_2,t\right)-\left(n_2+\text{d}n\right)\rho \left(n_1,n_2+\text{d}n,t\right)\right].$$

$${\partial }_t$$

$$+\frac{1}{2\Omega }\{{\partial }_{n_1}^2\left[\left({\kappa }_1+{\kappa }_2{n}_1^2\right)\rho \right]+{\partial }_{n_2}^2\left[\left(\frac{{\kappa }_2}{2}{n}_1^2+{\kappa }_3{n}_2\right)\rho \right]$$

$$+2{\partial }_{n_1}{\partial }_{n_2}\left[\left({\kappa }_2{n}_1^2\right)\rho \right]\}.$$

$${\partial }_t\rho =-{\partial }_{n_1}\left[\left({\kappa }_1-{\kappa }_2{n}_1^2\right)\rho \right]-{\partial }_{n_2}\left[\left(\frac{{\kappa }_2{n}_1^2}{2}-{\kappa }_3{n}_2\right)\rho \right].$$

$${\partial }_t\rho =-{\partial }_{n_1}\left[\left({\kappa }_1-{\kappa }_2{n}_1^2\right)\rho \right]-{\partial }_{n_2}\left[\left(\frac{{\kappa }_2}{2}{n}_1^2-{\kappa }_3{n}_2\right)\rho \right]$$

$$+\frac{1}{2\Omega \lambda }\{{\partial }_{n_1}^2\left[\left({\kappa }_1+{\kappa }_2{n}_1^2\right)\rho \right]+{\partial }_{n_2}^2\left[\left(\frac{{\kappa }_2}{2}{n}_1^2+{\kappa }_3{n}_2\right)\rho \right]$$

$$+2{\partial }_{n_1}{\partial }_{n_2}\left[\left({\kappa }_2{n}_1^2\right)\rho \right]\}.$$

$${\lambda }_r\left(t\right)=\frac{1}{\max \left\{1,\lfloor \frac{N_{\min }^r\left(t\right)}{N_c}\rfloor \right\}},$$

$$\delta N_i=\sum _{r\in I}\delta N_{i,r}.$$

$$\delta N_i^{\text{scaled}}=\sum _{r\in I^{\prime }}\frac{\delta N_{i,r}}{{\lambda }_r}+\sum _{r\in I\backslash I^{\prime }}\delta N_{i,r}.$$

$$\delta N_{i,r}\sim$$

$$\delta N_{i,r}^{\star }\sim$$

$$\mathbb{E}\left[\delta N_{i,r}^{\star }\right]=$$

$$\text{var}\left[\delta N_{i,r}^{\star }\right]=$$

$$\Delta \overline{{\mathbf{N}}^{(\lambda ,N_c)}}\left(t\right)\equiv \overline{{\mathbf{N}}^{(\lambda ,N_c)}}\left(t\right)-\overline{{\mathbf{N}}^{(0)}}\left(t\right).$$

$$\varepsilon^{i}\left(t\right)\equiv\left\{\begin{array}[]{l}\frac{\bar{N}_{i}^{\left(\lambda,\bar{N}_{c}\right)}\left(t\right)-\bar{N}_{i}^{\left(0\right)}\left(t\right)}{\bar{N}_{i}^{\left(0\right)}\left(t\right)}\text{, if }\bar{N}_{i}\left(t\right)\geq 1,\\ \bar{N}_{i}^{\left(\lambda,\bar{N}_{c}\right)}\left(t\right)-\bar{N}_{i}^{\left(0\right)}\left(t\right)\text{, otherwise}\end{array}\right.$$

$$⟨{\epsilon }^i⟩=\frac{\Delta t}{T}\sum _{k=0}^{T/\Delta t}{\epsilon }^i\left(k\Delta t\right)$$

$$⟨\epsilon ⟩\equiv \sum _{i=1}^M⟨{\epsilon }^i⟩$$

$$C{V}_{\text{VV}}\left(t\right)\equiv \sqrt{\frac{\text{tr}\left(\Sigma \left(t\right)\right)}{{\left||\bar{\mathbf{N}}\left(t\right)|\right|}_2^2}},$$

$$⟨C{V}_{\text{VV}}⟩=\frac{\Delta t}{T}\sum _{k=0}^{T/\Delta t}C{V}_{\text{VV}}\left(k\Delta t\right).$$

Full text

I Introduction

Kinetic Monte Carlo (KMC) procedures Voter (2007), such as the well-known direct and next-reaction methods of GillespieGillespie (2007), are commonly used to study stochastic chemical kinetics, especially in biochemical systems. These procedures have been widely used to study systems in which population fluctuations arise from intrinsic molecular noise and the discreteness of the chemical species of interest McAdams and Arkin (1999). Population fluctuations are characteristic of many genetic regulatory circuits, as populations in these systems tend to be discrete.

KMC procedures are also useful for studying systems, such as cell signaling networks, that have large state spaces Suderman et al. (2018) arising from the combinatorial number of chemical species that can be generated by biomolecular interactions of interest Chylek et al. (2014a). For such systems, it may be impracticable, even with the aid of a computer, to enumerate the chemical species that are potentially populated. Nonetheless, simulations can be performed by formulating rules to represent biomolecular interactionsChylek et al. (2014a) and then using these rules as event generators in a so-called network-free simulation algorithm Suderman et al. (2018), such as that implemented in the NFsim software package Yang et al. (2008); Sneddon, Faeder, and Emonet (2011). In a network-free simulation (of cell signaling dynamics), stochastic effects may be negligible Creamer et al. (2012). Populations in cell signaling networks tend to be large, such that population densities are nearly continuous variables.

In any exact KMC procedure, system state is updated (and time is advanced) only one reaction event at a time, which can be a severe limitation when the number of events per unit time is large. In many applications of KMC aimed at studying noisy system behavior, this difficulty is not encountered because populations of interest are small, which makes for efficient simulation. However, inefficiency may arise if some population sizes are large and/or some reactions are fast. Both of these features of a reaction system can introduce a large, unmanageable number of reaction events per unit time. A variety of approaches have therefore been developed for addressing this problem.

In the near-continuum limit where populations of chemical species are uniformly large (but not infinite, such that concentrations can be appropriately treated as continuous variables), the diffusive stochastic differential equation (SDE) obtained from the diffusion approximation Van Kampen (1992); Risken (1996) provides an efficient solution for simulating stochastic chemical kinetics. Gillespie’s $\tau$-leaping methodGillespie (2001) provides an efficient solution for another regime of behavior, where the discreteness of populations is still relevant. In this method, a time window $\tau$ is prescribed. Importantly, $\tau$ should be (1) large enough such that multiple events occur within this window and (2) small enough such that the system’s configuration and state transition rates (also called reaction propensity functions) do not change significantly as a result of the events occurring within the window. Unfortunately, these opposing requirements limit applicability of the method. Finally, hybrid methods have been developed and applied to solve specific problemsSalis and Kaznessis (2005); Newby, Bressloff, and Keener (2013); Bokes et al. (2013); Bressloff and Newby (2014); Bressloff (2015); Lin and Doering (2016); Lin and Galla (2016); Bressloff (2017a, b); Lin and Buchler (2018). A characteristic feature of these methods is the division of state variables into two distinct sets, the discrete variables and the continuous variables. The continuous variables are evolved forward in time via numerical integration of ordinary differential equations (ODEs) while the discrete variables are simultaneously evolved forward in time via KMC. A general framework for hybrid simulation is provided by the formalism of discretely switching Langevin dynamicsMao and Yuan (2006).

All of the above-mentioned methods are designed to trade accuracy for efficiency. In other words, acceleration gains are obtained through simplifying assumptions that introduce approximations. The diffusive SDE approach is based on the assumption that only the first and second moments of distributions are important; higher order moments are ignored. The $\tau$-leaping approach is based on the assumption of Poissonian statistics in the time window $\tau$. Hybrid approaches are based on the assumption that fluctuations in certain populations can be ignored entirely over the full time window of simulation. For most, if not all, of these methods, there is no way to determine a priori if their approximations are tolerable, and unacceptable errors can arise unexpectedly during the course of a simulation. SDE and hybrid approaches become inappropriate if populations of all or a subset of chemical species, respectively, are not sufficiently large throughout a simulation. As we will see, for some systems, population sizes are distributed across multiple scales, and the size of a population can change qualitatively (e.g., from large to discrete to large again) during the course of a simulation, as when a system exhibits oscillatory behavior. The $\tau$-leaping method is also problematic: for a given system, there is no guarantee that a time window $\tau$ having the necessary properties exists. For all of these methods, approximation accuracy depends in a non-obvious way on parameters (initial conditions and rate constants). Thus, their use in any analysis involving variation of parameter values, such as a curve fitting procedure, requires careful consideration.

An alternative, simpler idea for accelerating stochastic simulation is to consider only a sub-volume of the system of interest, taking this sub-volume to be representative of the whole system. Given a system with volume $\Omega$, if we simulate reactions in only a sub-volume $\lambda\Omega$, where $0<\lambda<1$, we can expect to reduce simulation costs by the factor $1/\lambda$ because the number of reaction events per unit time is reduced by this factor. As we will discuss later in detail, this homogeneous and static scaling approach can speed up estimation of first moments of the stochastic dynamics, and estimates are unbiased so long as populations are not scaled too aggressively to size of order 1 or below. A drawback of scaling is that, with achievement of any acceleration, there is an unavoidable bias in estimation of higher moments. Another drawback is that the acceleration attainable with reasonable accuracy may be limited by the small population size of one or more critical chemical species. If system behavior is influenced by a chemical species with a small population size, as in stochastic switchingKepler and Elston (2001), the setting for $\lambda$ cannot be such that the population of the chemical species falls near or below 1, which would fundamentally change or eliminate its influence on behavior.

Scaling is commonly used as part of a (network-free) stochastic simulation approach when the dynamics of interest are nearly deterministic, as recommended, for example, by Faeder et al.Faeder, Blinov, and Hlavacek (2009). In such cases, the goal of simulation is simply to calculate first moments. The inability of scaling to yield unbiased estimates of higher moments is of little concern. However, in these applications, there remains a need to avoid overaggressive scaling of critical populations. As we will see through examples, this constraint can severely limit the usefulness of scaling via the standard (homogeneous and static) approach.

Here, we present a new approach for accelerating stochastic simulation through scaling, which we term partial scaling. In partial scaling, for each enumerated individual reaction considered in a system, its reaction rate is scaled by a factor that effectively makes the smallest reactant or product population, if greater than a critical population size $N_{c}$, equivalent or close to this critical population size. Stoichiometric coefficients are scaled by the inverse of the scaling factor for reaction rate. No scaling is applied if the smallest reactant or product population is already at or below $N_{c}$. The scaling is heterogeneous (vs. homogeneous) because each reaction has its own scaling factor, which is determined by $N_{c}$ and the smallest population of the chemical species involved in the reaction. Furthermore, reaction rates are scaled on the fly as population sizes change over time, i.e., scaling is adaptive (vs. static). Because scaling is not performed for a reaction when the smallest population of the participating chemical species is smaller than $N_{c}$, partial scaling permits not only unbiased estimation of first moments but also less biased estimation of second moments compared to the standard scaling method. Although partial scaling entails heterogeneous and adaptive scaling, whereas standard scaling entails homogeneous and static scaling, both methods involve only a single algorithmic parameter: $\lambda$ in the case of standard scaling and $N_{c}$ in the case of partial scaling.

The remainder of this report is organized as follows. In Sec. II, we introduce notation. In Sec. III, we review how scaling is currently used to accelerate stochastic simulation. Through analysis of a simple reaction network, we explain why this approach yields unbiased estimates of first moments while improving simulation efficiency at the expense of biased estimates of higher moments (e.g., overestimates of variance). We also illustrate the consequences of overaggressive scaling. In Sec. IV, to address the limitations of standard scaling, we introduce partial scaling. To complement the description of the method, we provide pseudocode for a partial scaling algorithm. In Sec. V, using published models for three complex biochemical reaction networksKochańczyk et al. (2017); Rubenstein et al. (2007); Lipniacki et al. (2008), we evaluate the performance of partial scaling relative to standard scaling. In Sec. VI, we call attention to our general-purpose implementation of partial scaling in the BioNetGen software packageHarris et al. (2016). We conclude with a discussion of results and possible future directions in Sec. VII.

II Notation

We consider a closed, well-mixed, and isothermal reaction system having constant volume111The parameter $\Omega$ characterizes system size. Here, we take it to represent volume but it can alternatively be interpreted as a population that defines a population scale. With this interpretation, (dimensionless) population densities replace concentrations as the state variables in the continuum limit with no change in the mathematical form of the governing equations. $\Omega$. The system contains a (dilute) mixture of up to $M$ chemical species, $X_{i},\ldots,X_{M}$, that react within a network of $R$ reactions. We assume that each reaction $r\in\{1,\ldots,R\}$ obeys mass-action kinetics with rate constant $\kappa_{r}$. We use $N_{i}$ to denote the non-negative and discrete population of $X_{i}$. The discrete state of the system is defined by the vector $\mathbf{N}\equiv\left(N_{1},\ldots,N_{M}\right)$. We use $\mathbb{P}_{\mathbf{N}}\left(t\right)$ to denote the probability that the system is in state $\mathbf{N}$ at time $t$. In the continuum limit, reached as $\Omega\rightarrow\infty$ and $N_{i}\rightarrow\infty$ while $N_{i}/\Omega$ remains constant for all $i$, the concentration of each $X_{i}$, defined as $n_{i}\equiv N_{i}/\Omega$, is a continuous variable, and system state is defined by the vector $\mathbf{n}\equiv(n_{1},\ldots,n_{M})$. We use $\rho\left(\mathbf{n},t\right)$ to denote the probability density of state $\mathbf{n}$ at time $t$.

III Static and homogeneous scaling

In this section, which can be bypassed by readers who are expert in stochastic modeling, we provide a brief review of scaling as a tool to accelerate stochastic simulation, mainly through the exercise of analyzing a simple reaction system. Our goal is to introduce background facts important for appreciating the benefits and limitations of scaling, which should strictly be considered only for a system with large populations, such that (all) concentrations are nearly continuous. As we will see later, many biochemical systems have concentrations that are nearly continuous. The background facts of concern here are as follows. First, the number of events generated per unit time in an exact stochastic simulation is an extensive quantity. Accordingly, the efficiency of an exact stochastic simulation algorithm (SSA) scales with system size. Second, if after scaling, the system of interest remains near the continuum limit, exact simulations of the scaled system yield unbiased estimates of the first moments of the random variables that are being sampled in the simulations. Third, scaling can introduce discreteness in populations that is inappropriate. In other words, overaggressive scaling can introduce systematic errors.

An exact SSA generates sample paths of the stochastic dynamics of a (well-mixed) reaction system by executing a series of reaction events (in a Monte Carlo fashion) that each brings about discrete changes of the populations of the reacting chemical species. State transition rates are characterized by propensity functions that, by construction, are extensive quantities Gillespie (1977); Van Kampen (1992). Thus, as we will see shortly for a specific example, the values of these functions scale linearly with system size $\Omega$, and as $\Omega$ increases, an SSA generates more events per unit time. Events occur at an overall rate proportional to $\Omega$ and inversely proportional to the average waiting time between events. Consequently, as $\Omega$ increases, execution of an SSA eventually becomes too computationally expensive to be practical.

When system size is large, the joint probability distribution of the discrete-state random process sampled by an SSA can be approximated by various methods. Ordered by degree of granularity of approximation, these methods include (1) $\tau$-leaping, which approximates the state-dependent transition kernel of an SSA by a constant (Poissonian) transition kernel, which is used over a small time window $\tau$Gillespie (2001); Cao, Gillespie, and Petzold (2005, 2006); (2) the diffusion approximation, which coarse-grains the discrete state space into a continuous state space via the Kramers–Moyal expansion, the end result of which is a Fokker–Planck equation describing the joint probability density of the system’s state characterized in terms of intensive state variables, such as concentrations or population densitiesKramers (1940); Moyal (1949); Van Kampen (1992); Gardiner et al. (1985); Doering, Sargsyan, and Sander (2005); Lin, Kim, and Doering (2012, 2012, 2015a, 2015b); and (3) the linear-noise approximation (LNA), which linearizes the multiplicative noise kernel of a Fokker–Planck equation, such that the resulting dynamics involves only Gaussian white noiseKampen (1961); Van Kampen (1992).

Another approach for coping with large system size is to use an SSA to simulate the stochastic dynamics not in the whole system of interest but rather in a suitably small sub-volume $\lambda\Omega$, where $0<\lambda<1$. The dynamics of the sub-volume are taken to represent the dynamics of the whole system. This approach, which we will refer to as static and homogeneous scaling or standard scaling, is seemingly justified for a well-mixed system, and from our discussion above, it can be expected to reduce the cost of simulation by the factor $1/\lambda$: the smaller the value of $\lambda$, the greater the efficiency gain. However, the approach can produce erroneous and misleading results if scaling is overaggressive, as when scaling reduces the population of a critical chemical species (i.e., one that influences overall system behavior even if its population size is relatively small) to 0 or $\mathcal{O}(1)$. In these cases, the dynamics in the sub-volume will differ qualitatively from the dynamics in the whole system. For example, scaling to a population size $\sim 1$ will introduce bursty behavior, which has unique statistical propertiesLin and Doering (2016); Lin and Galla (2016); Lin and Buchler (2018). Thus, scaling should be limited, but it is not clear how to impose an appropriate lower bound on $\lambda$. A constraint such as $\lambda\geq\min_{1\leq i\leq M}N_{i}\left(t\right)$ requires a trial-and-error procedure to select $\lambda$ because the bounding term is dependent on both time and parameters (initial conditions and rate constants). Another concern, even with an appropriate setting for $\lambda$, is how scaling affects the statistics of sample paths. As we will see, estimates for first moments are unbiased but scaling leads to biased overestimates of second moments.

III.1 A simple reaction network

To illustrate features of static and homogeneous scaling, we will consider the following reaction network, consisting of zeroth-, first-, and second-order elementary reactions:

[TABLE]

Here, $M=2$ and $R=3$. The chemical species $X_{1}$ is injected via a zeroth-order reaction into the system (volume $\Omega$) with rate constant $\kappa_{1}$. In a second-order reaction, two copies of $X_{1}$ react irreversibly to form $X_{2}$ with rate constant $\kappa_{2}$. The chemical species $X_{2}$ is removed from the system via a first-order reaction with rate constant $\kappa_{3}$.

III.2 Deterministic and stochastic chemical kinetics

For the reaction network of Eqs. (1), we can write the following ODEs for mass-action kinetics:

[TABLE]

where $n_{1}(t)$ is the time-dependent concentration of $X_{1}$ and $n_{2}(t)$ is the time-dependent concentration of $X_{2}$. It should be noted that the $1/2$ factor in Eq. 2b accounts for the symmetry of the left-hand side of Eq. (1b)222The symmetry factor is not strictly required in the context of ODE modeling of chemical kinetics (because a constant times a constant is still a constant), but it emerges from the CME description of chemical kinetics. and that all terms on the right-hand side of Eqs. (2) incorporate the stoichiometric coefficients of the reactions in Eqs. (1). The terms $\kappa_{1}$, $\kappa_{1}n_{1}^{2}$, $(1/2)\kappa_{2}n_{1}^{2}$, and $\kappa_{3}n_{2}$ are the mass-action rate laws that follow from reactions $r=1,2,3$ in Eqs. (1); these terms give the rates of $X_{1}$ injection, $X_{1}$ consumption, $X_{2}$ generation, and $X_{2}$ removal, respectively. The ODEs of Eqs. (2) provide a deterministic description of the chemical kinetics, with each $n_{i}$ corresponding to $N_{i}/\Omega$ in the continuum limit.

To obtain a stochastic description of the chemical kinetics Gillespie (1977); Van Kampen (1992), let us enumerate the possible state transitions when the system is in state $\mathbf{N}=\left(N_{1},N_{2}\right)$. From Eqs. (1), the possible state transitions and their stochastic transition rates (given above the arrows) are as follows:

[TABLE]

Here, we use $h_{r}\left(\mathbf{N}\right)$ to denote the number of ways reaction $r$ can take place and $\bar{\kappa}_{r}$ to denote the transition rate for reaction $r$ when this reaction can take place in one and only one way333The limit of $h_{r}\bar{\kappa}_{r}\Delta t$, evaluated at time $t$, as $\Delta t\rightarrow 0$ equals the probability that a reaction $r$ takes place somewhere in the system within a time window of $t$ to $t+\Delta t$. A good approximation of the probability is obtained with finite $\Delta t$ so long as $\Delta t$ is small enough such that the probability of two or more reactions of any kind occurring within the time window $[t,t+\Delta t)$ is much smaller than the probability of just one reaction.. For reaction $r=1$ (Eq. (1a)), there is only one way to inject a copy of $X_{1}$ into the system. Thus, $h_{1}=1$. For reaction $r=2$ (Eq. (1b)), the number of ways that two copies of $X_{1}$ can react to form $X_{2}$ is $N_{1}$ choose 2. Thus, $h_{2}=C^{N_{1}}_{2}=N_{1}(N_{1}-1)/2$. Finally, for reaction $r=3$ (Eq. (1c)), there are $N_{2}$ choose 1 ways to remove a copy of $X_{2}$ from the system. Thus, $h_{3}=C^{N_{2}}_{1}=N_{2}$.

The transition rates in Eqs. (3) are related to the mass-action rate laws of Eq. (2). The relationships are revealed by considering the continuum limit, where Eqs. (3) must be consistent with Eqs. (2). Eq. (3a) indicates that injection of $X_{1}$ increases the concentration of $X_{1}$ (i.e., $N_{1}/\Omega$), on average, by $\bar{\kappa}_{1}/\Omega$ per unit time, which is the rate characterized by the rate law $\kappa_{1}$ in Eq. (2a). Thus, we find that $\bar{\kappa}_{1}=\kappa_{1}\Omega$. Eq. (3b) indicates that conversion of $X_{1}$ to $X_{2}$ decreases the concentration of $X_{1}$ (i.e., $N_{1}/\Omega$), on average, by $2\times(1/2)\bar{\kappa}_{2}N_{1}(N_{1}-1)/\Omega$ per unit time (noting that each firing of the reaction in Eq. (3b) consumes two copies of $X_{1}$), which is the rate characterized by the rate law $\kappa_{2}n_{1}^{2}$ in Eq. (2a). Thus, taking $N_{1}\ll N_{1}^{2}$, which is appropriate in the continuum limit, we find that $\bar{\kappa}_{2}=\kappa_{2}/\Omega$. Eq. (3c) indicates that removal of $X_{2}$ decreases the concentration of $X_{2}$ (i.e., $N_{2}/\Omega$), on average, by $\bar{\kappa}_{3}N_{2}/\Omega$ per unit time, which is the rate characterized by the rate law $\kappa_{3}n_{2}$ in Eq. (2b). Thus, we find that $\bar{\kappa}_{3}=\kappa_{3}$. In generalVan Kampen (1992), for a $k^{\text{th}}$ order reaction $r$, $\bar{\kappa}_{r}$ is equal to $\kappa_{r}\Omega^{1-k}$.

III.3 Efficiency, accuracy and precision of standard scaling

The utility of scaling derives from the fact that transition rates (i.e., propensity functions) are extensive, meaning that they scale with system size. To verify this claim, consider the continuum limit, which is approached as $\Omega\rightarrow\infty$ and $N_{i}\rightarrow\infty$ for all $i$, but with each $N_{i}/\Omega$ held constant. Thus, in the finite but large population limit, $N_{i}\propto\Omega$ $\forall i$. Accordingly, the stochastic transition rate of any $k^{\text{th}}$ order reaction $r$ has magnitude given by the expression $\bar{\kappa}_{r}\mathcal{O}(\Omega^{k})$. However, from the relationship discussed above, $\bar{\kappa}_{r}\propto\kappa_{r}\Omega^{1-k}$, we can see that the expression simplifies to $\mathcal{O}(\Omega)$, indicating that the number of reaction events per unit time is proportional to the system size. By reducing system size $\Omega\rightarrow\lambda\Omega$, a speed up by a factor of $1/\lambda$ can be obtained as events occur $1/\lambda$ times less frequently.

A key idea of scaling is to take the stochastic dynamics within in a sub-volume of a system, or scaled system, to be representative of the dynamics in the whole system, or unscaled system. What this means in practice is as follows. In simulation of a scaled system, when a reaction event causes a change in the population of species $X_{i}$ by $\Delta N_{i}$, we take this change to correspond to a change of $\Delta N_{i}/\lambda$ in the unscaled state space. Thus, the method of static and homogeneous scaling can be seen as scaling down the values of all stochastic transition rates by a factor $\lambda$ and scaling up the values of all stoichiometric coefficients by a factor $1/\lambda$. As we will see later, this interpretation can be generalized to make the scaling of transition rates and stoichiometric coefficients adaptive and heterogeneous. However, let us first establish the accuracy and precision of static and homogeneous scaling, which we will pursue through an analysis of the reaction network of Eqs. (1).

We start by considering the chemical master equationGardiner et al. (1985); Van Kampen (1992) (CME) for the reaction network of Eqs. (1). For arbitrary $N_{1}$ and $N_{2}$, the CME includes an ODE of the following form:

[TABLE]

Recall that $\mathbb{P}_{\mathbf{N}}\left(t\right)$, written here as $\mathbb{P}_{N_{1},N_{2}}$, denotes the probability of system state $\mathbf{N}=(N_{1},N_{2})$ at time $t$.

If the system has large populations (i.e., the probability mass concentrates in regions where both $N_{1}$ and $N_{2}$ are large), we can perform the Kramers–Moyal expansionKramers (1940); Moyal (1949); Van Kampen (1992) to transform the CME into a partial differential equation describing a diffusive process. This procedure begins by replacing extensive variables (discrete populations) by intensive ones (concentrations/population densities). We introduce the intensive variables $n_{1}\equiv N_{1}/\Omega$ and $n_{2}\equiv N_{2}/\Omega$. We also introduce $\text{d}n\equiv\text{d}n_{1}\equiv\text{d}n_{2}\equiv 1/\Omega$, which is the size of the grid in the state space of $n_{1}$ and $n_{2}$. Finally, we replace discrete probabilities with probability densities such that $\mathbb{P}_{N_{1},N_{2}}(t)=\rho(n_{1},n_{2},t)\text{d}n_{1}\text{d}n_{2}=\rho(n_{1},n_{2},t)/\Omega^{2}$. Thus, after changing variables, we obtain

[TABLE]

At this stage, Eq. (5) is simply a rewritten form of Eq. (4), obtained after changing from extensive to intensive variables. As such, $\rho(n_{1},n_{2},t)$ can only be evaluated on a lattice $(n_{1},n_{2})\in\left\{0,1/\Omega,2/\Omega,\ldots\right\}^{2}$.

The first key idea of the Kramers–Moyal expansion for the CME is to extend the domain of the function $\rho$ from $(n_{1},n_{2})\in\left\{0,1/\Omega,2/\Omega,\ldots\right\}^{2}$ to $\mathbb{R}_{\geq 0}^{2}$ while assuming (5) holds true for any non-negative and real-valued $n_{1}$ and $n_{2}$. This transforms Eq. (5) into a difference equation for $\rho$. The second key idea is that under the assumption that $\Omega\gg 1$, the difference equation can be approximated by local derivatives.

Under the assumption that $\Omega\gg 1$ and $\text{d}n\ll 1$ in Eq. (5), which is viewed now as a difference equation, the Kramers–Moyal expansion procedure continues by Taylor expanding the terms $\rho\left(n_{1}\pm\text{d}n,n_{2}\pm\text{d}n\right)$ in Eq. (5). After all terms in each Taylor expansion higher than second order are dropped444Terms higher than second order are dropped because of the Pawula theoremRisken (1996), which states that any higher-order truncation ($\mathcal{O}\left(1/\Omega^{2}\right)$) fails to preserve the positivity of $\rho$., we obtain the following Fokker–Planck equation:

[TABLE]

The solution to Eq. (6) approximates the stochastic dynamics of the system of interest (Eq. (1)) when $\Omega$ is large but finite (i.e., $\Omega\gg 1$).

It is worth pointing out that in the continuum limit $\Omega\rightarrow\infty$, the Fokker–Planck equation (6) reduces to a Liouville equation,

[TABLE]

It can be shown that when the initial condition of the probability density function is a Dirac $\delta$-distribution, $\rho\left(x,y,t=0\right)=\delta(x-x_{0})\delta\left(y-y_{0}\right)$, the temporal evolution of the probability density $\rho\left(x,y,t\right)$ remains as a $\delta$-distribution and the peak of the distribution corresponds to the solution of Eqs. (2) for any given time $t>0$. Thus, Eq. (7) is a slightly generalized version of the equations for mass-action kinetics (Eqs. (2)). Eq. (7) captures the same information as Eqs. (2), but Eq. (7) permits a probabilistic distribution as an initial condition.

When $\Omega$ is large but finite, as is well-known, the dynamics described by a Fokker-Planck equation such as Eq. (6) are largely determined by the equation’s advection or drift terms, i.e., the terms in Eq. (6) that are preserved in the Liouville equation (Eq. (7)). However, additional terms, i.e., the terms with the operators $\partial_{n_{1}}^{2}$, $\partial_{n_{2}}^{2}$ and $\partial_{n_{1}}\partial_{n_{2}}$, are present in Eq. (6) and these terms introduce diffusion in the state space $\left(n_{1},n_{2}\right)$. It is these terms that characterize density fluctuations that arise from stochastic and discrete reaction events. Because the scale of diffusion is of order $\mathcal{O}\left(1/\Omega\right)$ (i.e., the variance), sample paths $\left(n_{1}\left(t\right),n_{2}\left(t\right)\right)$ exhibit fluctuations of order $\mathcal{O}\left(1/\Omega^{1/2}\right)$ (i.e., the standard deviation).

Now let us consider the scaled system, i.e., let us replace $\Omega$ with $\lambda\Omega$ in the above analysis. If $\lambda\Omega$ remains much larger than 1, and we perform the Kramers–Moyal expansion of the CME for the scaled system, we obtain

[TABLE]

This equation is nearly the same as Eq. (6). The only difference is the factor multiplying the diffusive terms. The factor is $1/(2\Omega)$ in Eq. (6) and $1/(2\Omega\lambda)$ in Eq. (8). This result generalizes to any chemical reaction system with large species populations (i.e., $N_{i}\gg 1$ for all $i$).

The significance of our comparison of Eqs. (6), (7), and (8) is threefold. First, it is known that first moments are dominated by drift (i.e., the terms with the operators $\partial_{n_{1}}$ and $\partial_{n_{2}}$) when $\Omega\gg 1$, and drift is identical in Eqs. (6), (7), and (8). Thus, we can expect simulations of scaled and unscaled systems to yield comparable estimates of first moments when $\Omega$ and $\lambda\Omega$ are large enough such that effects of multiplicative noise can be ignored. Second, we can expect scaling to yield overestimates of second moments, because diffusivity (i.e., variance) is enhanced by a factor of $1/\lambda$ in the scaled system relative to that in the unscaled system. Third, given that standard error of the mean (SEM) is calculated as $\sigma/\sqrt{\mathcal{N}}$, where $\sigma$ is the sample standard deviation and $\mathcal{N}$ is the number of samples, and given that we can expect the sample variance $\sigma^{2}$ to be amplified by a factor of $1/\lambda$ in estimates of first moments based on scaled simulations, the SEM for each first moment estimate is amplified by a factor of $1/\sqrt{\lambda}$ when using the same number of sample paths from scaled simulations as for estimates based on unscaled simulations. This loss of precision can be overcome by increasing the number of sample paths used. However, the number of sample paths must be increased by a factor of $1/\lambda$ to match the precision of estimates from unscaled simulations. This factor exactly matches the efficiency gain of a single scaled simulation. Thus, scaling yields acceleration only when a loss of precision in first moment estimates is permitted.

III.4 Limitation by small but critical populations

The Kramers–Moyal expansion of the scaled system is only valid when $\mathcal{O}\left(\lambda N_{i}\right)\gg 1$ for each $i$ and $\text{d}n\ll 1$. In the context of scaling, these conditions are satisfied only when $\lambda\gg 1/\min_{1\leq i\leq M}\left\{N_{i}\right\}$. When $\lambda$ is too small, one or more populations become discrete. Discrete populations have distinct dynamics, which are not captured by a Fokker–Planck equation Lin and Doering (2016); Lin and Galla (2016); Lin and Buchler (2018). Thus, the analysis of the previous section eventually becomes invalid and errors are introduced as $\lambda$ decreases and becomes too small. For example, for the system of Eq. (1), $\lambda=3\times 10^{-4}$ is evidently too small, as indicated by the errors in the estimates of the first moments that can be seen in Fig. 1(e) and (g). These errors arise because the scaled system is no longer representative of the unscaled system.

Scaling a system so that the smallest population vanishes (i.e., becomes 0) or reduces to order 1 might be acceptable in practice but requires careful consideration. On one hand, the population of concern may be of little significance for system behavior, and the erroneous behavior predicted for it can be safely ignored. On the other hand, if the population is critical for system behavior, then scaling will produce seriously misleading results.

Specialized methods are available for coping with discrete populations that influence system dynamics so long as these populations are always discrete and well-separated from other populations Bokes et al. (2013); Lin and Doering (2016); Lin and Galla (2016); Lin and Buchler (2018). However, for many systems, populations are not only distributed over a spectrum of scales but also dynamically changing scale, as illustrated in Fig. 2. We suspect that this behavior is generic, at least for biochemical systems, which limits the applicability of static and homogeneous scaling. Below, we present a new scaling approach that overcomes this limitation.

IV Adaptive and heterogeneous scaling

There are two barriers to applying static and homogeneous scaling. First, population scales may be scattered across a broad spectrum (Fig. 2). Second, populations may evolve dynamically across multiple scales (Fig. 2). To address these problems, we propose a new scaling approach, which we term partial scaling, that, in contrast with standard scaling, is adaptive and heterogeneous, as we will see. Both of these features of partial scaling arise from tying scaling to population sizes, on the fly. Scaling is heterogeneous because population sizes are heterogeneous, and scaling is adaptive because populations are dynamic. A particular, intentionally simple implementation of partial scaling is outlined as pseudocode in Algorithm 1. We call this algorithm the partial scaling algorithm (PSA).

With partial scaling, there is no longer a global scaling factor $\lambda$, as in standard scaling. Instead, at any given time $t$, there is a scaling factor $\lambda_{r}(t)$ for each reaction $r=1,\ldots,R$. These scaling factors are not set in advance of a simulation, nor statically. Rather, they are assigned values according to a dynamic update schedule (e.g., a schedule synchronized with time updates). This update schedule ensures that

[TABLE]

where $N_{c}$ is a threshold or critical population—this quantity is the one parameter of the method—and $N_{\min}^{r}\left(t\right)$ is the smallest population among those of the reactants and products of reaction $r$ at time $t$. It should be noted that the brackets around the ratio $N^{r}_{\min}/N_{c}$ in Eq. (9) denote the floor function and, furthermore, that there is no scaling (i.e., $\lambda_{r}=1$) for reaction $r$ if any reactant or product population is smaller than $N_{c}$.

The scaling factors $\{\lambda_{r}(t)\}_{r=1}^{R}$ defined by Eq. (9) are used to modify an SSA, such as Gillespie’s direct method Gillespie (1977), as follows. In the calculation of stochastic transition rates (e.g., as for Eqs. (3) as we have discussed), the rate for each reaction $r$ is calculated as usual but then scaled (down) by a factor $\lambda_{r}$. When reaction $r$ with rate scaled by $\lambda_{r}<1$ is selected to fire, to compensate for its reduced rate of firing, its stoichiometric coefficients are each temporarily scaled (up) by the factor $1/\lambda_{r}$. The amplified stoichiometric coefficients are then used to update the system state $\mathbf{N}$. In other words, $\mathbf{N}$ is modified in accordance with scaled population changes triggered by the selected reaction event. It should be noted that $1/\lambda_{r}\in\mathbb{N}$ by construction (Eq. (9)). Thus, scaled stoichiometric coefficients, like unscaled coefficients, are natural numbers. A specific example of how the scaling factors $\{\lambda_{r}(t)\}_{r=1}^{R}$ are used to modify an SSA is provided in Algorithm 1. It should be noted that Algorithm 1 reduces to Gillespie’s direct method if we statically assign each $\lambda_{r}$ a value of 1. Furthermore, standard scaling is achieved by replacing each reaction-specific adaptive scaling factor $\lambda_{r}$ in Algorithm 1 with a universal static scaling factor $\lambda$.

With partial scaling, population sizes are not directly scaled but the populations of reactants and products participating in any reaction $r$ for which $\lambda_{r}<1$ are effectively $\mathcal{O}(N_{c})$ or larger. Thus, a choice of $N_{c}\gg 1$ guarantees that stochastic dynamics do not become inappropriately bursty. Furthermore, any critical species $X_{i}$ with a small population $N_{i}<N_{c}$ will never be scaled out of the system. In other words, $N_{i}$ is never effectively scaled such that the population vanishes because $\lambda_{r}=1$ for any reaction having a reactant or product population less than $N_{c}$ (Eq. (9)).

IV.1 Accuracy and precision of partial scaling

As we will see in this section, provided that $N_{c}\gg 1$, partial scaling preserves means but not variances, similar to standard scaling.

Let us consider a time interval $\Delta t\ll(\sum_{r=1}^{R}h_{r}\bar{\kappa}_{r})^{-1}$, i.e., an interval that is much shorter than the expected waiting time to the next reaction event. Furthermore, consider an arbitrary species $X_{i}$, $i\in\{1,\ldots,M\}$, which participates, as a reactant or product, in a set of reactions with indices $I\subseteq\left\{1,\ldots,R\right\}$. We will use $I^{\prime}\subseteq I$ to denote the subset of these indices corresponding to reactions associated with a scaling factor less than 1. Thus, for each $r\in I^{\prime}$, $\lambda_{r}<1$, and for each $r\in I\backslash I^{\prime}$, $\lambda_{r}=1$. For an unscaled simulation, we will use $\delta N_{i,r}$ to denote the change of population $N_{i}$ brought about by a reaction $r\in I$ occurring within the time interval $\Delta t$. Given that multiple reaction events can occur within the time window $\Delta t$ (albeit with low probability, by construction), the total change in $N_{i}$ over this time window, which we will denote as $\delta N_{i}$, is given by

[TABLE]

Similarly, for a (partially) scaled simulation, the total population change is given by

[TABLE]

This expression is the same as that given above except for the introduction of scaling factors for stoichiometric coefficients.

The population changes considered above (i.e., the $\delta N_{i}$ and $\delta N_{i,r}$ terms in Eqs. 10 and 11) are random variables. Because $\Delta t$ is small (relative to the expected waiting time to the next reaction), we can take each of these random variables to be drawn from a Poisson distribution. Thus,

[TABLE]

where $\delta N_{i,r}^{\star}\equiv\delta N_{i,r}/\lambda_{r}$ denotes a population change triggered by firing of a reaction $r$ with a scaled rate $\lambda_{r}h_{r}\bar{\kappa}_{r}$.

For constants $c$, $\mu>0$, recall that $\mathbb{E}\left[c\cdot\text{Poisson}\left(\mu\right)\right]=c\mu$ and $\text{var}\left[c\cdot\text{Poisson}\left(\mu\right)\right]=c^{2}\mu$. Thus, for chemical species $X_{i}$ and reaction $r\in I^{\prime}$,

[TABLE]

It should now be clear that partial scaling yields unbiased estimates of first moments and overestimates of second moments (in a time interval $\Delta t\ll 1$), as is also the case for standard scaling. However, the bias in estimates of second moments, overall, is less with partial scaling than with standard scaling, because the processes in a (partially) scaled simulation that are unaffected by scaling have the same statistics as the processes in an unscaled simulation.

In partial scaling, each time a reaction $r$ (with scaled rate $\lambda_{r}h_{r}\bar{\kappa}_{r}$) is fired, a deterministic multiplier, $1/\lambda_{r}$, is applied to the reaction’s stochiometric coefficients for the purpose of appropriately updating the system state (Algorithm 1, line 24). As indicated by the analysis above, this multiplier is essential for obtaining unbiased estimates of first moments. Estimates of first moments would still be unbiased if the multiplier were instead chosen randomly with an expected value of $1/\lambda_{r}$; however, estimates of variance would be more biased for a random multiplier than for a deterministic multiplier666Let us use $XY$ to denote the population change of a reaction event (in a partially scaled simulation), where $X$ is a possibly random multiplier and $Y\sim\text{Poisson}\left(\lambda_{r}\bar{\kappa}_{r}h_{r}\Delta t\right)$. If the multiplier is distributed such that $\mathbb{E}\left[X\right]=1/\lambda_{r}$ and $X$ and $Y$ are independent, $\mathbb{E}\left[XY\right]=h_{r}\bar{\kappa}_{r}\Delta t$ (as desired for preservation of first moments) and $\text{var}\left[XY\right]=\text{var}\left[X\right]\lambda_{r}h_{r}\bar{\kappa}_{r}\Delta t\left(1+\lambda_{r}h_{r}\bar{\kappa}_{r}\Delta t\right)+h_{r}\bar{\kappa}_{r}\Delta t/\lambda_{r}$. Because all terms in this expression are positive, it follows that $\text{var}\left[XY\right]$ is least when $\text{var}\left[X\right]=0$, i.e., when $X$ is a deterministic variable. Recall that for two independently distributed random variables $X$ and $Y$, $\mathbb{E}\left[XY\right]=\mathbb{E}\left[X\right]\mathbb{E}\left[Y\right]$ and $\text{var}\left[XY\right]=\text{var}\left[X\right]\,\text{var}\left[Y\right]+\text{var}\left[X\right]\,\mathbb{E}^{2}\left[Y\right]+\mathbb{E}^{2}\left[X\right]\,\text{var}\left[Y\right]$..

One may wonder if partial scaling yields unbiased first-moment estimates in the short term but not in the long run. We expect that the answer is model dependent. If a model permits noise-induced phenomena involving chemical species participating in reactions with scaled rates, then estimates may become biased because of demographic fluctuations. In other words, systematic errors may arise with partial scaling from behavior that depends on second or higher moments. Below, in part to address this issue, we turn our attention to numerical experiments, including simulations for a model that permits bistable stochastic switching, the TCR model Lipniacki et al. (2008). As we will see, partial scaling yields unbiased estimates of first moments (and also higher-order moments), even for this model.

V Comparison of scaling methods

To facilitate benchmarking of partial scaling, we implemented a variation of Gillespie’s direct method, an additional variation that incorporates standard scaling, and a variation of PSA (Algorithm 1). Importantly, we implemented these methods so as to eliminate all unnecessary implementation differences. We used our code to perform stochastic simulations based on three published models Kochańczyk et al. (2017); Rubenstein et al. (2007); Lipniacki et al. (2008) with the goal of evaluating the relative efficiency of each of the three methods (no scaling, standard scaling and partial scaling) and the relative accuracy of standard and partial scaling. In simulations, we used a model-specific pre-generated dependency graph to aid in updates of propensity functions but we did not use sorting of the propensity functions to accelerate simulations as in the Gibson-Bruck method Gibson and Bruck (2000). In our simulations, we considered different settings for the algorithmic parameters $\lambda$ (the static scaling factor in standard scaling) and $N_{c}$ (the critical population size used in partial scaling).

The models that we considered, all for biological systems, were chosen because each was originally analyzed using exact stochastic simulation and each challenges the application of standard scaling. The first of the three models characterizes activation of ERK Kochańczyk et al. (2017). The dynamics of this model are oscillatory. Populations span multiple orders of magnitude and their scales change dynamically (Fig. 2(a)). The second of the three models characterizes prion protein aggregation Rubenstein et al. (2007), which involves polymerization-like reactions and seeded nucleation of aggregates. As before, the populations considered in this model evolve over multiple scales, and the initial abundance of prion protein is discrete (Fig. 2(b)). The third of the three models characterizes TCR signaling Lipniacki et al. (2008). The behavior characterized by this model includes stochastic bistable switching. In other words, for a particular regime of behavior, which we considered in our simulations (as is evident from inspection of individual sample paths), intrinsic noise is capable of driving the system from one stable fixed point to another (and back). This type of behavior has been intensely studied in the context of genetic regulatory circuits Walczak, Sasai, and Wolynes (2005); Roma et al. (2005); Warren and ten Wolde (2005); Assaf, Roberts, and Luthey-Schulten (2011); Strasser, Theis, and Marr (2012); Bokes et al. (2013); Lu et al. (2013); Lin and Galla (2016). For the TCR model, unlike for many models of gene regulation that have been considered, many of the populations are large (i.e., near the continuum limit) (Fig. 2(c)).

In Figs. 3 and 4, we plot predicted time courses for selected chemical species considered in the ERK modelKochańczyk et al. (2017) and the prion modelRubenstein et al. (2007), respectively. Each point in each time course corresponds to the mean population for the indicated chemical species estimated on the basis of 500 sample paths. (This number of sample paths was deemed sufficient for accurate estimation of statistical quantities of interest.) In each figure, the top row of panels shows results of calculations based on either exact simulation or standard scaling (for different values of $\lambda$), as indicated in the legend at right. Similarly, in each figure, the bottom row of panels shows results of calculations based on either exact simulation or partial scaling (for different values of $N_{c}$), as indicated in the legend at right. As can been, standard scaling yields estimates that deviate markedly from those based on exact simulation. In contrast, even with aggressive scaling ($N_{c}=300$), partial scaling yields first-moment estimates that are barely distinguishable from those based on exact simulation.

In Fig. 5, we plot time-dependent marginal probability distributions (estimated on the basis of $10^{4}$ sample paths) for two selected chemical species considered in the TCR modelLipniacki et al. (2008). (We considered $10^{4}$ vs. 500 sample paths because here we are estimating distributions vs. properties of distributions.) In this figure, we focus on the initial transient shown in Fig. 2(c). The dynamics of this model are inherently stochastic and poorly represented by first moments. As can be seen, standard scaling yields marginalized probability distributions that differ markedly from those obtained from exact simulations (top row), whereas partial scaling yields approximate results that are very close to the exact results (bottom row).

Additional accuracy results, based on error measures introduced in the Appendix, and timing results for the three benchmark problems are summarized in Fig. 6 and Table 1. At least for these problems, partial scaling allows for more significant speed ups than standard scaling with the introduction of much less error in the estimates of first and second moments.

In Fig. 7, for one particular setting for $N_{c}$, we have plotted $\lambda_{r}$ as a function of time $t$ for every reaction included in each of the three models considered above. These plots show that partial scaling is highly dynamic and multiscale.

VI Implementation

We added a partial scaling feature to the SSA implemented in BioNetGenHarris et al. (2016), an open-source, general-purpose simulation package used by biological modelers. The SSA implemented in BioNetGen is an efficient variation of Gillespie’s direct method Gillespie (2007), which incorporates various ideas for optimizing simulation efficiency, such as on-the-fly generation (vs. pre-generation) of the list of reactions in which (populated) chemical species can participate Lok and Brent (2005); Faeder et al. (2005). BioNetGen is designed for compatibility with models defined using BNGL Faeder, Blinov, and Hlavacek (2009), a language for specifying deterministic and stochastic models for well-mixed (bio)chemical reaction kinetics and for specifying simulations based on such models. BioNetGen also supports models defined using SBMLHucka et al. (2018), such as those available in BioModels DatabaseChelliah, Laibe, and Le Novère (2013). To invoke partial scaling using BioNetGen’s simulate command, the method argument should be set to ssa or psa and a new poplevel argument, which is used only with partial scaling, should be assigned a non-negative integer value. For accuracy, care should be taken to assign a value much larger than 1. The setting for poplevel determines the value of $N_{c}$ and thereby the aggressiveness of scaling. Examples of usage are provided in BNGL files that we have made available online Lin (2019a, b); Lin, Y T (2019). These files define the ERK activationKochańczyk et al. (2017), prion protein aggregationRubenstein et al. (2007), and TCR signalingLipniacki et al. (2008) models considered in Figs. 3–5. It should be noted that our implementation of partial scaling in BioNetGen (version 2.5 or higher) is a so-called generate-first method, meaning that it requires an enumeration of chemical species and the individual reactions in which these species are able to participate. In contrast, so-called network-free methods do not require an explicit enumeration of chemical species or reactions Suderman et al. (2018).

VII Discussion and conclusion

The procedure that we have called standard scaling is commonly used for accelerating stochastic simulations of systems with large population sizes. However, this method has limited applicability and it may be challenging to provide an appropriate/optimal setting for the method’s one parameter, a universal static scaling factor $\lambda$. A simulation speed up can only be attained by sacrificing the ability to calculate unbiased estimates of moments higher than the first, and even first-moment information can be calculated incorrectly with overly aggressive scaling. Furthermore, the speed up that is attainable is constrained by the smallest population size measured over the entire time window of simulation. Thus, in practice, careful application of the method may require trial-and-error numerical experiments to find a suitable best setting for the scaling factor $\lambda$. For this reason, the method is especially problematic when one wishes to use it as part of a parameter identification procedure, or any procedure involving parameter variation, because the best $\lambda$ setting depends in a non-obvious way on parameter values (initial conditions and rate constants).

As we have seen for three non-trivial benchmark problems (Figs. 3–6, Table 1), partial scaling significantly outperforms standard scaling in two important ways. First, partial scaling yields greater acceleration for the same or better accuracy (Fig. 6 and Table 1). Second, as illustrated in Fig. 5, partial scaling better preserves second-moment information (Table 1). Another attractive feature of partial scaling is the better ability of the user to avoid overly aggressive scaling without any requirement for numerical experiments. With partial scaling, overly aggressive scaling is only possible if the setting for $N_{c}$ is near or below 1, which is not recommended. In contrast, with standard scaling, any setting for $\lambda$ must be tested in numerical experiments. As noted earlier, a trial-and-error procedure may be required to find a suitable setting for $\lambda$ and the suitability of a setting may change with a change of parameter values. Thus, standard scaling seems especially disadvantaged in comparison to partial scaling for simulations within the context of a parameter identification procedure (i.e., a fitting procedure).

Partial scaling is similar to the probability-weighted dynamic Monte Carlo (PW-DMC) methodResat, Wiley, and Dixon (2001). In this method, reaction rates are scaled dynamically according to user-specified rules. The rules define scaling factors for rates that fall into specified ranges. Relative to partial scaling, PW-DMC has more parameters: a set of ranges and the associated scaling factors vs. partial scaling’s single parameter, a critical population size ($N_{c}$). Unlike $N_{c}$, the PW-DMC parameters must be set through a trial-and-error process. Because rates are scaled without regard to reactant population sizes, PW-DMC can introduce errors where partial scaling does not. These errors arise, for example, when a bimolecular reaction has reactants with disparate population sizes, such that one is discrete but the other is large enough to yield a rate that qualifies for scaling. Scaling in this scenario can cause the discrete population to become negative.

We have not directly compared partial scaling against hybrid methods or $\tau$-leaping methods. However, partial scaling has clear advantages over these methods. Partial scaling is far easier to implement than hybrid methods, and partial scaling is more broadly applicable than $\tau$-leaping methods. Partial scaling is useful when population sizes are distributed smoothly across multiple scales (i.e., without clear separation of discrete and continuous populations), because there is no a priori requirement for classification of population sizes as either discrete or continuous. Such a requirement can be highly problematic when the classification changes during the time window of simulation, as could be the case for a system that exhibits oscillatory behavior. The chief disadvantage of $\tau$-leaping methods is the lack of guarantee of the existence of a time window $\tau$ having the necessary properties. In our experience, hybrid and $\tau$-leaping methods are useful in restricted circumstances. We expect partial scaling to be useful for a wider array of circumstances.

Although partial scaling is an approximate SSA (vs. an exact SSA), such methods have important applications. For example, in parameter identification procedures, which typically entail numerous simulation runs, success is highly dependent on simulation efficiency and inexactness is a lesser concern, especially when parameter estimation is based on noisy data. The approximations of partial scaling are very likely to be tolerable, even when statistical distribution data are being used in a fitting procedure. As the optimization algorithm in a fitting procedure converges, numerous samples around a local optimum are inevitably generated. By averaging over a moving window of trial parameter sets, for example, one can, in principle, compensate for the noisy (but unbiased) estimates of first moments when determining goodness of fit.

An attractive future use of partial scaling is acceleration of network-free simulationSuderman et al. (2018), which is an important outstanding problem, as we have discussed elsewhereLin et al. (2018). One reason to pursue this goal is the high cost of fitting when network-free simulation is necessaryThomas et al. (2016), as is the case for the model of Chylek et al.Chylek et al. (2014b). Unfortunately, application of partial scaling requires system state to be tracked in terms of populations and no currently available, general-purpose implementation of a network-free simulation algorithm tracks system state in this way. Rather, state is followed in terms of the states of individual biomolecular sitesSuderman et al. (2018), which reflects an agent-based or particle-based approach to simulation. However, as discussed by Liu et al. Liu et al. (2010), it is feasible to develop a network-free simulation algorithm in which system state is tracked in terms of populations, just those of the chemical species with non-zero population sizes. (An exhaustive enumeration of potentially populated chemical species is typically impracticable whenever network-free simulation is under consideration.) Liu et al. Liu et al. (2010) presented a specialized algorithm that uses this approach to state tracking, which they termed the full-scale SSA method. However, this method is specialized for a subset of rule-based models—it doesn’t consider systems in which molecules interact to form assemblies (complexes). Thus, an interesting future research direction would be developing a more general version of the full-scale SSA method that is compatible with partial scaling.

Acknowledgments

We thank Danny Perez and Arthur F. Voter for helpful discussions. We also thank the research group of James R. Faeder, which maintains the BioNetGen code base with support from the National Institute of General Medical Sciences (NIGMS) of the National Institutes of Health (NIH) (grant no. P41GM103712). W.S.H acknowledges support from NIGMS/NIH (grant no. R01GM111510). Y.T.L. and S.F. acknowledge support from the Center for Nonlinear Studies (CNLS). CNLS is funded by the Laboratory-Directed Research and Development program at Los Alamos National Laboratory, which is operated by Triad National Security, LLC for the National Nuclear Security Administration of the U.S. Department of Energy (contract no. 89233218CNA000001).

Appendix

First-moment errors of scaled simulations per species and overall are reported in Fig. 6 and Table 1, respectively. These errors were calculated as described here. Below, we also describe the means used to compare the aggressiveness of standard and partial scaling fairly (via calculation of an effective scaling factor, $\lambda_{\rm eff}$) and the means used to characterize second moments obtained from either an unscaled or scaled simulation in terms of a multivariate coefficient of variation.

For simulations based on any of the three models of interest involving scaling (either standard scaling or partial scaling), the mean population of chemical species $X_{i}$ (for $i=1,\ldots,M$) at report time $t$ was calculated from $n_{s}$ sample paths. Similarly, on the basis of $n_{s}$ sample paths, we calculated the corresponding mean obtained via exact simulation. In this way, for either scaling scheme of interest (standard scaling with $0<\lambda<1$ or partial scaling with $N_{c}\gg 1$), we obtain a deviation vector

[TABLE]

where $\overline{\mathbf{N}^{(\lambda,N_{c})}}\left(t\right)$ denotes the time-dependent $M$-dimensional vector of sample path means calculated from scaled simulations and $\overline{\mathbf{N}^{(0)}}(t)$ denotes the time-dependent $M$-dimensional vector of sample path means calculated from exact simulations. In using the notation $\overline{\mathbf{N}^{(\lambda,N_{c})}}$, we adopt the convention that replacement of $(\lambda,N_{c})$ with $(1,\infty)$, $(\lambda,\infty)$ and $(1,N_{c})$ indicates no scaling, standard scaling, and partial scaling, respectively.

To quantify the accuracy of first moments obtained from scaled simulations, we define the following error measure:

[TABLE]

Here, $\bar{N}_{i}\left(t\right)$ denotes the averaged population of species $i$ at time $t$ determined from a set of unscaled simulations. Similarly, $\bar{N}_{i}^{\left(\lambda,\bar{N}_{c}\right)}\left(t\right)$ denotes the corresponding quantity determined from an equal number of scaled simulations. We note that $\varepsilon_{1}^{i}$ is continuous with respect to $\bar{N}_{i}^{\left(\lambda,N_{c}\right)}$ for arbitrary $\bar{N}_{i}\geq 0$. The derivative of $\varepsilon_{1}^{i}$ with respect to $\bar{N}_{i}^{\left(\lambda,N_{c}\right)}$ is also continuous. With the definition given above, error does not diverge to infinity as $\bar{N}_{i}$ approaches 0.

For each numerical experiment of interest (aimed at quantifying the errors introduced by scaling), we calculated $\varepsilon_{1}^{i}$ for each $i$ at a set of evenly spaced report times, from 0 to $T$. We then calculated the temporal average $\left\langle\varepsilon_{1}^{i}\right\rangle$, which is given by

[TABLE]

where $\Delta t$ is the spacing between report times. Quantities found using this equation are plotted in Fig. 6.

The aggressiveness of partial scaling is related to $N_{c}$, and the aggressiveness of standard scaling is related to $\lambda$. These quantities are not directly comparable. To fairly compare the aggressiveness of partial and standard scaling, we define $\lambda_{\text{eff}}$ as the normalized value of $\kappa$ in a scaled simulation (Algorithm 1, line 14) averaged over all discrete report times. We take the normalization constant to be the temporally averaged value of $\kappa$ in a corresponding unscaled simulation. Thus, in the case of standard scaling, $\lambda_{\text{eff}}$ is exactly $\lambda$. In Fig. 6, we indicate the values of $\lambda_{\text{eff}}$ calculated for partially scaled simulations with different values of $N_{c}$.

In Table 1, for each numerical experiment, we report a summary statistic for overall error, which is defined as

[TABLE]

To characterize second moments, we calculated the Van Valen multivariate coefficient of variationVan Valen (1974) of the $M$-dimensional vector $\mathbf{N}(t)$ describing the species popualtions at each report time:

[TABLE]

where $\left||\bar{\mathbf{N}}\left(t\right)|\right|_{2}^{2}$ is the square of the $L^{2}$ norm of the averaged populations and $\Sigma\left(t\right)$ is the covariance matrix. To report a second-moment summary statistic, we averaged over report times:

[TABLE]

Values of $\left\langle CV_{\text{VV}}\right\rangle$ are reported in Table 1 for both scaled and unscaled simulations. Second-moment errors attributable to scaling are reflected in the difference between the values of $\left\langle CV_{\text{VV}}\right\rangle$ given for a scaled simulation and the corresponding unscaled simulation.

Values given in Table 1 for $\left\langle\varepsilon_{1}\right\rangle$ and $\left\langle CV_{\text{VV}}\right\rangle$ are based on the following choices for $\Delta t$, $T$, and $n_{s}$, the number of sample paths. For the ERK activation modelKochańczyk et al. (2017), $\Delta t=8.64$ sec, $T=8,640$ sec, and $n_{s}=500$. For the prion protein aggregation modelRubenstein et al. (2007), $\Delta t=0.01$ day, $T=300$ days, and $n_{s}=500$. For the TCR signaling modelLipniacki et al. (2008), $\Delta t=10$ sec, $T=10^{4}$ sec, and $n_{s}=10^{4}$.