Constructing stochastic models from deterministic process equations by propensity adjustment

Abstract

Background

Gillespie's stochastic simulation algorithm (SSA) for chemical reactions admits three kinds of elementary processes, namely, mass action reactions of 0^th, 1^st or 2^nd order. All other types of reaction processes, for instance those containing non-integer kinetic orders or following other types of kinetic laws, are assumed to be convertible to one of the three elementary kinds, so that SSA can validly be applied. However, the conversion to elementary reactions is often difficult, if not impossible. Within deterministic contexts, a strategy of model reduction is often used. Such a reduction simplifies the actual system of reactions by merging or approximating intermediate steps and omitting reactants such as transient complexes. It would be valuable to adopt a similar reduction strategy to stochastic modelling. Indeed, efforts have been devoted to manipulating the chemical master equation (CME) in order to achieve a proper propensity function for a reduced stochastic system. However, manipulations of CME are almost always complicated, and successes have been limited to relative simple cases.

Results

We propose a rather general strategy for converting a deterministic process model into a corresponding stochastic model and characterize the mathematical connections between the two. The deterministic framework is assumed to be a generalized mass action system and the stochastic analogue is in the format of the chemical master equation. The analysis identifies situations: where a direct conversion is valid; where internal noise affecting the system needs to be taken into account; and where the propensity function must be mathematically adjusted. The conversion from deterministic to stochastic models is illustrated with several representative examples, including reversible reactions with feedback controls, Michaelis-Menten enzyme kinetics, a genetic regulatory motif, and stochastic focusing.

Conclusions

The construction of a stochastic model for a biochemical network requires the utilization of information associated with an equation-based model. The conversion strategy proposed here guides a model design process that ensures a valid transition between deterministic and stochastic models.

Background

Most stochastic models of biochemical reactions are based on the fundamental assumption that no more than one reaction can occur at the exact same time. A consequence of this assumption is that only elementary chemical reactions can be converted directly into stochastic analogues [1]. These include: 1) zero-order reactions, such as the generation of molecules at a constant rate; 2) first-order reactions, with examples including elemental chemical reactions as well as transport and decay processes; and 3) second-order reactions, which include heterogeneous and homogeneous bimolecular reactions (dimerization). Reactions with integer kinetic orders other than 0, 1 and 2 are to be treated as combinations of sequential elementary reactions. The advantage of the premise of non-simultaneous reaction steps is that the stochastic reaction rate can be calculated from a deterministic, equation-based model with some degree of rigor, even though the derivation is usually not based on first physical principles but instead depends on other assumptions and on macroscopic information, such as a fixed rate constant in the equation-based model. The severe disadvantage is that this rigorous treatment is not practical for modelling larger biochemical reaction systems. The reasons include the following. First, in many cases, elementary reaction rates are not available. Secondly, even in the case that all reaction parameters are available, the computational expense is very significant when the system involves many species and reactions, and this fact ultimately leads to a combinatorial explosion of required computations. Within a deterministic modelling framework, the common practice in this situation is to fit the transient and steady-state experimental data with a phenomenological, (differential) equation-based model, which explicitly or implicitly eliminates or merges some intermediate species and reactions. The best-known examples are probably Michaelis-Menten and Hill rate laws, which are ultimately explicit, but in truth approximate a multivariate system of underlying chemical processes.

Similar model reduction efforts have been carried out for stochastic modelling. For instance, the use of a complex-order function (which corresponds to a reduced equation-based model) was shown to be justified for some types of stochastic simulations. A prominent example is again the Michaelis-Menten rate law, which can be reduced from a system of elementary reactions to an explicit function by means of the quasi-steady-state assumption (see Result section and [2, 3]). However, model reduction within the stochastic framework has proven to be far more difficult than in the deterministic counterpart. The difficulties are mainly due to the fact that the reduction must be carried out on the chemical master equation (CME). This process is nontrivial and has succeeded only in simple cases.

In general, the construction of a stochastic model for a large biochemical network requires the use of information available from an equation-based model. In the past, several strategies have been proposed for this purpose and within the context of Gillespie's exact stochastic simulation algorithm (SSA; [1]) and its variants [4]. For example, Tian and Burrage [5] proposed that a stochastic model could be directly formulated from the deterministic model through a Poisson leaping procedure. However, a rigorous mathematical justification for such a conversion is lacking. Typical moment-based approaches [6–8] derive ODEs for the statistical moments of the stochastic model from an equation-based model where the 0^th, 1^st and 2^nd order reactions follow mass action rate laws. More recently the moment method was extended to cover models consisting of rational rate laws [9]. Moreover, it was realized that the moment method is complementary to, but cannot fully replace, stochastic simulations, because it does not cover situations like genetic switches [6, 10].

In this article, we explore the mathematical connection between deterministic and stochastic frameworks for the pertinent case of Generalized Mass Action (GMA) systems, which are frequently used in Biochemical Systems Theory (BST; [11–13]). Specifically, we address two questions: First, under what conditions can a deterministic, equation-based model be converted directly into a stochastic simulation model? And second, what is a proper way of implementing this conversion? We will develop a method to answer these questions and demonstrate it for functions in the canonical power-law format of GMA systems. However, the results are applicable to other functions and formats as well, as we will demonstrate with several examples.

Representations of systems of biochemical reactions

Consider a well-stirred biochemical reaction system with constant volume and temperature, where N _s different chemical species ${\left\{S_s\right\}}_{s=1}^{N_s}$, interact through N _r unidirectional reaction channels ${\left\{R_r\right\}}_{r=1}^{N_r}$. Each reaction channel can be characterized as

$$R_r:{\underset{-}{v}}_{r1}{S}_1+\dots +{\underset{-}{v}}_{r{N}_s}{S}_{N_s}\underset{}{\overset{k_r}{\to }}{\bar{v}}_{r1}{S}_1+\dots +{\bar{v}}_{r{N}_s}{S}_{N_s},$$

where ${\underset{-}{v}}_{rs}$ and ${\bar{v}}_{rs}$ are the counts of molecular species S _s consumed and produced due to reaction R _r , respectively, and k _r is the rate constant. The changed amount of $S_s{v}_{rs}={\bar{v}}_{rs}-{\underset{-}{v}}_{rs}$, which is due to the firing of reaction R _r , defines the stoichiometric coefficient of S _s in R _r . The stoichiometric coefficients of all species can be summarized according to each reaction R _r in the stoichiometric vector

$${\mathbf{v}}_r\triangleq \left\{\begin{array}{c}\hfill v_{r1}\hfill \\ \hfill ⋮\hfill \\ \hfill v_{r{N}_s}\hfill \end{array}\right\}\in {\mathbb{Z}}^{N_s}.$$

The stoichiometric vectors of all reactions can further be arranged as the stoichiometric matrix of the system

$$V\triangleq \left[{\mathbf{v}}_1,\dots ,{\mathbf{v}}_{N_r}\right]\in {\mathbb{Z}}^{N_s\times N_r}.$$

The size of the system is defined as Φ = AU, where A is the Avogadro number and U is the reaction volume.

The modelling of biochemical reaction networks typically uses one of two conceptual frameworks: deterministic or stochastic. In a deterministic framework, the state of the system is given by the a non-negative vector $\left[\mathbf{X}\left(t\right)\right]={\left[\left[X_1\left(t\right)\right],\dots ,\left[X_{N_s}\left(t\right)\right]\right]}^T\in {ℝ}^{N_s}$, where component [X _s (t)] represents the concentration of species S _s , measured in moles per unit volume. The temporal evolution of the state of the system is modelled by a set of ordinary differential equations, which in our case are assumed to follow a generalized mass action (GMA) kinetic law. By contrast, in a stochastic framework, the state of the systems is characterized by a vector $\mathbf{x}\left(t\right)={\left[x_1\left(t\right),\dots ,x_{N_s}\left(t\right)\right]}^T\in {\mathbb{Z}}^{N_s}$, whose values are non-negative integers. Specifically, x _s (t) = Φ [X _s (t)] is the count of S _s molecules, which is a sample value of the random variable X _s (t). The system dynamics of this process is typically described with the chemical master equation (CME). Both GMA and CME will be discussed in detail in the following sections.

Motivation for the power-law formalism: reactions in crowded media

Power-law functions with non-integer kinetics have proven very useful in biochemical systems analysis, and forty years of research have demonstrated their wide applicability (e.g., see [11–13]). Generically, this type of description of a biochemical reaction can be seen either as a Taylor approximation in logarithmic space or as a heuristic or phenomenological model that has been applied successfully hundreds of times and in different contexts, even though it is difficult or impossible in many situations to trace it back to first mechanistic principles. A particularly interesting line of support for the power-law format can be seen in the example of a bimolecular reaction occurring in a spatially restricted environment. Savageau demonstrated that the kinetics of such a reaction can be validly formulated as a generalization of the law of mass action, where non-integer kinetic orders are allowed [14, 15]. Neff and colleagues [16–18] showed with careful experiments that this formulation is actually more accurate than alternative approaches.

Within the conceptual framework of power-law representations, the rate of the association reaction between molecules of species S₁ and S₂ is given as $k{\left[X_1\left(t\right)\right]}^{f_1}{\left[X_2\left(t\right)\right]}^{f_2}$. Here, k is the rate constant and f₁ and f₂ are real-valued kinetic orders, which are no longer necessarily positive integers as it is assumed in a mass action law. As an example, consider the reversible bimolecular reaction $S_1+S_2\underset{k_b}{\overset{k_f}{\rightleftharpoons }}{S}_3$. Like Neff and colleagues [17], we begin by formulating a discrete update function for the population of S₃ molecules as

$$x_3\left(t+\mathrm{\Delta }t\right)\mathsf{\text{ - }}{x}_3\left(t\right)=f\left(\left[X_1\right],\left[X_2\right]\right)\mathrm{\Delta }t{x}_1{x}_2\mathsf{\text{ - }}g\left(\left[X_3\right]\right)\mathrm{\Delta }t{x}_3.$$

(1)

The first term on the right-hand side of this equation, f ([X₁], [X₂])Δt x₁x₂, describes the production of S₃: it depends on the totality of possible collisions x₁ x₂ and also on some fraction f ([X₁], [X₂])Δt that actually reacts and forms the product. In a dilute environment, f ([X₁], [X₂]) equals a traditional rate constant, and the reaction obeys the law of mass action, while in a spatially restricted environment, such as the cytoplasm, one needs to take crowding effects into account. As shown in Savageau [14, 15], the desired fraction of a reaction in a crowded environment becomes a rate function that depends on the current concentrations of S₁ and S₂. The second term, g ([X₃]) Δtx₃, describes the fraction g ([X₃]) Δt of species S₃ that dissociates back into S₁ and S₂. This fraction may depend on some functional form of [X₃] because in a crowded environment the complex may not be able to dissociate effectively. Thus, rate constants in the generalized mass action setting become rate functions (cf. [17]).

By taking the limit Δt → 0, one obtains the differential equation

$$\frac{d{x}_3}{dt}=f\left(\left[X_1\right],\left[X_2\right]\right)x_1{x}_2\mathsf{\text{ - }}g\left(\left[X_3\right]\right)x_3.$$

(2)

Savageau used Taylor series expansion to approximate the functions f and g in the logarithmic space (log [X₁], log[X₂]) around some operating point (a, b). The result for f is

