Determining a Reasonable Range of Relative Numerical Tolerance Values for Simulating Deterministic Models of Biochemical Reactions

Adeterministic model of a biological process in the form of differential equations presents mathematical relationships among variables of precise values, but in simulation of such a model, only approximate values of the variables can be employed. The approximation is governed by either an absolute or a relative numerical tolerance value. An absolute numerical tolerance value is the absolute difference between the actual and expected values. Dividing such an absolute difference by the expected value gives rise to a relative numerical tolerance value. Because a model of a biological system may produce different simulation outcomes with different ranges of numerical tolerance values, it is important to determine a reasonable range of numerical tolerance values. Currently a reasonable range of numerical tolerance values for modeling a biological process such as that involving biochemical reactions has not been explored.

Numerical tolerance issues have been studied before, which led to the general conclusion that decreasing numerical tolerance values to a sufficiently small and yet practical level may effectively correct erroneous outcomes in computer simulation (Byrne and Thompson, 2013). The erroneous outcomes arise from numerical stiffness, a phenomenon of unacceptably large errors in dependent variables as a result of small errors in independent variables (DiStefano III, 2015). Here we present a case that a stable numerical outcome may be achieved only when the numerical tolerance values approach zero. This case was previously described in an attempt to model oscillations of biochemical reactions in a negative feedback system (Yang and Yang, 2015). In this case, the magnitudes of the oscillations of chemicals x and y decrease along with the reduction in numerical tolerance values, but chemical concentrations remain oscillatory even when the relative and absolute numerical tolerance values are reduced to 10⁻¹⁰ and 10⁻¹², respectively (Fig. 1). It seems that the chemical concentrations will remain oscillatory unless the numerical tolerance values approach zero or are much smaller than those shown in Figure 1c. A question then arises: which simulation outcome should be accepted in this case, an outcome with a certain range of numerical tolerance values far from zero or one with the numerical tolerance values approaching or close to zero?

OPEN IN VIEWER

FIG. 1.

Simulation outcomes of a negative feedback system of biochemical reactions at three numerical tolerance levels. This is a two-ODE system numerically integrated as described for figure 2B in Yang and Yang (2015). (a) Relative numerical tolerance = 10⁻³ and absolute numerical tolerance = 10⁻⁶. (b) Relative numerical tolerance = 10⁻⁵ and absolute numerical tolerance = 10⁻⁸. (c) Relative numerical tolerance = 10⁻¹⁰ and absolute numerical tolerance = 10⁻¹².

We would argue that the answer to the mentioned question is an outcome with the numerical tolerance values far from zero. We think that at the fundamental level, chemical concentrations cannot be precisely determined and the inherent variance of a chemical concentration is governed by the uncertainly principle. The uncertainty principle was typically applied to a subatomic particle. Because each molecule can be considered as a population of subatomic particles, the uncertainty level of such a population in position and energy state should be equal to or more than that of a subatomic particle, depending on the complexity of the molecule. The uncertainty level of a subatomic particle, therefore, can serve as a guidance value for the estimation of the lower limit of the uncertainty level of molecules. Precisely calculating the uncertainty level of a molecule, or a population of molecules in the case of determining the concentration of a molecular species, is challenging due to the complexity arising from the large number of subatomic particles involved and the dynamics of the molecule dictated by its internal structure and external environment.

Toward finding the lower limit of the uncertainty level of a population of molecules, we start from the uncertainty level of a subatomic particle, which is typically presented as the following: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \sigma _x } { \sigma _p } \ge \frac { h } { { { \rm { 4 } } \pi } } , \tag { 1 } \end{align*} \end{document}

where σ_x and σ_p are the standard deviations of the position and momentum of a particle, respectively, and h is the Planck constant. Now let us consider a biochemical reaction as a representative case. The kinetics of a chemical reaction is typically described by an equation with variables of the concentrations of the reactant(s) and product(s), and the standard deviation of the concentration σ_c is, in effect, the relative numerical tolerance value, which exerts its influence on the simulation when the concentration values are reasonably large (e.g., >0.1). The standard deviation of the concentration of a chemical, σ_c, can be calculated as the following: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \sigma _c } = \sqrt { \frac { { \sum \limits_ { i = 1 } ^n \displaystyle { \bigg( { \frac { { w_i } } { { v_i } } } } - { \frac { { w_0 } } { { v_0 } } } { \bigg) ^2 } } } { n } } , \tag { 2 } \end{align*} \end{document}

