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)
0th-order reaction kinetics
Consider a very simple equation-based model of the type
(47)
for all s = 1, ..., N
s
, f
rs
= 0. According to Equations (40)-(44), one obtains
Thus, for a 0th-order reaction, its rate equation can be taken directly as the propensity function in stochastic simulations.
-
2)
1st-order reaction kinetics
Direct application of Equations (40)-(44) yields
(48)
f
rs
= δ
sj
, i, j = 1, ..., N
s
. Therefore, according to Equations (40)-(44)
Thus, for 1st-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
(49)
f
rj
≠ 0, f
rs
= 0, for any s ≠ j, s = 1, ..., N
s
. Equations (40)-(44) lead to
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)
2nd-order reaction kinetics
This type of reaction can be expressed as
(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)
Thus, the proper propensity function for 2nd-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
(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
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
(52)
Here S3 feeds back to inhibit the forward reaction and S1 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 x1, x2 and x3 contribute to the forward flux and two variables x1 and x3 contribute to the backward flux . 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,
(53)
Here μ = (μ1, μ2, μ3)T, , .
Moreover, for r = 1, 2 and m, n = 1, 2, 3, , α" = (α1", α2")T, , , α"⊙σ ≜ (α1"⊙ σ, α2"⊙ σ)T, α' = (α1', α2'), , and .
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.
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 G1 codes for protein x1, whose dimer y1 subsequently represses the transcription of the gene G2. Similarly, y2, the dimer of gene G2's protein product x2, represses the transcription of gene G3, and y3, the dimer of gene G3's protein product x3, represses the transcription of gene G1. The corresponding differential equation model following mass action kinetics is given by [28]
(54)
where i = 1, 2, 3; j = 2, 3, 1; k = 3, 1, 2; the rate constants are explained in the diagram below
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
(55)
[28]. Here Φ = 1, , c
p
= κ+/κ-, c
d
= k+/k- and d = d0, 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 and , which yields
(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.
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 2nd and the 4th 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.
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:
for x
i
∈ [30, 300]. In Cartesian space, this line corresponds to the power-law function
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.
Enzymatic reaction using a quasi-steady state assumption (QSSA)
We consider an enzymatic reaction following the Michaelis-Menten mechanism:
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:
(57)
where the total amount of enzyme in the form of free enzyme and complex [E]0 ≜ [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 . As it has been discussed many times in the literature, QSSA reduces the system and leads to the approximate form
(58)
which is known as Michaelis-Menten kinetics [30]. The characterizing parameters are Vmax = k2[E]0 and K
m
= (k-1 + k2)/k1.
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
(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,
(60)
is exactly equivalent to the reduced system in (58) with the initial condition [S]0 and [T]0 = K
m
+ [S]0. The corresponding stochastic model has only one reaction channel and the propensity function is
(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
(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
(63)
This system can be interpreted as follows: the intermediate species I is produced at constant rate k1 from some source Φ and degrades with rate k4 through the catalysis with signalling molecule S; the end product P is converted from species I at rate k2 and degrades at rate k3; the signalling molecule S is produced and degrades at rates k5 and k6, respectively. Moreover, the value of k5 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].
(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 paf4(t) = exp (-f
I
f
S
cov [log I(t), log S(t)]) for reaction , we formulate equations (cf. (60)) for the moments as
(65)
Here μ = (μ
I
, μ
P
, μ
S
)T, , .
Moreover, for r = 1, ..., 6 and m, n = i, p, s, , α" = (α1", ..., α6")T, , , α"⊙ σ ≜ (α1 "⊙ σ, ..., α6"⊙ σ)T, α' = (α1', ..., α6'), , and the diagonal matrix Λ is defined by .
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 1st panel are predicted from the moment equations (65) and the blue error bars for μ
P
in the 2nd 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 4th panel.