$$\begin{array}{c}logf\left(\left[X_1\right],\left[X_2\right]\right)\triangleq F\left(log\left[X_1\right],log\left[X_2\right]\right)\\ ={\left.F\left(a,b\right)+\frac{\partial f\left(\left[X_1\right],\left[X_2\right]\right)}{\partial \left[X_1\right]}\right|}_{\left(a,b\right)}\left(log\left[X_1\right]-a\right)\\ +{\left.\frac{\partial f\left(\left[X_1\right],\left[X_2\right]\right)}{\partial \left[X_2\right]}\right|}_{\left(a,b\right)}\left(log\left[X_2\right]-b\right)+\mathsf{\text{HOT}}\\ \approx k_f+\alpha log\left[X_1\right]+\beta log\left[X_2\right],\end{array}$$

(3)

where k _f , α, and β are constants related to the chosen operating point (a, b). The final step is achieved by ignoring all higher order terms (HOT) beyond the constant and linear terms. Transformation back to the Cartesian space yields

$$f\left(\left[X_1\right],\left[X_2\right]\right)\approx k_a{\left[X_1\right]}^{\alpha }{\left[X_2\right]}^{\beta },k_a=e^{k_f}.$$

(4)

The same procedure leads to the power-law expression for the degradation term: g ([X₃]) ≈ k _d [X₃]^γ. By combining constants we arrive at a power-law representation for the dynamics of species S₃ as

$$\begin{array}{cc}\hfill \frac{d\left[X_3\right]}{dt}& =k_a{\left[X_1\right]}^{\alpha }{\left[X_2\right]}^{\beta }\left[X_1\right]\left[X_2\right]-k_d{\left[X_3\right]}^{\gamma }\left[X_3\right]\hfill \\ =k_a{\left[X_1\right]}^a{\left[X_2\right]}^b-k_d{\left[X_3\right]}^c,\hfill \end{array}$$

(5)

where a = α + 1, b = β + 1, and c = γ + 1. As long as k _f , k _d , a, b and c remain more or less constant throughout a relevant range, the power-law model is mathematically well justified. In actual applications, the values of rate constants and kinetic orders can be estimated from experimental data [19]. When the functions f and g are originally not in power-law format, they can be locally approximated by power-law functions with a procedure similar to the one shown above (Equations (3) to (5)). An illustration will be given in the example section.

The Generalized Mass Action (GMA) format

In the GMA format within Biochemical Systems Theory, each process is represented as a univariate or multivariate power-law function. GMA models may be developed de novo or as an approximation of some other nonlinear rate laws. GMA models characterize the time evolution of the system state given that the system was in the state X (t₀) at some initial time t₀. Generically, the state of the system is changed within a sufficiently small time interval by one out of the N _r possible reactions that can occur in the system. The reaction velocity through reaction channel R _r is:

$$\left[\frac{{\left[X_1\left(t\right)\right]}^{\prime }}{v_{r1}}=\dots =\frac{{\left[X_{N_s}\left(t\right)\right]}^{\prime }}{v_{r{N}_s}}\right]=k_r\prod _{s=1}^{N_s}{\left[X_s\left(t\right)\right]}^{f_{rs}}$$

(6)

for those $v_{rs}={\bar{v}}_{rs}-{\underset{-}{v}}_{rs}\ne 0$, s = 1, ..., N _s . As shown in the example of a bimolecular reaction, the kinetic order f _rs associated with species S _s captures the effects of both reactant properties (such as the stoichiometric coefficient v _rs ) and environmental influences (such as temperature, pressure, molecular crowding effects, etc.). Therefore f _rs does not necessary equal an integer v _rs , which is assumed to be the case in mass action kinetics, but is possibly real-valued and may be negative. Summing up the contributions of all reactions, one obtains a GMA model describing the dynamics of S _s as

$$\frac{d}{dt}\left[X_s\left(t\right)\right]=\sum _{r=1}^{N_r}{v}_{rs}{k}_r\prod _{s=1}^{N_s}{\left[X_s\left(t\right)\right]}^{f_{rs}}$$

(7)

for every s = 1, ..., N _s . Each reaction contributes either a production flux or a degradation flux to the dynamics of a certain species. Positive terms (v _rs > 0) represent the production of S _s , while negative terms (v _rs < 0) describe degradation. If f _rs is positive, then S _s accelerates the reaction R _r ; a negative value represents that S _s inhibits the reaction, and f _rs = 0 implies that S _s has no influence on the reaction. The rate constant k _r for reaction R _r , is either positive or zero. Both, the rate constant and the kinetic order, are to be estimated from data.

Proper use of equation-based functions for stochastic simulations

The fundamental concept of a stochastic simulation is the propensity function α(X), and α(X)dt describes the probability that a reaction will change the value of a system variable within the next (infinitesimal) time interval (t, t +dt). While a formal definition will be given later (Equation 18), it is easy to intuit that the propensity function is in some sense analogous to the rate in the corresponding deterministic model. In fact, the propensity function is traditionally assumed to be α(X) = f _s (X), if the deterministic model is X _s ' = f _s (X, t), s = 1, ..., N _s . However, a proper justification for this common practice is by and large missing. Indeed, we will show that the direct use of a rate function as the propensity function in a stochastic simulation algorithm requires that at least one of the following assumptions be true:

• 1) f is a linear function;

• 2) the reaction is monomolecular;

• 3) all X _i in the system are noise-free variables, i.e., without (or with ignorable) fluctuations, which implies that the covariance of any two participating reactants is zero (or close to zero).

Each of these assumptions constitutes a sufficient condition for the direct use of a rate function as the propensity function and applies, in principle, to GMA as well as other systems. The validity of these conditions will be discussed later. Specifically, the first condition will be addressed in the Results section under the headings "0^th-order reaction kinetics" and "1^st-order reaction kinetics, " while the second condition will be discussed under the heading "Real-valued order monomolecular reaction kinetics." The third condition will be the focus of Equations (29-36) and their associated explanations.

In reality, the rates of reactions in biochemical systems are commonly nonlinear functions of the reactant species, and fluctuations within each species are not necessarily ignorable. Therefore, to the valid use of an equation-based model in a stochastic simulation mandates that we know how to define a proper propensity function. The following section addresses this issue. It uses statistical techniques to characterize estimates for both the mean and variance of the propensity function, and these features will allow an assessment of the validity of the assumption α(X) = f _s (X) and prescribe adjustments if the assumption is not valid.

Methods

Deriving the mean and variance of a power-law function of random variables

Consider a generic power-law function of random variables X _s with the format $PL\left(\mathbf{X}\right)=k\prod _{s=1}^{N_s}{X}_s^{f_s}$. Estimates of its mean μ _PL and variance σ _PL are given as

$${\mu }_{PL}\approx k\prod _{s=1}^{N_s}{\mu }_s^{f_s}exp\left(\sum _{i<j}^{N_s}{f}_i{f}_j\mathsf{\text{cov}}\left[log{X}_i,log{X}_j\right]\right)$$

(8)

$${{\sigma }_{PL}}^2\approx {{\mu }_{PL}}^2\mathrm{\Omega }$$

(9)

(for details, see Additional file 1). Here,

$$\mathrm{\Omega }=\sum _{s=1}^{N_s}{f}_s{\mu }_s^{-2}{\sigma }_s^2+2\sum _{i<j}^{N_s}{f}_i{f}_j\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right]$$

(10)

and μ _s = E[X _s ] and ${\sigma }_s^2=\mathsf{\text{E}}\left[{\left(X_s-{\mu }_s\right)}^2\right]$ are the mean and variance of random variable X _s , respectively. If we choose to express cov [logX _i , logX _j ] as a function of μ _s , σ _s ² and covariance σ _ij = cov [X _i , X _j ], using a Taylor approximation, we obtain

$${\mu }_{PL}\approx k\prod _{s=1}^{N_s}{\mu }_s^{f_s}exp\left(-\frac{1}{2}\sum _{s=1}^{N_s}{f}_s{\sigma }_s^2/{\mu }_s^2+\frac{1}{2}\mathrm{\Omega }\right)$$

(11)

$${{\sigma }_{PL}}^2\approx {{\mu }_{PL}}^2\mathrm{\Omega },$$

(12)

where

$$\begin{array}{c}\mathrm{\Omega }\approx \sum _{s=1}^{N_s}{f}_s{\left({\sigma }_s/{\mu }_s\right)}^2+2\sum _{i<j}^{N_s}{f}_i{f}_j\left\{{\sigma }_{ij}/\left({\mu }_i{\mu }_j\right)\right.\\ +\frac{1}{2}log\left({\mu }_i\right){\left({\sigma }_j/{\mu }_j\right)}^2+\frac{1}{2}log\left({\mu }_j\right){\left({\sigma }_i/{\mu }_i\right)}^2\\ \left.-\frac{1}{4}{\left({\sigma }_i/{\mu }_i\right)}^2{\left({\sigma }_j/{\mu }_j\right)}^2\right\}.\end{array}$$

(13)

Since many biochemical variables approximately follow a log-normal distribution [20–22], it is valuable to consider the special situation where (X₁, ..., X _s )is log-normally distributed (i.e., (logX₁, ..., logX _s ) is normally distributed). In such a case, a simpler alternative way to calculate cov [logX _i , logX _j ] is

$$\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right]=log\left(1+\frac{{\sigma }_{ij}}{{\mu }_i{\mu }_j}\right).$$

(14)

[23]. By substituting this result into (8)-(10), one obtains

$${\mu }_{PL}\approx k\prod _{s=1}^{N_s}{\mu }_s^{f_s}\prod _{i<j}^{N_s}{\left(1+\frac{{\sigma }_{ij}}{{\mu }_i{\mu }_j}\right)}^{f_i{f}_j}$$

(15)

$${{\sigma }_{PL}}^2\approx {{\mu }_{PL}}^2\mathrm{\Omega },$$

(16)

where

$$\mathrm{\Omega }=\sum _{s=1}^{N_s}{f}_s{\left(\frac{{\sigma }_s}{{\mu }_s}\right)}^2+2\sum _{i<j}^{N_s}{f}_i{f}_{j}log\left(1+\frac{{\sigma }_{ij}}{{\mu }_i{\mu }_j}\right).$$

(17)

The approximation formulae for μ _PL and σ _PL ² in eqns. (8)-(10) provide an easy numerical implementation if observation data are available to estimate cov [logX _i , logX _j ]. Furthermore, Equations (11)-(13) demonstrate how μ _PL and σ _PL ² are related to μ _s , σ _s ² and σ _ij ; however, the price of this insight is paid by the possible inaccuracy introduced through the Taylor approximation. Equations (15)-(17) also provide a functional dependence of μ _PL and σ _PL ² on (μ _s , σ _s ², σ _ij ), but it is only valid if the additional assumption of log-normality is acceptable.

Deriving proper propensity functions for stochastic simulations from differential equation-based models

Assuming that the GMA model faithfully captures the average behaviour of a biochemical reaction system and recalling $\left[\mathbf{X}\left(t\right)\right]={\left(\left[X_1\left(t\right)\right],\dots ,\left[X_{N_s}\left(t\right)\right]\right)}^T$, the expected metabolite numbers are defined as the expectation

$$E\left[\mathbf{X}\right]=\left[\mathbf{X}\right]\mathrm{\Phi },$$

(18)

where Φ is the system size as defined above.

To describe the reaction channel R _r stochastically, one needs the state update vector v _r and must characterize the quantity of molecules flowing through of reaction channel R _r during a small time interval. The key concept of this type of description is the propensity function α _r (x), which is defined as

