Deriving the mean and variance of a powerlaw function of random variables
Consider a generic powerlaw 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}\phantom{\rule{0.3em}{0ex}}\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}\phantom{\rule{0.3em}{0ex}}\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
}^{2} 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(\right)close="\}">\frac{1}{4}{\left({\sigma}_{i}/{\mu}_{i}\right)}^{2}{\left({\sigma}_{j}/{\mu}_{j}\right)}^{2}& .\end{array}\n
(13)
Since many biochemical variables approximately follow a lognormal distribution [20–22], it is valuable to consider the special situation where (X_{1}, ..., X_{
s
})is lognormally distributed (i.e., (logX_{1}, ..., 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
}^{2} 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
}^{2} are related to μ_{
s
}, σ_{
s
}^{2} 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
}^{2} on (μ_{
s
}, σ_{
s
}^{2}, σ_{
ij
}), but it is only valid if the additional assumption of lognormality is acceptable.
Deriving proper propensity functions for stochastic simulations from differential equationbased 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}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{probability}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{that}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{exactly}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{one}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{reaction}}\\ {R}_{r}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{will}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{occur}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{some}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{where}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{inside}}\phantom{\rule{2.77695pt}{0ex}}U\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{within}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{infinitesmal}}\\ \mathsf{\text{interval}}\phantom{\rule{2.77695pt}{0ex}}\left(t,t+dt\right),\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{given}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{current}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{state}}\phantom{\rule{2.77695pt}{0ex}}\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},\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{given}}\phantom{\rule{2.77695pt}{0ex}}\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 \\ \phantom{\rule{1em}{0ex}}{\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 wellstirred 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),\phantom{\rule{1em}{0ex}}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}}\phantom{\rule{2.77695pt}{0ex}}{x}_{s}\ge {\underset{}{v}}_{rs}>0\hfill \\ \hfill 0,\phantom{\rule{2.77695pt}{0ex}}\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 positivevalued differentiable function can be approximated locally by a powerlaw 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 realvalued. 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],\phantom{\rule{2.77695pt}{0ex}}s=1,\dots ,\phantom{\rule{2.77695pt}{0ex}}{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}}},\phantom{\rule{1em}{0ex}}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 "particlebased" 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}}},\phantom{\rule{1em}{0ex}}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 lefthand side of equation (29):

1)
adopt a zerocovariance 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 0covariance propensity function from a second type of propensity in the next section.
With the zerocovariance 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 equationbased 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 equationbased 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 momentbased 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}\phantom{\rule{0.3em}{0ex}}\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}\phantom{\rule{0.3em}{0ex}}\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 valueorder 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 zerocovariance 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 equationbased 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 nonzero 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{\xaf}{v}}_{rs}}!}{\displaystyle \prod _{s=1}^{{N}_{s}}{X}_{s}^{{\underset{\xaf}{v}}_{rs}}{}^{+{\epsilon}_{rs}}}\right]\\ \approx \frac{{\kappa}_{r}}{{\displaystyle \prod _{s=1}^{{N}_{s}}{\underset{\xaf}{v}}_{rs}!}}{\displaystyle \prod _{s=1}^{{N}_{s}}E{\left[{X}_{s}\right]}^{{\underset{\xaf}{v}}_{rs}+{\epsilon}_{rs}}}\\ \times \mathrm{exp}\left({\displaystyle \sum _{i<j}^{{N}_{s}}\left({\underset{\xaf}{v}}_{ri}+{\epsilon}_{ri}\right)\left({\underset{\xaf}{v}}_{rj}+{\epsilon}_{rj}\right)}\text{}\text{cov}\left[1\text{og}{X}_{i},\mathrm{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}\phantom{\rule{2.77695pt}{0ex}}\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\phantom{\rule{2.77695pt}{0ex}}\left(\sum _{i<j}^{{N}_{s}}{f}_{ri}{f}_{rj}\phantom{\rule{2.77695pt}{0ex}}\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∈S}appear 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 welljustified 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\phantom{\rule{2.77695pt}{0ex}}\left(\sum _{i<j}^{{N}_{s}}{f}_{ri}{f}_{rj}\phantom{\rule{2.77695pt}{0ex}}\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}\phantom{\rule{2.77695pt}{0ex}}\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 zerocovariance, 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}}\phantom{\rule{1em}{0ex}}}}=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}}\phantom{\rule{1em}{0ex}}}}^{2}\approx {{\mu}_{{\alpha}_{{r}_{}\mathsf{\text{cov}}\phantom{\rule{1em}{0ex}}}}}^{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}\phantom{\rule{2.77695pt}{0ex}}\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_{1}, ..., X_{
s
}) are lognormally distributed [8, 23]. The consideration of a lognormal distribution is often justified by the fact that many biochemical data have indeed been observed to be lognormally 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 momentbased 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 \\ \phantom{\rule{1em}{0ex}}\left(\right)close="\}">+{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]\hfill \end{array}\n
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}^{\u2033}\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}^{\u2033}={\left({\alpha}_{r}^{{}^{\u2033}},.\phantom{\rule{2.77695pt}{0ex}}.\phantom{\rule{2.77695pt}{0ex}}.\phantom{\rule{2.77695pt}{0ex}},{\alpha}_{{N}_{r}}^{{}^{\u2033}}\right)}^{T}, {\left({\alpha}_{r}^{{}^{\u2033}}\right)}_{mn}=\frac{{\partial}^{2}{\alpha}_{r}\left(\mathbf{X}\right)}{\partial {X}_{m}\partial {X}_{n}}, {\alpha}_{r}^{{}^{\u2033}}\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}^{{}^{\u2033}}\odot \sigma \triangleq {\left({\alpha}_{1}^{{}^{\u2033}}\odot \sigma ,\dots ,{\alpha}_{{N}_{r}}^{{}^{\u2033}}\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 equationbased 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 momentbase 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 momentbased method or from n runs of stochastic simulations. The Repressilator example in the result section will serve as a demonstration.