Classical versus stochastic kinetics modeling of biochemical reaction systems

Abstract

We study fundamental relationships between classical and stochastic chemical kinetics for general biochemical systems with elementary reactions. Analytical and numerical investigations show that intrinsic fluctuations may qualitatively and quantitatively affect both transient and stationary system behavior. Thus, we provide a theoretical understanding of the role that intrinsic fluctuations may play in inducing biochemical function. The mean concentration dynamics are governed by differential equations that are similar to the ones of classical chemical kinetics, expressed in terms of the stoichiometry matrix and time-dependent fluxes. However, each flux is decomposed into a macroscopic term, which accounts for the effect of mean reactant concentrations on the rate of product synthesis, and a mesoscopic term, which accounts for the effect of statistical correlations among interacting reactions. We demonstrate that the ability of a model to account for phenomena induced by intrinsic fluctuations may be seriously compromised if we do not include the mesoscopic fluxes. Unfortunately, computation of fluxes and mean concentration dynamics requires intensive Monte Carlo simulation. To circumvent the computational expense, we employ a moment closure scheme, which leads to differential equations that can be solved by standard numerical techniques to obtain more accurate approximations of fluxes and mean concentration dynamics than the ones obtained with the classical approach.

FIGURE 1

Normalized dimer accumulation in the unidirectional dimerization reaction, given by Eq. 11, predicted by the SCKEs (solid lines) and CKEs (dotted lines). The dynamics obtained by the SCKEs have been computed by Monte Carlo simulation using the Gillespie algorithm, whereas the dynamics obtained by the CKEs have been computed analytically from Eq. 38. The system is initialized with (A) one molecule P and one molecule Q, (B) 10 molecules P, and 10 molecules Q. The associated normalized flux and mesoscopic forcing term dynamics are depicted as well. Parameters used are c₁ = 10⁻³ s⁻¹, V = 2 pL, and K = 6000 cells.

FIGURE 2

Protein accumulation in the quadratic autocatalator, given by Eq. 12, for the case when κ₅ = κ₂d, initialized with 10 molecules S (concentration of 8.30 pM), two molecules D (the number of DNA copies of a particular gene per eukaryotic cell), and zero molecules P and Q, predicted by the SCKEs (solid lines) and CKEs (dotted lines). The dynamics obtained by the SCKEs have been computed by Monte Carlo simulation using the Gillespie algorithm, whereas the dynamics obtained by the CKEs have been computed numerically. The flux dynamics of the third and fourth reactions are depicted as well. Parameters used are c₁ = 0.002 s⁻¹, c₂ = 0.001 s⁻¹, c₃ = 0.005 s⁻¹, c₄ = 0.004 s⁻¹, c₅ = 0.002 s⁻¹, c₆ = 0.05 s⁻¹, V = 2 pL, and K = 10,000 cells. Although the steady-state concentration of Q predicted by the CKEs is theoretically identical to the one predicted by the SCKEs, this is not true for the concentration of P. The dashed lines indicate the mean concentration and flux dynamics predicted by the second-order SCKEs discussed in this article.

FIGURE 3

Protein accumulation and flux dynamics in the quadratic autocatalator, given by Eq. 12, for the case when κ₅ > κ₂d. The parameters used are the same as in Fig. 2, but now c₅ = 0.006 s⁻¹. In this case, the steady-state concentrations of P and Q predicted by the CKEs (dotted lines) are both different than the actual steady-state concentrations predicted by the SCKEs (solid lines). The dashed lines indicate the mean concentration and flux dynamics predicted by the second-order SCKEs discussed in this article.

FIGURE 4

The input flux k₁s versus the stationary mesoscopic forcing term , the stationary concentration of P as a function of s, predicted by the SCKEs (solid lines) and CKEs (dotted lines), and the ratio associated with the quadratic autocatalator, given by Eq. 12. The steady-state values obtained by the SCKEs have been computed by Monte Carlo simulation using the Gillespie algorithm, whereas the values obtained by the CKEs have been computed analytically. The system is initialized with two molecules D (the number of DNA copies of a particular gene per eukaryotic cell), and zero molecules P and Q. Parameters used are (A) c₁ = 0.002 s⁻¹, c₂ = 0.0005 s⁻¹, c₃ = 0.005 s⁻¹, c₄ = 0.004 s⁻¹, c₅ = 0.004 s⁻¹, c₆ = 0.05 s⁻¹, and (B) c₁ = 0.0004 s⁻¹, c₂ = 0.02 s⁻¹, c₃ = 0.05 s⁻¹, c₄ = 0.04 s⁻¹, c₅ = 0.01 s⁻¹, and c₆ = 0.05 s⁻¹. Moreover, V = 2 pL and K = 8000 cells. The heavy bold line in the middle figure depicts the steady-state response curve of P, calculated by Monte Carlo simulation using the Gillespie algorithm.

FIGURE 5

CV dynamics in the quadratic autocatalator, given by Eq. 12, associated with intrinsic stochastic fluctuations in the concentrations of P and Q, for the case when κ₅ = κ₂d, in panel A, and κ₅ > κ₂d, in panel B, predicted by the exact SCKEs (solid lines) and second-order SCKEs (dashed lines). The dynamics obtained by the exact SCKEs have been computed by Monte Carlo simulation using the Gillespie algorithm, whereas the dynamics obtained by the second-order SCKEs have been computed numerically. The parameters used are the same as in Figs. 2 and 3.

FIGURE 6

Absolute relative errors in the steady-state concentrations of P and Q associated with the quadratic autocatalator, given by Eq. 12, as a function of the input substrate concentration. (Solid lines) Second-order SCKEs with respect to the exact SCKEs; (dotted lines) CKEs with respect to the exact SCKEs; and (dashed lines) CKEs with respect to the second-order SCKEs. The steady-state values obtained by the exact and second-order SCKEs have been respectively computed by Monte Carlo simulation using the Gillespie algorithm and numerically, whereas the values obtained by the CKEs have been computed analytically from Eq. 44. The system is initialized with two molecules D (the number of DNA copies of a particular gene per eukaryotic cell), and zero molecules P and Q. Parameters used are c₁ = 0.002 s⁻¹, c₂ = 0.001 s⁻¹, c₃ = 0.005 s⁻¹, c₄ = 0.004 s⁻¹, c₆ = 0.05 s⁻¹, V = 2 pL, and K = 8000 cells. Moreover, c₅ = 0.002 s⁻¹ in panel A, and c₅ = 0.006 s⁻¹ in panel B.