$$\begin{array}{c}{\alpha }_r\left(\mathbf{x}\right)dt=\mathsf{\text{the}}\mathsf{\text{probability}}\mathsf{\text{that}}\mathsf{\text{exactly}}\mathsf{\text{one}}\mathsf{\text{reaction}}\\ R_r\mathsf{\text{will}}\mathsf{\text{occur}}\mathsf{\text{some}}\mathsf{\text{where}}\mathsf{\text{inside}}U\mathsf{\text{within}}\mathsf{\text{infinitesmal}}\\ \mathsf{\text{interval}}\left(t,t+dt\right),\mathsf{\text{given}}\mathsf{\text{current}}\mathsf{\text{state}}\mathbf{X}\left(t\right)=\mathbf{x}.\end{array}$$

(19)

[1]. Because of the probabilistic nature of the propensity function, X(t) is no longer deterministic, and the result is instead stochastic and based on the transition probability

$$P\left(\mathbf{x},t|{\mathbf{x}}_0,t_0\right)=\mathsf{\text{Prob}}\left\{\mathbf{X}\left(t\right)=\mathbf{x},\mathsf{\text{given}}\mathbf{X}\left(t_0\right)={\mathbf{x}}_0\right\},$$

(20)

which follows the chemical master equation (CME)

$$\begin{array}{cc}\hfill \frac{\partial P\left(\mathbf{x},t|{\mathbf{x}}_0,t_0\right)}{\partial t}& =\sum _{r=1}^{N_r}[{\alpha }_r\left(\mathbf{x}+{\mathbf{v}}_r\right)P\left(\mathbf{x}+{\mathbf{v}}_r,t|{\mathbf{x}}_0,t_0\right)\hfill \\ -{\alpha }_r\left(\mathbf{x}\right)P\left(\mathbf{x},t|{\mathbf{x}}_0,t_0\right)]\hfill \end{array}$$

(21)

Updating CME requires knowledge of every possible combination of all species counts within the population, which immediately implies that it can be solved analytically for only a few very simple systems and that numerical solutions are usually prohibitively expensive [24]. To address the inherent intractability of CME, Gillespie developed an algorithm, called the Stochastic Simulation Algorithm (SSA), to simulate CME models [1]. SSA is an exact procedure for numerically simulating the time evolution of a well-stirred reaction system. It is rigorously based on the same microphysical premise that underlies CME and gives a more realistic representation of a system's evolution than a deterministic reaction rate equation represented by ODEs. SSA requires knowledge of the propensity function, which however is truly available only for elementary reactions. These reactions include: 1) 0^th order reactions, exemplified with the generation of a molecule at a constant rate; 2) 1^st order monomolecular reactions, such as an elemental chemical conversion or decay of a single molecule; 3) 2^nd order bimolecular reactions, including reactive collisions between two molecules of the same or different species. The reactive collision of more than two molecules at exactly the same time is considered highly unlikely and modelled as two or more sequential bimolecular reactions.

For elementary reactions, the propensity function of reaction R _r is computed as the product of a stochastic rate constant c _r and the number h _r of distinct combinations of reactant molecules, i.e.

$${\alpha }_r\left(\mathbf{x}\right)=c_r{h}_r\left(\mathbf{x}\right),r=1,\dots ,N_r.$$

(22)

Here $h_r\left(\mathbf{x}\right)=\left\{\begin{array}{c}\hfill \prod _{s=1}^{N_s}\left(\begin{array}{c}\hfill x_s\hfill \\ \hfill {\underset{-}{v}}_{rs}\hfill \end{array}\right)\approx \frac{\prod _{s=1}^{N_s}{x}_s^{{\underset{-}{v}}_{rs}}}{\prod _{s=1}^{N_s}{\underset{-}{v}}_{rs}!},\mathsf{\text{for}}{x}_s\ge {\underset{-}{v}}_{rs}>0\hfill \\ \hfill 0,\mathsf{\text{otherwise}}\hfill \end{array}\right.$, where x _s is the sample value of random variable X _s . The approximation is invoked when x _s is large and (x _s - 1), ..., (x _s - v _rs + 1) are approximately equal to x _s .

In Gillespie's original formulation [1] c _r is a constant that only depends on the physical properties of the reactant molecules and the temperature of the system, and c _r dt is the probability that a particular combination of reactant molecules will react within the next infinitesimally small time interval (t, t + dt). The constant c _r can be calculated from the corresponding deterministic rate constants, if they are known.

Since the assumption of mass action kinetics is not valid generally, especially in spatially restricted environments and in situations dominated by macromolecular crowding, we address the broader scenario where c _r is not a constant but a function of the reactant concentrations. Thus, we denote c _r as a stochastic rate function, while retaining the definition of h _r as above. Knowing that any positive-valued differentiable function can be approximated locally by a power-law function, we assume the functional form of the stochastic rate function as

$$c_r\left(\mathbf{x}\right)={\kappa }_r\prod _{s=1}^{N_s}{x}_s{\left(t\right)}^{{\epsilon }_{rs}}.$$

(23)

Here, κ _r and ε _rs are constants that will be specified in the next section, and r = 1, ..., N _r . Note that ε _rs are now real-valued. Once the stochastic rate function is determined (see below), the propensity function can be calculated as

$${\alpha }_r\left(\mathbf{x}\right)=c_r\left(\mathbf{x}\right)h_r\left(\mathbf{x}\right)=\frac{{\kappa }_r}{\prod _{s=1}^{N_s}{\underset{-}{v}}_{rs}!}\prod _{s=1}^{N_s}{x}_s^{{\underset{-}{v}}_{rs}+{\epsilon }_{rs}}.$$

(24)

In order to identify the functional expression for a stochastic rate function, and thus the propensity function, we consider the connection between the stochastic and the deterministic equation models. By multiplying CME with x and summing over all x, we obtain

$$\frac{d}{dt}E\left[\mathbf{X}\left(t\right)\right]=\sum _{r=1}^{N_r}{\mathbf{v}}_{r}E\left[{\alpha }_r\left(\mathbf{X}\left(t\right)\right)\right].$$

(25)

Similarly, the expectation for any species X _s (t) is given as

$$\frac{d}{dt}E\left[X_s\left(t\right)\right]=\sum _{r=1}^{N_r}{v}_{rs}E\left[{\alpha }_r\left(\mathbf{X}\left(t\right)\right)\right],s=1,\dots ,N_s.$$

(26)

The details of these derivations are shown in Additional file 1.

We can use these results directly to compute the propensity function for a stochastic GMA model, assuming that its deterministic counterpart is well defined. Specifically, we start with the deterministic GMA equation for X _s ,

$$\frac{d}{dt}\left[X_s\left(t\right)\right]=\sum _{r=1}^{N_r}{v}_{rs}{k}_r\prod _{s^{\prime }=1}^{N_s}{\left[X_{s^{\prime }}\left(t\right)\right]}^{f_{r{s}^{\prime }}},s=1,\dots ,N_s,$$

(27)

where v _rs , k _r and f_rs' are again the stoichiometric coefficients, rate constants, and kinetic orders, respectively. By substituting $\left[X_s\right]=\frac{E\left[X_s\right]}{\mathrm{\Phi }}$ from Equation (18) into this GMA model, we obtain a "particle-based" equation of the format

$$\frac{d}{dt}\left(\frac{E\left[X_s\right]}{\mathrm{\Phi }}\right)=\sum _{r=1}^{N_r}{v}_{rs}{k}_r\prod _{s^{\prime }=1}^{N_s}{\left(\frac{E\left[X_{s^{\prime }}\right]}{\mathrm{\Phi }}\right)}^{f_{r{s}^{\prime }}},s=1,\dots ,N_s.$$

Elementary operations allow us to rewrite this equation as

$$\frac{d}{dt}\left(E\left[X_s\right]\right)=\sum _{r=1}^{N_r}{v}_{rs}{k}_r{\mathrm{\Phi }}^{1-F_r}\prod _{s^{\prime }=1}^{N_s}E{\left[X_{s^{\prime }}\right]}^{f_{r{s}^{\prime }}},s=1,\dots ,N_s,$$

(28)

where $F_r=\sum _{s^{\prime }=1}^{N_s}{f}_{r{s}^{\prime }}$. In this formulation, the differential operator is justified only when large numbers of molecules are involved. The assumption that the deterministic equations precisely capture the average behaviour of the biochemical reaction system directly equates the stochastic CME (25) to the deterministic equation based model (28)

$$E\left[{\alpha }_r\left(\mathbf{X}\left(t\right)\right)\right]=k_r{\mathrm{\Phi }}^{1-F_r}\prod _{s^{\prime }=1}^{N_s}E{\left[X_{s^{\prime }}\right]}^{f_{r{s}^{\prime }}}.$$

(29)

Now we have two choices for approximating the expectation of the propensity function on left-hand side of equation (29):

• 1) adopt a zero-covariance assumption as was done in [25], which implies ignoring random fluctuations within every species as well as their correlations. This assumption is only justified for some special cases such as monomolecular and bimolecular reactions under the thermodynamic limit (cf. [4, 6]), but is not necessary valid in generality. Here the thermodynamic limit is defined as a finite concentration limit which the system reaches when both population and volume approach infinity. Under this assumption, the left hand side of (29) becomes

  $$\begin{array}{cc}\hfill E\left[{\alpha }_r\left(\mathbf{x}\right)\right]& =E\left[\frac{{\kappa }_r}{\prod _{s=1}^{N_s}{\underset{-}{v}}_{rs}!}\prod _{s=1}^{N_s}{x}_s^{{\underset{-}{v}}_{rs}+{\epsilon }_{rs}}\right]\hfill \\ =\frac{{\kappa }_r}{\prod _{s=1}^{N_s}{\underset{-}{v}}_{rs}!}\prod _{s=1}^{N_s}E{\left[X_s\right]}^{{\underset{-}{v}}_{rs}+{\epsilon }_{rs}}\hfill \end{array}$$

  (30)

for every r = 1, ..., N _r , and Equation (24) yields

$$\begin{array}{cc}\hfill {\epsilon }_{rs}& =f_{rs}-{\underset{-}{v}}_{rs}\hfill \\ \hfill {\kappa }_r& =k_r{\mathrm{\Phi }}^{1-F_r}\prod _{s=1}^{N_s}{\underset{-}{v}}_{rs}!\hfill \\ \hfill c_r\left(\mathbf{x}\right)& =k_r{\mathrm{\Phi }}^{1-F_r}\prod _{s=1}^{N_s}{\underset{-}{v}}_{rs}!x_s^{{\epsilon }_{rs}}\hfill \end{array}$$

(31)

and

$${\alpha }_{r_{-}0}\left(\mathbf{x}\right)=k_r{\mathrm{\Phi }}^{1-F_r}\prod _{s=1}^{N_s}{x}_s^{f_{rs}}.$$

(32)

Here, the index r_0 is used to distinguish this 0-covariance propensity function from a second type of propensity in the next section.

With the zero-covariance assumption, one can substitute (32) back into the equation for the expectation for each species, which yields

$$\frac{d}{dt}E\left[X_s\left(t\right)\right]=\sum _{r=1}^{N_r}{v}_{rs}{k}_r{\mathrm{\Phi }}^{1-F_s}\prod _{s=1}^{N_s}{\mu }_s^{f_{rs}}$$

(33)

for every s = 1, ..., N _s .. Note that this result is exactly equivalent to the equation-based model (27).

Equation (33) is based on assumption that both the fluctuations within species and their correlations are ignorable, which is not necessarily true in reality. If one uses it in simulations where the assumptions are not satisfied, it is possible that the means for the molecular species are significantly different from the corresponding equation-based model values. This discrepancy arises because the evolution of each species in the stochastic simulation is in truth affected by the covariance which is not necessarily zero, as it was assumed. This phenomenon was observed by Paulsson and collaborators [26] and further discussed in different moment-based approaches [6, 7]. To assess the applicability limit of the propensity defined by (32), we can apply approximation techniques as shown in eqns. (8)-(10) on the functional expression of α_{r_0} and obtain mean and variance as