where w₀, v₀, and n are the mean of the chemical weight measurement, mean of the volume measurement, and the number of the measurements, respectively.

To simplify the analysis, it is assumed that the weight measurements are quite accurate, that is, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${w_i} \approx {w_0}$$ \end{document} , therefore, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \sigma _c } \approx { w_0 } \sqrt { \frac { { \sum \limits_ { i = 1 } ^n \displaystyle { { { \bigg( { \frac { { v_i } - { v_0 } } { { v_i } { v_o } } } \bigg) } ^2 } } } } { n } } . \tag { 3 } \end{align*} \end{document}

For a given shape of the volume, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {v_i} = ax_i^3. \tag{4} \end{align*} \end{document}

Therefore, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \sigma _c } \approx \frac { { { w_0 } } } { a } \sqrt { \frac { { \sum \limits_ { i = 1 } ^n \displaystyle { { { \bigg( { \frac { x_i^3 - x_0^3 } { x_i^3x_0^3 } } \bigg) } ^2 } } } } { n } } . \tag { 5 } \end{align*} \end{document}

Because \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$x_i^3 - x_0^3 \approx 3{x_i}{x_0} ( {x_i} - {x_0} ) \approx 3x_0^2 ( {x_i} - {x_0} )$$ \end{document} when \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${x_i} \approx {x_0}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \sigma _c } \approx { \frac { 3 { w_0 } } { ax_0^4 } } \sqrt { \frac { { \sum \limits_ { i = 1 } ^n { { { ( { x_i } - { x_0 } ) } ^2 } } } } { n } } = { \frac { 3 { w_0 } } { ax_0^4 } } { \sigma _x } . \tag { 6 } \end{align*} \end{document}

For a subcellular compartment, assuming its shape is a sphere, the mentioned equation is then simplified as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \sigma _c } \approx { \frac { 9 { w_0 } { \sigma _x } } { 4 \pi x_0^4 } } . \tag { 7 } \end{align*} \end{document}

It is proposed that in Equation (1), when \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \sigma _x} = { \sigma _p}$$ \end{document} , it strikes a balance between the measurement accuracy of the position and that of the momentum of a particle. Under this condition, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \sigma _x } \ge \sqrt { \frac { h } { { 4 \pi } } } . \tag { 8 } \end{align*} \end{document}

Assuming that σ_x in Equation (7) is equivalent to σ_x in Equation (8) since they both represent a standard deviation of distance measurements, and based on the early justification that the uncertainty level of a population of molecules is equal to or more than that of a subatomic particle, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \sigma _c } \ge { \frac { 9 { w_0 } } { 8x_0^4 } } \sqrt { \frac { h } { { { \pi ^3 } } } } . \tag { 9 } \end{align*} \end{document}

To calculate values of σ_c in Equation (9), biologically relevant values of w₀ and x₀ should be determined. Based on the study by Milo and Phillips (2016), it is assumed that, in a subcellular compartment, w₀ = 4.9817 × 10⁻¹⁸ g (the approximate weight of 100 molecules of a protein of 30,000 Da) and x₀ = 1 μm = 10⁻⁶ m, then \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \sigma _c} \ge 8.1928 \times {10^{ - 13}}{ \rm{g}} / {{ \rm{m}}^3} = 8.1928 \times {10^{ - 4}}{ \rm{ng / L}}{ \rm{.}} \end{align*} \end{document}

The mentioned results indicate that when the in vivo concentration values are expressed in ng/L (likely being reasonably large absolute numerical values) in simulation of biochemical reactions, the relative numerical tolerance can be reasonably set at or around 10⁻⁴. For chemical concentrations expressed in a unit other than ng/L, the corresponding relative numerical tolerance values in simulation may be adjusted based on the σ_c value. The proposed range of relative numerical tolerance values suggests that the simulation results shown in Figure 1a and b are likely closer to reality than that shown in Figure 1c.

Based on the mentioned analysis, reducing relative numerical tolerance values for chemical concentration values in simulation will inevitably be met with increased standard deviation values in determining the energy (momentum) state of the molecules. This inverse relationship between the chemical concentrations and the momentum of the molecules likely means that the mathematical relationship among variables in a deterministic model will increasingly suffer inaccuracy with the increase in the accuracy of the concentration values. Furthermore, what is presented here may also be relevant in general to mathematical models that need to account for positions of molecules.