A Comparison of Deterministic and Stochastic Modeling Approaches for Biochemical Reaction Systems: On Fixed Points, Means, and Modes

Abstract

In the mathematical modeling of biochemical reactions, a convenient standard approach is to use ordinary differential equations (ODEs) that follow the law of mass action. However, this deterministic ansatz is based on simplifications; in particular, it neglects noise, which is inherent to biological processes. In contrast, the stochasticity of reactions is captured in detail by the discrete chemical master equation (CME). Therefore, the CME is frequently applied to mesoscopic systems, where copy numbers of involved components are small and random fluctuations are thus significant. Here, we compare those two common modeling approaches, aiming at identifying parallels and discrepancies between deterministic variables and possible stochastic counterparts like the mean or modes of the state space probability distribution. To that end, a mathematically flexible reaction scheme of autoregulatory gene expression is translated into the corresponding ODE and CME formulations. We show that in the thermodynamic limit, deterministic stable fixed points usually correspond well to the modes in the stationary probability distribution. However, this connection might be disrupted in small systems. The discrepancies are characterized and systematically traced back to the magnitude of the stoichiometric coefficients and to the presence of nonlinear reactions. These factors are found to synergistically promote large and highly asymmetric fluctuations. As a consequence, bistable but unimodal, and monostable but bimodal systems can emerge. This clearly challenges the role of ODE modeling in the description of cellular signaling and regulation, where some of the involved components usually occur in low copy numbers. Nevertheless, systems whose bimodality originates from deterministic bistability are found to sustain a more robust separation of the two states compared to bimodal, but monostable systems. In regulatory circuits that require precise coordination, ODE modeling is thus still expected to provide relevant indications on the underlying dynamics.

Keywords: bimodality; bistability; chemical master equation; gene expression; ordinary differential equations; protein bursts.

Figure 1

Influence of bursting and of nonlinear feedback on the protein distribution. From left to right, the feedback characteristics are varied. The top row shows the analytical results. The deterministic fixed points can be read from the value of n at the intersection of $\frac{f(.)}{\delta }$, marked in red, and the identity line in green. The systems are monostable except for column (C), where it is bistable. The intersection points of the red line and the blue lines .+μ^* yield the locations of the extrema. Three different values for μ^* are shown: μ^*₁ = 1 (no burst, dark blue), μ^*₂ = 6 (medium-size burst, mid blue), μ^*₃ = 11 (strong burst, light blue). Through bursts, the modes are shifted toward smaller numbers of protein molecules. The second row shows the histograms of the protein distribution obtained from 5·10⁴ protein time-course simulations using the Gillespie algorithm. The distribution is shown for each burst size (same color as above). The location of the extrema corresponds well to the analytical results. The average values (dashed lines) match the deterministic steady state if f is linear, which is only the case in panel (A). In (C), large bursts generate a unimodal distribution (marked in light blue), although the system is bistable. In (B,D), medium-size bursts lead to bimodality in spite of deterministic monostability (mid-blue line). Parameters are given in the Supplementary Table 1.

Figure 2

Influence of the system size on the correspondence between deterministic and stochastic modeling results. Two systems with differing sizes are compared: The volume V₁ of system 1 (graphs in light blue) is chosen 50-fold smaller than the volume V₂ of system 2 (graphs in dark blue), while the protein concentrations at the deterministic fixed points are identical. The intersections of the blue lines and the red line in the upper plot mark the analytical locations of the extrema in the protein probability mass function. The extrema of the larger system nearly coincide with the deterministic fixed points, since the expression $\frac{\mu *}{V_2}$ is almost negligible. The distributions in the bottom plot (obtained using the Gillespie algorithm) confirm these results: the larger system shows a clear bimodal distribution whose modes match the stable deterministic fixed points, while the modes of the small system are shifted, and the distribution is much broader. The dashed lines show that the analytical determination of the modes fits well to the simulations. Parameters are given in the Supplementary Table 2.

Figure 3

Robustness of bimodality in different regulatory systems with feedback. In each column, the robustness of the modes in two regulatory systems with varying burst sizes and varying functions f are compared based on the protein distributions and exemplary protein time-courses. (A) Comparison of two non-cooperative feedback regulations with differing burst sizes μ^*₁ < μ^*₂ and accordingly shifted functions f₁ and f₂ with identical shape. (B) Comparison of non-cooperative and cooperative regulation with identical burst size μ^*. (C) Comparison of two non-cooperative regulations with identical basal protein expression and under differing burst sizes μ^*₆ < μ^*₅. In all cases, the system marked in dark colors (system 1, 3, and 5, respectively) exhibits a sharper distribution and a better separation of the modes is visible in the protein time-course simulations. Further explanations are given in the main text. Parameter values are listed in the Supplementary Table 3.