$$\begin{array}{cc}\hfill {\mu }_{{\alpha }_{r_{-}0}}& =E\left[{\alpha }_{r_{-}0}\left(\mathbf{X}\left(t\right)\right)\right]\hfill \\ \approx k_r{\mathrm{\Phi }}^{1-F_r}\prod _{s^{\prime }=1}^{N_s}E{\left[X_{s^{\prime }}\right]}^{f_{r{s}^{\prime }}}exp\left(\sum _{i<j}^{N_s}{f}_{ri}{f}_{rj}\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right]\right)\hfill \end{array}$$

(34)

$${{\sigma }_{{\alpha }_{r_{-}0}}}^2\approx {{\mu }_{{\alpha }_{r_{-}0}}}^2{\mathrm{\Omega }}_r,$$

(35)

where

$${\mathrm{\Omega }}_r=\sum _{s=1}^{N_s}{f}_{rs}{\mu }_s^{-2}{\sigma }_s^2+2\sum _{i<j}^{N_s}{f}_{ri}{f}_{rj}\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right],$$

(36)

for every s = 1, ..., N _s . These expressions demonstrate that even with large numbers of molecules the mean of CME does not always converge to the GMA model. Indeed, the convergence is only guaranteed in one of the following special situations: 1) the reaction is of 0^th order; 2) the reaction is a real value-order monomolecular reaction, with 1^st order reaction as a special case; 3) the covariance contribution in (34) is sufficiently small to be ignored for all participating reactant species of a particular reaction channel. Except for these three special situations, the covariance as shown in (34) significantly affects the mean dynamics. Therefore, stochastic simulations using zero-covariance propensity functions will in general yield means different from what the deterministic GMA model produces. How large these differences are cannot be said in generality. Under the assumption that the GMA model correctly captures the mean dynamics of every species, this conclusion means that α_{r_0} is not necessarily an accurate propensity function for stochastic simulations, and the direct conversion of the equation-based model into a propensity function must be considered with caution.

Moreover, there is no theoretical basis to assume that there are no fluctuations in the molecular species or that these are independent. Therefore, we need to consider the second treatment of the expectation of the propensity function and study the possible effects of a non-zero covariance.

• 2) We again assume that the GMA model is well defined, which implies that information regarding the species correlations and fluctuations has been captured in the parameters of the GMA model on the left hand size of Equations (7) and (28). To gain information regarding correlations, we use Taylor expansion to approximate the propensity function (see Additional file 1 for details):

  $$\begin{array}{c}E\left[{\alpha }_r(X(t))\right]=E\left[\frac{{\kappa }_r}{{\displaystyle \prod _{s=1}^{N_s}{\underset{¯}{v}}_{rs}}!}{\displaystyle \prod _{s=1}^{N_s}{X}_s^{{\underset{¯}{v}}_{rs}}{}^{+{\epsilon }_{rs}}}\right]\\ \approx \frac{{\kappa }_r}{{\displaystyle \prod _{s=1}^{N_s}{\underset{¯}{v}}_{rs}!}}{\displaystyle \prod _{s=1}^{N_s}E{\left[X_s\right]}^{{\underset{¯}{v}}_{rs}+{\epsilon }_{rs}}}\\ \times \exp \left({\displaystyle \sum _{i<j}^{N_s}\left({\underset{¯}{v}}_{ri}+{\epsilon }_{ri}\right)\left({\underset{¯}{v}}_{rj}+{\epsilon }_{rj}\right)}\text{cov}\left[1\text{og}{X}_i,\log X_j\right]\right)\end{array}$$

  (37)

After substitution of (37) in (29), one obtains

$$\begin{array}{cc}\hfill {\kappa }_r& =k_r{\mathrm{\Phi }}^{1-F_r}\prod _{s=1}^{N_s}{\underset{-}{v}}_{rs}!exp\left(-\sum _{i<j}^{N_s}{f}_{ri}{f}_{rj}\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right]\right)\hfill \\ \hfill {\epsilon }_{rs}& =f_{rs}-{\underset{-}{v}}_{rs}.\hfill \end{array}$$

Given the state x of the system at time t, the stochastic rate function of reaction R _r is

$$\begin{array}{cc}\hfill c_r\left(\mathbf{x}\right)& ={\kappa }_r\prod _{s=1}^{N_s}{x}_s^{{\epsilon }_{rs}}\hfill \\ =k_r{\mathrm{\Phi }}^{1-F_r}\prod _{s=1}^{N_s}{\underset{-}{v}}_{rs}!\hfill \\ \times exp\left(-\sum _{i<j}^{N_s}{f}_{ri}{f}_{rj}\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i\left(t\right),log{X}_j\left(t\right)\right]\right)\prod _{s=1}^{N_s}{x}_s^{f_{rs}-{\underset{-}{v}}_{rs}}.\hfill \end{array}$$

(38)

Here it is important to understand that although the random variables {X _s }_s∈Sappear in the expression c _r (x), c _r (x) is not a function of random variables but a deterministic function. The reason is that the cov [logX _i (t), logX _j (t)] in the composition of c _r (x), which as the numerical characteristic of the random variables {X _s }_s∈S, is deterministic. Therefore, the stochastic rate function c _r (x) is a well-justified deterministic function that is affected by both the state of the system $\left[x_1,\dots ,x_{N_s}\right]$ and cov [logX _i (t), logX _j (t)], the numerical characteristic of fluctuations in the random variables {X _s }_s∈S.

Given the expression c _r (x), the propensity function is

$$\begin{array}{cc}\hfill {\alpha }_r\left(\mathbf{x}\right)& =c_r\left(\mathbf{x}\right)h_r\left(\mathbf{x}\right)\hfill \\ =k_r{\mathrm{\Phi }}^{1-F_r}\prod _{s=1}^{N_s}{x}_s^{f_{rs}}\hfill \\ \times exp\left(-\sum _{i<j}^{N_s}{f}_{ri}{f}_{rj}\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i\left(t\right),log{X}_j\left(t\right)\right]\right).\hfill \end{array}$$

(39)

These results are based on the assumption that there are large numbers of molecules for all reactant species participating in reaction R _r . For simplicity of discussion, we define the propensity adjustment factor (paf) of reaction R _r as

$$paf\left(t\right)\triangleq exp\left(-\sum _{i<j}^{N_s}{f}_{ri}{f}_{rj}\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i\left(t\right),log{X}_j\left(t\right)\right]\right).$$

(40)

paf is a function of time t and represents the contribution of the reactants to correlations among species in the calculation of the propensity function for reaction R _r . We denote the propensity function in (39), which accounts for the contribution of the covariance, as α_{r_cov}, in order to distinguish it from the propensity function α_{r_0} (32), which is based on the assumption of zero-covariance, i.e.,

$${\alpha }_{r_{-}\mathsf{\text{cov}}}\left(\mathbf{x}\right)=paf\left(t\right)k_r{\mathrm{\Phi }}^{1-F_r}\prod _{s=1}^{N_s}{x}_s^{f_{rs}}.$$

(41)

Remembering that cov [logX _i (t), logX _j (t)], which is a component in both the stochastic rate function c _r (x) and now in the function paf(t), is a deterministic function rather than a function of random variables, paf(t) is a deterministic correction to the kinetic constant k _r in the construction of α_{r_cov} in (41), which corrects the stochastic simulation toward the correct average.

In contrast to the propensity function α_{r_0}, α_{r_cov} leads to accurate stochastic simulations. To illustrate this difference, we analyze $\frac{d}{dt}E\left[X_s\left(t\right)\right]$ as follows: We apply the approximation techniques in eqns. (9)-(11) in order to obtain the mean and variance of the propensity function α_{r_cov}:

$${\mu }_{{\alpha }_{r_{-}\mathsf{\text{cov}}}}=E\left[{\alpha }_{r_{-}\mathsf{\text{cov}}}\left(\mathbf{X}\left(t\right)\right)\right]\approx k_r{\mathrm{\Phi }}^{1-F_r}\prod _{s^{\prime }=1}^{N_s}E{\left[X_{s^{\prime }}\right]}^{f_{r{s}^{\prime }}}$$

(42)

$${{\sigma }_{r_{-}\mathsf{\text{cov}}}}^2\approx {{\mu }_{{\alpha }_{r_{-}\mathsf{\text{cov}}}}}^2{\mathrm{\Omega }}_r.$$

(43)

Here

$${\mathrm{\Omega }}_r=\sum _{s=1}^{N_s}{f}_{rs}{\mu }_s^{-2}{\sigma }_s^2+2\sum _{i<j}^{N_s}{f}_{ri}{f}_{rj}\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right].$$

(44)

By substituting (42) back into the derivation of CME (26), one obtains

$$\begin{array}{c}\frac{d}{dt}E\left[X_s\left(t\right)\right]\\ =\sum _{r=1}^{N_r}{v}_{rs}E\left[{\alpha }_{r_{-}\mathsf{\text{cov}}}\left(\mathbf{X}\left(t\right)\right)\right]\\ \approx \sum _{r=1}^{N_r}{v}_{rs}{k}_r{\mathrm{\Phi }}^{1-F_r}\prod _{s^{\prime }=1}^{N_s}{\mu }_{s^{\prime }}^{f_{r{s}^{\prime }}}\end{array}$$

(45)

for every s = 1, ..., N _s , which is equivalent in approximation to the GMA model (28). In the other words, the mean of every molecular species obtained by using α_{r_cov} in the CME derived equation (27) is approximately identical to the corresponding macroscopic variable in the GMA model.

Calculation of cov [logX _i (t), logX _j (t)]

When data in the form of multiple time series for all the reactants are available, it is possible to compute cov [logX _i (t), logX _j (t)] directly from these data. Once this covariance is known, the function paf, α_{r_cov} and the mean dynamics can all be assessed. Alas, the availability of several time series data for all reactants under comparable conditions is rare, so that cov [logX _i (t), logX _j (t)] must be estimated in a different manner.

If one can validly assume that the covariance based on α_{r_0} does not differ significantly from the covariance based on α_{r_cov}, one may calculate cov [logX _i (t), logX _j (t)] by one of following methods.

Method 1:

One uses α_{r_0} to generate multiple sets of time series data of all reactants and then computes cov [logX _i (t), logX _j (t)].

Method 2:

First, cov [logX _i (t), logX _j (t)] is expressed as a function of mean and covariance in one of the following ways; either as

$$\begin{array}{c}\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right]\approx {\sigma }_{ij}/\left({\mu }_i{\mu }_j\right)+\frac{1}{2}log\left({\mu }_i\right){\left({\sigma }_j/{\mu }_j\right)}^2\\ +\frac{1}{2}log\left({\mu }_j\right){\left({\sigma }_i/{\mu }_i\right)}^2-\frac{1}{4}{\left({\sigma }_i/{\mu }_i\right)}^2{\left({\sigma }_j/{\mu }_j\right)}^2\end{array}$$

(46)

or as Equation (14):

$$\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right]=log\left(1+\frac{{\sigma }_{ij}}{{\mu }_i{\mu }_j}\right).$$

The first functional expression of cov [logX _i (t), logX _j (t)] is achieved by Taylor approximation, whereas the second expression is obtained by the additional assumption that the concentrations (X₁, ..., X _s ) are log-normally distributed [8, 23]. The consideration of a log-normal distribution is often justified by the fact that many biochemical data have indeed been observed to be log-normally distributed (e.g., [20–22]).

Second, one uses α_{r_0} to approximate the mean and covariance either by direct simulation, as shown in method 1, or by a moment-based approach, which is explained in Additional file 2, and which yields the differential equations

$$\begin{array}{cc}\hfill \frac{\partial {\mu }_s}{\partial t}& \approx \sum _{r=1}^{N_r}{v}_{r,s}\left\{{\alpha }_{r_{-}0}\left(\mathbf{\mu }\right)+\frac{1}{2}\sum _{m,n=1}^{N_s}\frac{{\partial }^2{\alpha }_{r\text{\_}0}\left(\mathbf{\mu }\right)}{\partial X_m\partial X_n}{\sigma }_{mn}\right\}\hfill \\ \hfill \frac{\partial {\sigma }_{ij}}{\partial t}& \approx \sum _{r=1}^{N_r}\left\{v_{r,i}\sum _{s=1}^{N_s}\frac{\partial {\alpha }_{r\text{\_}0}\left(\mathbf{\mu }\right)}{\partial X_s}{\sigma }_{js}+v_{r,j}\sum _{s=1}^{N_s}\frac{\partial {\alpha }_{r\text{\_}0}\left(\mathbf{\mu }\right)}{\partial X_s}{\sigma }_{is}\right.\hfill \\ \left.+v_{r,i}{v}_{r,j}\left[{\alpha }_{r_{-}0}\left(\mathbf{\mu }\right)+\frac{1}{2}\sum _{m,n=1}^{N_s}\frac{{\partial }^2{\alpha }_{r\text{\_}0}\left(\mathbf{\mu }\right)}{\partial X_m\partial X_n}{\sigma }_{mn}\right]\right\}\hfill \end{array}$$

For convenience of computational implementation, the above equations can be written in matrix format

$$\begin{array}{cc}\hfill \frac{\partial \mu }{\partial t}& \approx V^T\left(\alpha +\frac{1}{2}{\alpha }^{″}\odot \sigma \right)\hfill \\ \hfill \frac{\partial \sigma }{\partial t}& \approx {\left(\sigma \left({\alpha }^{\prime }V\right)\right)}^T+\sigma \left({\alpha }^{\prime }V\right)+V^T\mathrm{\Lambda }V.\hfill \end{array}$$

Here for r = 1, ..., N _r , and s, m, n = 1, ..., N _s , $\mu ={\left({\mu }_1,\dots ,{\mu }_{N_s}\right)}^T$, (V) _rs = v _rs , $\alpha ={\left({\alpha }_1,...,{\alpha }_{N_r}\right)}^T$, ${\alpha }^{″}={\left({\alpha }_r^{{}^{″}},...,{\alpha }_{N_r}^{{}^{″}}\right)}^T$, ${\left({\alpha }_r^{{}^{″}}\right)}_{mn}=\frac{{\partial }^2{\alpha }_r\left(\mathbf{X}\right)}{\partial X_m\partial X_n}$, ${\alpha }_r^{{}^{″}}\odot \sigma =\sum _{m,n=1}^{N_s}\frac{{\partial }^2{\alpha }_r\left(\mathbf{X}\right)}{\partial X_m\partial X_n}{|}_{\mathbf{X}=\mu }{\sigma }_{mn}$, ${\alpha }^{{}^{″}}\odot \sigma \triangleq {\left({\alpha }_1^{{}^{″}}\odot \sigma ,\dots ,{\alpha }_{N_r}^{{}^{″}}\odot \sigma \right)}^T$, ${\alpha }^{\prime }=\left({{\alpha }_1}^{\prime },\dots ,{{\alpha }_{N_r}}^{\prime }\right)$, ${{\alpha }_r}^{\prime }={\left(\frac{\partial {\alpha }_r\left(\mathbf{\mu }\right)}{\partial X_1},\dots ,\frac{\partial {\alpha }_r\left(\mathbf{\mu }\right)}{\partial X_{N_s}}\right)}^T$, and Λ is a diagonal matrix with ${\left(\mathrm{\Lambda }\right)}_{rr}={\alpha }_r\left(\mu \right)+\frac{1}{2}\sum _{m,n=1}^{N_s}\frac{{\partial }^2{\alpha }_r\left(\mu \right)}{\partial X_m\partial X_n}{\sigma }_{mn}$.

Statistical criteria for propensity adjustment

Suppose an equation-based model captures the average behavior of a stochastic system and one intends to find the propensity function for a stochastic simulation that will reproduce that means. One can use the 95% confidence interval to evaluate the need for a propensity adjustment. Specifically, for stable systems that will reach a steady state, we use the reversible reaction model as an example. If the steady state of the ODE x _st is within the 95% confidence interval of n runs of stochastic simulations, i.e. $x_{st}\in \left[{\mu }_{st}-1.96\frac{{\delta }_{st}}{\sqrt{n}},{\mu }_{st}+1.96\frac{{\delta }_{st}}{\sqrt{n}}\right]$, then the rate function in the original ODEs can be used as the propensity without adjustment; otherwise propensity adjustment is needed. Here μ _st and δ _st can be attained from either a moment-base method or from n independent runs of stochastic simulations using propensity without adjustment. An example discussing a reversible reaction with feedback controls can be found in the results section.

For other systems that do not reach a steady state, but where instead transient characteristics are of the highest interest, one can judge the need of propensity adjustment by whether the pertinent characteristics of the ODEs are within the 95% confidence interval of the corresponding characteristic, which is given by a prediction from the moment-based method or from n runs of stochastic simulations. The Repressilator example in the result section will serve as a demonstration.

Results

Generic special cases

It is generally not valid to translate a rate from a deterministic biochemical model into a propensity function of the corresponding stochastic simulation without adjustment (see Equations. (34)-(36)). However, in some situations, the propensity adjustment (e.g., Equations (40)-(44)) is not needed, and in some other cases it becomes relatively simple.

• 1) 0^th-order reaction kinetics

Consider a very simple equation-based model of the type

$$\frac{d\left[X_s\left(t\right)\right]}{dt}=k_r\mathsf{\text{or}}\frac{dE\left[X_s\left(t\right)\right]}{dt}=k_r\mathrm{\Phi },$$

(47)

for all s = 1, ..., N _s , f _rs = 0. According to Equations (40)-(44), one obtains

$$\begin{array}{c}{\mathrm{\Omega }}_r=0\\ {{\sigma }_{{\alpha }_r}}^2\approx 0\\ {\mu }_{{\alpha }_r}\approx exp\left(log\left(k_r\mathrm{\Phi }\right)\right)=k_r\mathrm{\Phi }\\ \mathsf{\text{i}}\mathsf{\text{.e}}\mathsf{\text{.}}E\left[{\alpha }_r\left(\mathbf{X}\right)\right]\approx {\alpha }_r\left(E\left[\mathbf{X}\right]\right)\\ {\alpha }_{r\text{\_}\mathsf{\text{cov}}}\approx k_r\mathrm{\Phi }={\alpha }_{r\text{\_}0}.\end{array}$$

Thus, for a 0^th-order reaction, its rate equation can be taken directly as the propensity function in stochastic simulations.

• 2) 1^st-order reaction kinetics

Direct application of Equations (40)-(44) yields

$$\begin{array}{cc}\hfill \frac{d\left[X_i\left(t\right)\right]}{dt}& =k_r\left[X_j\left(t\right)\right]\hfill \\ \hfill \mathsf{\text{or}}\frac{dE\left[X_i\left(t\right)\right]}{dt}& =k_{r}E\left[X_j\left(t\right)\right],\hfill \end{array}$$

(48)

f _rs = δ _sj , i, j = 1, ..., N _s . Therefore, according to Equations (40)-(44)

$$\begin{array}{c}{\mathrm{\Omega }}_r={\left({\sigma }_j/{\mu }_j\right)}^2\\ {\left({\sigma }_{{\alpha }_r}/{\mu }_{{\alpha }_r}\right)}^2={\left({\sigma }_j/{\mu }_j\right)}^2\\ {\mu }_{{\alpha }_r}\approx exp\left(log\left(k_r{\mu }_j\right)\right)=k_r{\mu }_j\\ \mathsf{\text{i}}\mathsf{\text{.e}}\mathsf{\text{.}}E\left[{\alpha }_r\left(\mathbf{X}\right)\right]\approx {\alpha }_r\left(E\left[\mathbf{X}\right]\right)\\ {\alpha }_{r_{-}\mathsf{\text{cov}}}\left(\mathbf{X}\right)\approx k_r{X}_j={\alpha }_{r_{-}0}\left(\mathbf{X}\right).\end{array}$$

Thus, for 1^st-order reactions, the rate equation can again be taken directly as the propensity function in stochastic simulations.

• 3) Real-valued order monomolecular reaction kinetics

Consider a reaction with kinetics of the type

$$\begin{array}{cc}\hfill \frac{d\left[X_i\left(t\right)\right]}{dt}& =k_r{\left[X_j\left(t\right)\right]}^{f_{rj}}\hfill \\ \hfill \mathsf{\text{or}}\frac{dE\left[X_i\left(t\right)\right]}{dt}& =k_r{\mathrm{\Phi }}^{1-f_{rj}}E{\left[X_j\left(t\right)\right]}^{f_{rj}},\hfill \end{array}$$

(49)

f _rj ≠ 0, f _rs = 0, for any s ≠ j, s = 1, ..., N _s . Equations (40)-(44) lead to

$$\begin{array}{c}{\mathrm{\Omega }}_r={\left({\sigma }_j/{\mu }_j\right)}^2\\ {\left({\sigma }_{{\alpha }_r}/{\mu }_{{\alpha }_r}\right)}^2={\left({\sigma }_j/{\mu }_j\right)}^2\\ {\mu }_{{\alpha }_r}\approx k_r{\mathrm{\Phi }}^{1-f_{rj}}{\mu }_j^{f_{rj}}\mathsf{\text{i}}\mathsf{\text{.e}}\mathsf{\text{.}}E\left[{\alpha }_r\left(\mathbf{X}\right)\right]\approx {\alpha }_r\left(E\left[\mathbf{X}\right]\right)\\ {\alpha }_{r_{-}\mathsf{\text{cov}}}\left(\mathbf{X}\right)\approx k_r{\mathrm{\Phi }}^{1-f_{rj}}{X}_j^{f_{rj}}={\alpha }_{r_{-}\mathsf{\text{0}}}\left(\mathbf{X}\right).\end{array}$$

Thus, for reaction kinetics involving a single variable and a real-valued order, the rate equation can again be taken as the propensity function in stochastic simulations.

• 4) 2^nd-order reaction kinetics

This type of reaction can be expressed as

$$\begin{array}{cc}\hfill \frac{d}{dt}\left[X_s\left(t\right)\right]& =k_r\left[X_i\left(t\right)\right]\left[X_j\left(t\right)\right]\hfill \\ \hfill \mathsf{\text{or}}\frac{dE\left[X_s\left(t\right)\right]}{dt}& =k_r{\mathrm{\Phi }}^{-1}E\left[X_i\left(t\right)\right]E\left[X_j\left(t\right)\right],\hfill \end{array}$$

(50)

i, j ∈ {1, ..., N _s }, i ≠ j, f _ri = f _rj = 1, and f _rs = 0, for all s ≠ i, j. Therefore, according to Equations (40)-(44)

$$\begin{array}{cc}\hfill {\mathrm{\Omega }}_r& ={\left({\sigma }_i/{\mu }_i\right)}^2+{\left({\sigma }_j/{\mu }_j\right)}^2+2\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right]\hfill \\ ={\left({\sigma }_i/{\mu }_i\right)}^2+{\left({\sigma }_s/{\mu }_s\right)}^2+2\left\{\mathsf{\text{cov}}\left[X_i/{\mu }_i,X_j/{\mu }_j\right]\right.\hfill \\ \left.+\frac{1}{2}log\left({\mu }_i\right){\left({\sigma }_j/{\mu }_j\right)}^2+\frac{1}{2}log\left({\mu }_j\right){\left({\sigma }_i/{\mu }_i\right)}^2-\frac{1}{4}{\left({\sigma }_i/{\mu }_i\right)}^2{\left({\sigma }_j/{\mu }_j\right)}^2\right\}\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\left({\sigma }_{{\alpha }_r}/{\mu }_{{\alpha }_r}\right)}^2& ={\mathrm{\Omega }}_r\hfill \\ \hfill \partial & ={\left({\sigma }_i/{\mu }_i\right)}^2+{\left({\sigma }_j/{\mu }_j\right)}^2+2\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right]\hfill \\ \hfill {\mu }_{{\alpha }_r}& \approx k_r{\left(N_{A}V\right)}^{-1}{\mu }_i{\mu }_j\hfill \\ \hfill {\alpha }_{r_{-}\mathsf{\text{cov}}}\left(\mathbf{X}\right)& =k_r{\mathrm{\Phi }}^{-1}{X}_i{X}_{j}exp\left(-\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right]\right)\ne {\alpha }_{r_{-}0}\left(\mathbf{X}\right).\hfill \end{array}$$

Thus, the proper propensity function for 2^nd-order reactions is different from the rate equation. The difference can be ignored only if the contribution from the covariance is insignificant. In general, the rate equation yields only an approximate propensity function for stochastic simulations, and the approximation quality must be assessed on a case-by-case basis.

• 5) Bimolecular reaction with real-valued order kinetics

This type of reaction can be formulated as

$$\begin{array}{cc}\hfill \frac{d\left[X_s\left(t\right)\right]}{dt}& =k_r{\left[X_i\left(t\right)\right]}^{f_i}{\left[X_j\left(t\right)\right]}^{f_j}\hfill \\ \hfill \mathsf{\text{or}}\frac{dE\left[X_s\left(t\right)\right]}{dt}& =k_r{\mathrm{\Phi }}^{1-f_i-f_j}E{\left[X_i\left(t\right)\right]}^{f_i}E{\left[X_j\left(t\right)\right]}^{f_j},\hfill \end{array}$$

(51)

i, j ∈ {1, ..., N _s }, i ≠ j, f _ri , f _rj ≠ 0, and f _rs = 0, for all s ≠ i, j. According to Equations (40)-(44) we obtain

$$\begin{array}{cc}\hfill {\mathrm{\Omega }}_r& ={\left({\sigma }_i/{\mu }_i\right)}^2+{\left({\sigma }_j/{\mu }_j\right)}^2+2f_i{f}_j\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right]\hfill \\ ={\left({\sigma }_i/{\mu }_i\right)}^2+{\left({\sigma }_j/{\mu }_j\right)}^2+2f_i{f}_j\left\{\mathsf{\text{cov}}\left[X_i/{\mu }_i,X_j/{\mu }_j\right]\right.\hfill \\ \left.+\frac{1}{2}log\left({\mu }_i\right){\left({\sigma }_j/{\mu }_j\right)}^2+\frac{1}{2}log\left({\mu }_j\right){\left({\sigma }_i/{\mu }_i\right)}^2-\frac{1}{4}{\left({\sigma }_i/{\mu }_i\right)}^2{\left({\sigma }_j/{\mu }_j\right)}^2\right\}\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\left({\sigma }_{{\alpha }_r}/{\mu }_{{\alpha }_r}\right)}^2& ={\mathrm{\Omega }}_r\hfill \\ ={\left({\sigma }_i/{\mu }_i\right)}^2+{\left({\sigma }_j/{\mu }_j\right)}^2+2f_i{f}_j\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right]\hfill \\ \hfill {\mu }_{{\alpha }_r}& \approx k_r{\mathrm{\Phi }}^{f_i+f_j-1}{\mu }_i^{f_i}{\mu }_j^{f_j}\hfill \\ \hfill {\alpha }_{r_{-}\mathsf{\text{cov}}}\left(\mathbf{X}\right)& =k_r{\mathrm{\Phi }}^{f_i+f_j-1}{X}_i^{f_i}{X}_j^{f_j}exp\left(-f_i{f}_j\mathsf{\text{cov}}\left[1\mathsf{\text{og}}{X}_i,log{X}_j\right]\right)\hfill \\ \ne {\alpha }_{r_{-}0}\left(\mathbf{X}\right).\hfill \end{array}$$

For bimolecular reactions of complex order, the propensity function is different from the rate equation. The difference can be ignored only if the contribution from the covariance is insignificant.

Power-law representation of a reversible reaction with feedback controls

We consider a reversible reaction with feedback controls (see Figure 1) whose average behaviour is accurately described by the following GMA model

Figure 1

Scheme of reversible reaction with feedback controls. S₃ inhibits the forward reaction and S₁ activates the reverse reaction.

$$\begin{array}{c}\frac{d{x}_1}{dt}=\frac{d{x}_2}{dt}=-\frac{d{x}_3}{dt}\\ =-k_f{\mathrm{\Phi }}^{1-f_1-f_2-f_3}{x}_1^{f_1}{x}_2^{f_2}{x}_3^{f_3}+k_b{\mathrm{\Phi }}^{1-g_1-g_3}{x}_1^{g_1}{x}_3^{g_3}.\end{array}$$

(52)

Here S₃ feeds back to inhibit the forward reaction and S₁ feeds back on the reverse reaction and accelerates it. The task is to develop a stochastic model whose performance converges to that of the deterministic GMA model. We can see from equations (52) that three variables x₁, x₂ and x₃ contribute to the forward flux $k_f{\mathrm{\Phi }}^{1-f_1-f_2-f_3}{x}_1^{f_1}{x}_2^{f_2}{x}_3^{f_3}$ and two variables x₁ and x₃ contribute to the backward flux $k_b{\mathrm{\Phi }}^{1-g_1-g_3}{x}_1^{g_1}{x}_3^{g_3}$. Because several variables are involved, their covariance has the potential of affecting the forward and the backward propensity functions in a stochastic simulation. To obtain the covariance information, we formulate the moment equations (53) from the ODE model (52).

To simplify the calculation, as explained in detail in Additional file 2, we set the third central moment to zero and obtain a closed-form set of ODEs. Expressed differently, the rate of change in mean and covariance depends only on the functions of mean and covariance themselves, but not on higher-order moments. Thus,

$$\begin{array}{cc}\hfill \frac{\partial \mu }{\partial t}& \approx V^T\left(\alpha +\frac{1}{2}{\alpha }^{″}\odot \sigma \right)\hfill \\ \hfill \frac{\partial \sigma }{\partial t}& \approx {\left(\sigma \left({\alpha }^{\prime }V\right)\right)}^T+\sigma \left({\alpha }^{\prime }V\right)+V^T\mathrm{\Lambda }V.\hfill \end{array}$$

(53)

Here μ = (μ₁, μ₂, μ₃)^T, $V=\left[\begin{array}{ccc}\hfill -1\hfill & \hfill -1\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 1\hfill & \hfill -1\hfill \end{array}\right]$, $\alpha ={\left({\alpha }_1,{\alpha }_2\right)}^T=\left[\begin{array}{c}\hfill k_f{\mathrm{\Phi }}^{1-f_1-f_2-f_3}{x}_1^{f_1}{x}_2^{f_2}{x}_3^{f_3}\hfill \\ \hfill k_b{\mathrm{\Phi }}^{1-g_1-g_3}{x}_1^{g_1}{x}_3^{g_3}\hfill \end{array}\right]$.

Moreover, for r = 1, 2 and m, n = 1, 2, 3, ${\left({\alpha }_r^{{}^{″}}\right)}_{mn}=\frac{{\partial }^2{\alpha }_r\left(\mathbf{X}\right)}{\partial x_m\partial x_n}$, α" = (α₁", α₂")^T, $\sigma =\left[\begin{array}{ccc}\hfill {\sigma }_{11}\hfill & \hfill {\sigma }_{12}\hfill & \hfill {\sigma }_{13}\hfill \\ \hfill {\sigma }_{21}\hfill & \hfill {\sigma }_{22}\hfill & \hfill {\sigma }_{23}\hfill \\ \hfill {\sigma }_{31}\hfill & \hfill {\sigma }_{32}\hfill & \hfill {\sigma }_{33}\hfill \end{array}\right]$, ${\alpha }_r^{{}^{″}}\odot \sigma \triangleq \sum _{m,n=1}^3\frac{{\partial }^2{\alpha }_r\left(\mathbf{X}\right)}{\partial x_m\partial x_n}{|}_{\mathbf{X}=\mu }{\sigma }_{mn}$, α"⊙σ ≜ (α₁"⊙ σ, α₂"⊙ σ)^T, α' = (α₁', α₂'), ${{\alpha }_r}^{\prime }={\left(\frac{\partial {\alpha }_r\left(\mu \right)}{\partial x_1},\frac{\partial {\alpha }_r\left(\mu \right)}{\partial x_2},\frac{\partial {\alpha }_r\left(\mu \right)}{\partial x_3}\right)}^T$, and $\mathrm{\Lambda }=\left[\begin{array}{cc}\hfill {\alpha }_1\left(\mu \right)+\frac{1}{2}\sum _{m,n=1}^3\frac{{\partial }^2{\alpha }_1\left(\mu \right)}{\partial x_m\partial x_n}{\sigma }_{mn}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\alpha }_2\left(\mu \right)+\frac{1}{2}\sum _{m,n=1}^3\frac{{\partial }^2{\alpha }_2\left(\mu \right)}{\partial x_m\partial x_n}{\sigma }_{mn}\hfill \end{array}\right]$.

Two initial conditions are chosen for representative simulations; they differ by a factor of 20 in species populations and reaction volume between the upper and lower panels of Figure 2. The purpose is to observe the thermodynamic limit of the systems: both scenarios have the same initial concentrations, but the system in the lower panel case has a larger species populations and reaction volume and can thus be regarded as the thermodynamic limit sample of system in the upper panel. As demonstrated by the figures in the first column, the moment approach predicts that for both population sizes the average trajectories of the stochastic model (without propensity adjustment) dynamics is lower than that of the equation-based model: the differences are about 10% of the steady-state value of the equation-based model in the upper figure and 1% in the lower figure; for 100 runs of the stochastic simulation, the steady-state value of the equation-based model lies outside the 95% confidence interval in the upper figure, while it is inside the interval in the lower figure. Therefore, we can expect that the propensity adjustment will significantly contribute to the stochastic simulation for the upper case while not for the lower case. This expectation is confirmed by the simulation results in the third and fourth columns. With the common assumption that the deterministic equations precisely capture the system's average behaviour, the case in the upper panel represents the situation where propensity adjustment is needed, while the lower panel represents the situation that a propensity without adjustment is sufficient when the system approaches its thermodynamics limit. This example furthermore demonstrates that either the moment approach or the stochastic simulations without propensity adjustment can be used to estimate whether there is a need to construct a propensity adjustment function for stochastic simulations.

Figure 2

Comparative simulation results for a reversible reaction with feedback controls. In all panels, the x-axis denotes time in seconds and the y-axis represents the number of molecules of species S₁. The upper and lower panels use two different sets of initial numbers of molecules, namely: (x₁(0), x₂(0), x₃(0), U) = (5, 5, 6, 1μm³) and (x₁(0), x₂(0), x₃(0), U) = (100, 100, 120, 20μm³), respectively. Other simulation parameters are (f₁, f₂, f₃, g₁, g₃, k _f , k _g ) = (1.3, 1.8, -1, 1, 1, 0.5, 0.5). In both the upper and lower panels, the first column compares the time evolution of S₁ molecules by different methods: the black line shows the ODE solution of Equation (52) for x₁ ; the blue lines are the solutions of Equation (53) for μ₁ and for μ₁ ± σ₁, respectively. The red dotted lines framing the mean indicate the 95% confidence interval. The second column shows the propensity adjustment functions for the forward reaction (solid line) and the backward reaction (dashed line). The third column shows 100 independent stochastic simulations with propensity adjustment (blue means and error bars), in comparison with the ODE (Equation (52)) prediction (black line). The fourth column shows a second set of 100 independent stochastic simulations without propensity adjustment (blue means and error bars), in comparison with the ODE (Equation (52)) prediction (black line). The red dotted lines framing the mean in columns 3 and 4 again indicate the 95% confidence intervals.

Repressilator

Interestingly, a propensity function may even be obtained through power-law approximation of some function that describes complex transient behaviours of a reaction network. As an example, consider the so-called Repressilator [27], which is a three-component genetic circuit where each component represses its downstream neighbour. More specifically (as shown in Figure 3), gene G₁ codes for protein x₁, whose dimer y₁ subsequently represses the transcription of the gene G₂. Similarly, y₂, the dimer of gene G₂'s protein product x₂, represses the transcription of gene G₃, and y₃, the dimer of gene G₃'s protein product x₃, represses the transcription of gene G₁. The corresponding differential equation model following mass action kinetics is given by [28]

Figure 3

Reaction scheme of the Repressilator. Gene G₁ codes for protein x₁, whose dimer y₁ represses the transcription of gene G₂. Similarly, y₂, the dimer of gene G₂'s protein product x₂, represses the transcription of gene G₃, and y₃, the dimer of gene G₃'s protein product x₃, represses the transcription of gene G₁.

$$\begin{array}{cc}\hfill {x_i}^{\prime }& =-2{\kappa }_{+}{x}_i^2+2{\kappa }_{-}{y}_i+\sigma m_i-{\gamma }_p{x}_i\hfill \\ \hfill {y_i}^{\prime }& ={\kappa }_{+}{x}_i^2-{\kappa }_{-}{y}_i-k_{+}{y}_i{d}_{0,j}+k_{-}{d}_{r,j}\hfill \\ \hfill {d_{0,i}}^{\prime }& =-k_{+}{y}_k{d}_{0,i}+k_{-}{d}_{r,i}\hfill \\ \hfill {d_{r,i}}^{\prime }& =k_{+}{y}_k{d}_{0,i}-k_{-}{d}_{r,i}\hfill \\ \hfill {m_i}^{\prime }& =d_{0,i}-{\gamma }_m{m}_i,\hfill \end{array}$$

(54)

where i = 1, 2, 3; j = 2, 3, 1; k = 3, 1, 2; the rate constants are explained in the diagram below

$$\begin{array}{c}{x}_i+x_i\underset{{\kappa }_{-}}{\overset{{\kappa }_{+}}{\rightleftharpoons }}{y}_i\\ d_{0,i}+y_k\underset{{\kappa }_{-}}{\overset{{\kappa }_{+}}{\rightleftharpoons }}{d}_{r,i}\\ d_{0,i}\underset{}{\overset{\alpha }{\to }}{d}_{0,i}+m_i\\ m_i\underset{}{\overset{\alpha }{\to }}{m}_i+x_i\\ x_i\underset{}{\overset{{\gamma }_p}{\to }}\varphi \\ m_i\underset{}{\overset{{\gamma }_m}{\to }}\varphi \end{array}$$

Assuming that the reversible dimerization and the dissociation/association of a protein dimer from/to the promoter are much faster than other processes, the full systems can be reduced to

$$\begin{array}{cc}\hfill {x_i}^{\prime }& =\sigma p{\left(x_i\right)}^{-1}{m}_i-{\gamma }_{p}p{\left(x_i\right)}^{-1}{x}_i\hfill \\ \hfill {m_i}^{\prime }& =\frac{\alpha d}{1+c_d{c}_p{x}_k^2}-{\gamma }_m{m}_i\hfill \end{array}$$

(55)

[28]. Here Φ = 1, $p\left(x_i\right)=1+4c_p{x}_i+\frac{4c_d{c}_{p}d{x}_i}{{\left(1+c_d{c}_p{x}_i^2\right)}^2}$, c _p = κ₊/κ_-, c _d = k₊/k_- and d = d_{0, i}+ d _{r, i} for i = 1, 2, 3. It has been shown that the simplified ODEs rather accurately approximate the transient dynamics of the full system by retaining the original oscillation period and amplitude.

In [28], the system (55) is further rescaled by setting $\tilde{t}={\gamma }_{m}t,{\tilde{x}}_i=\sqrt{c_d{c}_p}{x}_i$ and $\widetilde{m_i}=\left(\sigma \sqrt{c_d{c}_p}{m}_i\right)/\left({\gamma }_m\beta \right)$, which yields

$$\begin{array}{cc}\hfill \frac{d{\tilde{x}}_i}{d\tilde{t}}& =\beta p{\left({\tilde{x}}_i\right)}^{-1}\widetilde{m_i}-\beta p{\left({\tilde{x}}_i\right)}^{-1}{\tilde{x}}_i\hfill \\ \hfill \frac{d\widetilde{m_i}}{d\tilde{t}}& =\frac{\kappa d^{\prime }}{1+{\tilde{x}}_k^2}-\widetilde{m_i}.\hfill \end{array}$$

(56)

Intriguingly, one makes the following observation. The scaled ODE system (56) is consistent with the original system (55) in oscillation amplitude and period. However, its corresponding stochastic model produces results that deviate substantially from the average responses. To see the effects of the transition from a deterministic to a stochastic model, we apply SSA to the scaled system (56). The main result is that the oscillation periods of both x _i and m _i are reduced to half (Figure 4). The reason is that, in the stochastic simulation, the oscillation period is very sensitive to the ratio of x _i and m _i , which has been altered by the scaling operation. Therefore, in general one needs to pay attention to how scaling may affect the stochastic performance when the model is generated through the conversion of an ODE model.

Figure 4

Scaling of the Repressilator equations changes the oscillation period in the stochastic simulation. Solid lines represent solutions of ODEs (56), while dotted lines are trajectories of a stochastic simulation; blue lines represent x₁ and black lines represent m₁.

We can see from equations (55) that two variables x _i and m _i contribute to the production of x _i ; hence, their covariance could affect the propensity function of x _i in the production reaction of a stochastic simulation. Similar to the example of a reversible reaction (Equation 52), it is therefore necessary to evaluate covariance effects and to judge whether the propensity function needs adjusting. Thus, we need to compare the difference between the dynamics of the phenomenological model (55) and the dynamics under the influence of covariance, which can be produced by either stochastic simulation or the moment approach.

The influence of the covariance on the dynamics of the stochastic simulation is relatively easy to assess: we simply use the terms on the right-hand side of the differential equations (54) as the propensity functions in SSA and obtain simulation results shown in the 2^nd and the 4^th panels of Figure 5. Obtaining the covariance-influenced dynamics with the moment-based approach is more complicated, and we need to discuss some implementation issues.

Figure 5

Power-law approximation of p ( x _i )^-1. Left panel: Approximation of y _i = p(x _i )^-1 by a straight line in log-log space. Right panel: Corresponding power-law function in Cartesian space. Both axes are unitless.

First, the moment-based approach requires information regarding the first and the second derivatives of p(x _i )^-1, which have rather complicated functional forms. To simplify the calculation, we replace the function p(x _i )^-1 with an approximating power-law function. Specifically, suppose the original parameter values are κ₊ = k₊ = 5, κ_- = k_- = 100 and d = 20. Plotting the data (x _i , p(x _i )^-1)in log-log space (Figure 5) indicates that the original function is represented well by a straight line:

$$log{y}_i=log3.5188-0.9384log{x}_i.$$

for x _i ∈ [30, 300]. In Cartesian space, this line corresponds to the power-law function

$$y_i=3.5188x_i^{-0.9384},$$

which models the original function very well (see Figure 5). For x _i ∈ [1, 30], this power-law function does not fit the original function precisely; the effect of this imprecision can be evaluated later at after we use this power-law function in the moment-based method. Moreover, using the truncated moment equations to estimate the mean and variance involves multiple approximations: First, the function p (x _i )^-1 on the right-hand side of (55) is replaced by a power-law function (see Figure 5). Second, the result is approximated by Taylor expansion to the second order. Third, similar to the example of a reversible reaction, the central moment of the third degree is assumed to be zero, which leads to a closed-form ODE for the first two moments.

Solving the technical issues as described, one obtains the corresponding moment-based model of (55) (not shown) with results shown in Figure 6. Suppose one is particularly interested in the period and amplitude of the oscillation within a time interval between 0 and 400 seconds. As shown in Figure 6, the GMA approximation (black dashed line) fits the original ODEs (55) (bold black solid line) very well at the beginning, but as time goes on, the approximation error accumulates. As seen in the time interval [350, 400], the GMA approximation deviates from the original ODEs significantly. However, this does not mean that the GMA approximation cannot be used as a propensity function for stochastic simulations; the moment-based method with the GMA approximation shows that, when the GMA approximation is used as propensity function (without adjustment) for stochastic simulations, the resulting mean (red solid line) consistently fits the trajectory of the original ODEs (bold solid black line) very well up to about t = 400 seconds. The oscillation period and amplitude in the stochastic simulation based on the GMA approximation (without adjustment) are almost identical to those of the original ODEs. Therefore, a propensity adjustment for the GMA approximation is not needed, and the GMA approximation can be used as a propensity function for stochastic simulations. In other words, a stochastic model for the Repressilator system can be generated by using the scheme in (32) without propensity adjustment. Moreover, the imprecision caused by the power-law approximation can be tolerated when its corresponding moment-based mean matches the original ODEs sufficiently well with respect to the features of highest interest.

Figure 6

Comparison of the dynamics of the Repressilator models using the original ODEs (55), the GMA approximation, and the moment approach based on the GMA approximation. The mean of the moment approach based on the GMA approximation fits the original ODEs (55) very well up to about t = 400 s. Black bold line: solution of the original ODEs (55); black dashed line: the GMA approximation; Red line: mean of the moment approach based on the GMA approximation; red dashed lines framing around the red line: mean ± standard deviation, which were produced with the moment approach. x-axis is time in second, y-axis is the number of x₁ molecules (unitless).

Enzymatic reaction using a quasi-steady state assumption (QSSA)

We consider an enzymatic reaction following the Michaelis-Menten mechanism:

$$E+S\underset{k_{-1}}{\overset{k_1}{\rightleftharpoons }}ES\underset{}{\overset{k_2}{\to }}P+E.$$

Here enzyme E reacts with substrate S through a reversible reaction to form complex ES, which can proceed to yield product P and to release the enzyme E. By assuming the law of mass action for the reaction kinetics we obtain a set of differential equations for the system dynamics:

$$\begin{array}{cc}\hfill \frac{d\left[S\right]}{dt}& =k_{-1}\left[ES\right]-k_1\left[S\right]\left({\left[E\right]}_0-\left[ES\right]\right)\hfill \\ \hfill \frac{d\left[ES\right]}{dt}& =k_1\left[S\right]\left({\left[E\right]}_0-\left[ES\right]\right)-\left(k_1+k_2\right)\left[ES\right]\hfill \\ \hfill \frac{d\left[P\right]}{dt}& =k_2\left[ES\right],\hfill \end{array}$$

(57)

where the total amount of enzyme in the form of free enzyme and complex [E]₀ ≜ [E] + [ES] is assumed to be constant. In addiction, by making the so-called quasi steady state assumption (QSSA) [29, 30], assuming that the complex ES is essentially in steady state, we can assert $\frac{d\left[ES\right]}{dt}\approx 0$. As it has been discussed many times in the literature, QSSA reduces the system and leads to the approximate form

$$\begin{array}{cc}\hfill \frac{d\left[S\right]}{dt}& =-\frac{V_{max}\left[S\right]}{K_m+\left[S\right]}\hfill \\ \hfill \frac{d\left[P\right]}{dt}& =\frac{V_{max}\left[S\right]}{K_m+\left[S\right]},\hfill \end{array}$$

(58)

which is known as Michaelis-Menten kinetics [30]. The characterizing parameters are V_max = k₂[E]₀ and K _m = (k_-1 + k₂)/k₁.

Applying QSSA, Rao and Arkin [2] were able to reduce the CME of S and ES to a CME only containing S. For the reduced CME, the propensity function for the overall reaction S → P is

$$\alpha \left(s\right)=\frac{V_{max}s}{K_m+s},$$

(59)

where the volume was scaled so that Φ = 1 and the lower-case letter s denotes the molecule count of species S. Instead of reviewing the relatively complicated manipulations with CME, we show in the following that the techniques described above lead directly from the equation-based model to the propensity function for the reduced system.

First, we recast the equation-based model into the GMA format [31], by introducing an auxiliary variable [T] ≜ K _m + [S]. The result,

$$\begin{array}{cc}\hfill \frac{d\left[S\right]}{dt}& =-V_{max}\left[S\right]{\left[T\right]}^{-1}\hfill \\ \hfill \frac{d\left[T\right]}{dt}& =\frac{d\left[S\right]}{dt}=-V_{max}\left[S\right]{\left[T\right]}^{-1}\hfill \end{array}$$

(60)

is exactly equivalent to the reduced system in (58) with the initial condition [S]₀ and [T]₀ = K _m + [S]₀. The corresponding stochastic model has only one reaction channel and the propensity function is

$$\alpha \left(s,t\right)=V_{max}s{t}^{-1}.$$

(61)

The propensity adjustment factor can be set to 1 because T is a function of s and its covariance with s is therefore 1. By applying t = K _m + s, the propensity function can be simplified as

$$\alpha \left(s,t\right)=V_{max}s{t}^{-1}=V_{max}s{\left(K_m+s\right)}^{-1}=\alpha \left(s\right).$$

(62)

Thus, we arrive at the propensity function for the reduced system, which is identical to the result of Rao and Arkin obtained through manipulations of CME.

In the above derivation, we used the simplest type of recasting, where a new, auxiliary variable simply consists of an old variable plus a constant. This reformulation of the Michaelis-Menten process as a pair of GMA equations is a special case of a much more general recasting technique that permits the equivalent conversion of any system of ordinary differential equations into a power-law format [31]. However, this equivalence transformation imposes constraints on the variables of the GMA equations, and it is at this point unclear whether there are mathematical warranties ensuring that the proposed transition from differential to stochastic equations in general preserves these constraints in all cases. This question will require further investigation.

Stochastic Focusing

Stochastic focusing [26] describes the phenomenon that the fluctuations of a chemical species can drive the system to reach a different steady state than what a deterministic ODE model predicts. To demonstrate the utility of propensity adjustment, we derive a stochastic model which produces consistent results with those of the deterministic model.

Following [32], we consider the following reactions system

$$\begin{array}{c}\varphi \underset{}{\overset{k_1}{\to }}I\underset{}{\overset{k_2}{\to }}P\underset{}{\overset{k_3}{\to }}\varphi \\ I+S\underset{}{\overset{k_4}{\to }}S\\ S\underset{k_5}{\overset{k_6}{\rightleftharpoons }}\varphi \end{array}$$

(63)

This system can be interpreted as follows: the intermediate species I is produced at constant rate k₁ from some source Φ and degrades with rate k₄ through the catalysis with signalling molecule S; the end product P is converted from species I at rate k₂ and degrades at rate k₃; the signalling molecule S is produced and degrades at rates k₅ and k₆, respectively. Moreover, the value of k₅ is reduced to half at a certain time point to achieve a significant divergence effect. In order to capture the average dynamics of the system accurately, we use a power-law model in GMA format instead of the mass action rate law in [32].

$$\begin{array}{cc}\hfill \frac{di}{dt}& =k_1-k_{2}i-k_4{i}^{f_I}{s}^{f_S}\hfill \\ \hfill \frac{dp}{dt}& =k_{2}i-k_{3}p\hfill \\ \hfill \frac{ds}{dt}& =k_5-k_{6}s\hfill \end{array}$$

(64)

The system size is set to 1. We can see from equations (64) that two variables i and s contribute to the degradation of I and that their covariance could therefore affect the propensity function of I in the degradation reaction of a stochastic simulation. To calculate the propensity adjustment function paf₄(t) = exp (-f _I f _S cov [log I(t), log S(t)]) for reaction $R_4:I+S\underset{}{\overset{k_4}{\to }}S$, we formulate equations (cf. (60)) for the moments as

$$\begin{array}{cc}\hfill \frac{\partial \mu }{\partial t}& \approx V^T\left(\alpha +\frac{1}{2}{\alpha }^{″}\odot \sigma \right)\hfill \\ \hfill \frac{\partial \sigma }{\partial t}& \approx {\left(\sigma \left({\alpha }^{\prime }V\right)\right)}^T+\sigma \left({\alpha }^{\prime }V\right)+V^T\mathrm{\Lambda }V.\hfill \end{array}$$

(65)

Here μ = (μ _I , μ _P , μ _S )^T, $V=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill \end{array}\right]$, $\alpha =\left[\begin{array}{c}\hfill {\alpha }_1\hfill \\ \hfill {\alpha }_2\hfill \\ \hfill {\alpha }_3\hfill \\ \hfill {\alpha }_4\hfill \\ \hfill {\alpha }_5\hfill \\ \hfill {\alpha }_6\hfill \end{array}\right]=\left[\begin{array}{c}\hfill k_1\hfill \\ \hfill k_{2}i\hfill \\ \hfill k_{3}p\hfill \\ \hfill k_4{i}^{f_I}{s}^{f_S}\hfill \\ \hfill k_5\hfill \\ \hfill k_{6}s\hfill \end{array}\right]$.

Moreover, for r = 1, ..., 6 and m, n = i, p, s, ${\left({\alpha }_r^{{}^{″}}\right)}_{mn}=\frac{{\partial }^2{\alpha }_r}{\partial m\partial n}$, α" = (α₁", ..., α₆")^T, $\sigma =\left[\begin{array}{ccc}\hfill {\sigma }_{11}\hfill & \hfill {\sigma }_{12}\hfill & \hfill {\sigma }_{13}\hfill \\ \hfill {\sigma }_{21}\hfill & \hfill {\sigma }_{22}\hfill & \hfill {\sigma }_{23}\hfill \\ \hfill {\sigma }_{31}\hfill & \hfill {\sigma }_{32}\hfill & \hfill {\sigma }_{33}\hfill \end{array}\right]$, ${\alpha }_r^{{}^{″}}\odot \sigma =\sum _{m,n=i,p,s}\frac{{\partial }^2{\alpha }_r\left(\mathbf{\mu }\right)}{\partial m\partial n}{\sigma }_{mn}$, α"⊙ σ ≜ (α₁ "⊙ σ, ..., α₆"⊙ σ)^T, α' = (α₁', ..., α₆'), ${{\alpha }_r}^{\prime }={\left(\frac{\partial {\alpha }_r}{\partial i},\frac{\partial {\alpha }_r}{\partial p},\frac{\partial {\alpha }_r}{\partial s}\right)}^T$, and the diagonal matrix Λ is defined by ${\left(\mathrm{\Lambda }\right)}_{rr}={\alpha }_r+\frac{1}{2}\sum _{m,n=i,p,s}\frac{{\partial }^2{\alpha }_r}{\partial m\partial n}{\sigma }_{mn}$.

The stochastic focusing model without propensity adjustment yields results quite different from those of the deterministic model, as is illustrated in Figure 7. In this figure, the blue lines in the 1^st panel are predicted from the moment equations (65) and the blue error bars for μ _P in the 2^nd panel are obtained from ten independent stochastic simulations. Both diverge systematically from the black line predicted by ODE model (64). By contrast, the stochastic model with propensity adjustment produces results consistent with the deterministic model, as shown by the 4^th panel.

Figure 7

Stochastic focusing. The first panel from the top compares the time evolution of product molecules P obtained with different methods: the black line represents the solution of ODE model (64) for P; the blue solid line and blue dashed lines are the solutions of the moment-based model (65) for μ _P and μ _P ± σ _P , respectively. The second panel indicates that the stochastic simulations without propensity adjustment (blue error bar) diverge from the prediction by the ODE model (64) (black line). The third panel shows the propensity adjustment function paf4(t) = exp(-f _I f _S cov [log I(t), log S(t)]) for the reaction $R_4:I+S\underset{}{\overset{k_4}{\to }}S$. The bottom panel demonstrates that the propensity adjustment function paf₄ achieves convergence between the stochastic simulation and the ODE model (64) (black line): the blue error bars were computed from 100 independent stochastic simulations with propensity adjustment paf₄. The simulation parameters are (i(0), p(0), s(0), k₁, k₂, k₃, k₄, k₅, k₆, f _I , f _S ) = (0, 10, 100, 10⁴, 10³, 1, 9.9 × 10³, 10⁴, 10³, 1.1, 0.9); at t = 0.1, the value of k₅ changes from 10⁴ to 0.5 ×10⁴.

Discussion

It is often implicitly assumed that the rate of a dynamic process can be directly taken as the propensity for a corresponding stochastic process. We have shown here that this is sometimes, but not always, true. Our results fall into three categories. The first develops conditions for a valid conversion of a rate to a propensity, the second presents a general conversion procedure, and the third discusses computational issues of propensity adjustment.

Conditions for the direct use of a rate constant (function) as propensity function

We have shown that the direct use of a rate constant or a rate function f as the propensity function in a stochastic simulation algorithm requires that at least one of the following assumptions be true:

• 1) f is a linear function; this assumption has been validated in the Results sections addressing 0^th-order and 1^st-order reaction kinetics.

• 2) the reaction is monomolecular; this assumptions was evaluated in the Results section describing real-valued order monomolecular reaction kinetics.

• 3) all X _i in the system are noise-free variables, i.e., without (or with ignorable) fluctuations; this assumption implies that the covariance of any two participating reactants is zero (or close to zero). This assumption is assessed in equations (29 - 36).

Each of these three conditions is a sufficient condition for the direct use of a rate function f as the propensity function. Moreover, these statements are valid for functions of a general format, not just for GMA. This is so because the functional formats in cases 1 and 2 above are special cases of the GMA format. For the third case, a formal proof is only given for functions in GMA format, because this structured format allows us to give an explicit estimation on how the covariance can affect the average behavior of a stochastic simulation through equation (34). For functions not in GMA format, the conclusion is still holds, although an analogous explicit estimation is lacking. The argument is as follows. The bimolecular reaction E[α _r (X(t))] contains at least one quadratic moment of the form E[X _i (t)X _j (t)] (cf. [4] and page 38). Therefore, by definition of the covariance, E[X _i (t)X _j (t)] = E[X _i (t)]E[X _i (t)] + cov (X _i (t), X _j (t)), we obtain

$$E\left[X_i\left(t\right)X_j(\right]$$