^{1,2}, Rhys Adams

^{3}, Gastone C. Castellani

^{4}and Harel Z. Shouval

^{1}

### Abstract

Many biochemical networks have complex multidimensional dynamics and there is a long history of methods that have been used for dimensionality reduction for such reaction networks. Usually a deterministic mass action approach is used; however, in small volumes, there are significant fluctuations from the mean which the mass action approach cannot capture. In such cases stochastic simulation methods should be used. In this paper, we evaluate the applicability of one such dimensionality reduction method, the quasi-steady state approximation (QSSA) [L. Menten and M. Michaelis, “Die kinetik der invertinwirkung,” Biochem. Z49, 333369 (1913)] for dimensionality reduction in case of stochastic dynamics. First, the applicability of QSSA approach is evaluated for a canonical system of enzyme reactions. Application of QSSA to such a reaction system in a deterministic setting leads to Michaelis-Menten reduced kinetics which can be used to derive the equilibrium concentrations of the reaction species. In the case of stochastic simulations, however, the steady state is characterized by fluctuations around the mean equilibrium concentration. Our analysis shows that a QSSA based approach for dimensionality reduction captures well the mean of the distribution as obtained from a full dimensional simulation but fails to accurately capture the distribution around that mean. Moreover, the QSSA approximation is not unique. We have then extended the analysis to a simple bistable biochemical network model proposed to account for the stability of synaptic efficacies; the substrate of learning and memory [J. E. Lisman, “A mechanism of memory storage insensitive to molecular turnover: A bistable autophosphorylating kinase,” Proc. Natl. Acad. Sci. U.S.A.82, 3055–3057 (1985)]. Our analysis shows that a QSSA based dimensionality reduction method results in errors as big as two orders of magnitude in predicting the residence times in the two stable states.

We would also like to thank our anonymous reviewers for suggestions that have considerably improved this paper, and R. Grima for sharing with us some of his methods. A.A. was supported by the UT graduate initiative, H.Z.S. was partially supported by NIH Grant No. 1RO1DA034970-01, and G.C.C. would like to acknowledge the Fondazione CarisBO.

I. INTRODUCTION

II. ENZYME REACTIONS

III. A BISTABLE KINASE-PHOSPHATASE MOLECULAR SWITCH

IV. DISCUSSION

### Key Topics

- Enzyme kinetics
- 20.0
- Enzymes
- 11.0
- Stochastic processes
- 10.0
- Biochemical reactions
- 9.0
- Mesoscopic systems
- 9.0

##### C12

## Figures

(a) Quasi-steady state approximation. Comparison of deterministic solution of the complete enzymatic system and the reduced system with the quasi-steady state assumption. The dashed line is the substrate trace for the 1D equation with quasi-steady state assumption. For the chosen parameters that satisfy the QSSA validity criterion, the dashed line is reasonably close to the substrate trace for the full model. The substrate trace for both the full and the reduced model settle down to the same steady state (parameters: k 1 = 0.01 [1/(sec* μM)], k −1 = 0.02, k 2 = 0.002, k 3 = 0.001 [1/sec], E T = 25 μM, S T = 100 μM). (b) Stochastic simulation of the chemical master equation for the complete set enzymatic reactions using the Gillespie algorithm. Overlaid on top, is the reduced stochastic simulation by implementing a Gillespie-like algorithm for the reduced chemical master equation, using the quasi-steady state assumption (Eq. (17) ). The mean of the stochastic simulation is well approximated by the deterministic simulation. The distribution of substrate concentration for the reduced stochastic simulation is different than the distribution obtained for the full stochastic model. The volume of the reaction mixture is 10^{−17} l.

(a) Quasi-steady state approximation. Comparison of deterministic solution of the complete enzymatic system and the reduced system with the quasi-steady state assumption. The dashed line is the substrate trace for the 1D equation with quasi-steady state assumption. For the chosen parameters that satisfy the QSSA validity criterion, the dashed line is reasonably close to the substrate trace for the full model. The substrate trace for both the full and the reduced model settle down to the same steady state (parameters: k 1 = 0.01 [1/(sec* μM)], k −1 = 0.02, k 2 = 0.002, k 3 = 0.001 [1/sec], E T = 25 μM, S T = 100 μM). (b) Stochastic simulation of the chemical master equation for the complete set enzymatic reactions using the Gillespie algorithm. Overlaid on top, is the reduced stochastic simulation by implementing a Gillespie-like algorithm for the reduced chemical master equation, using the quasi-steady state assumption (Eq. (17) ). The mean of the stochastic simulation is well approximated by the deterministic simulation. The distribution of substrate concentration for the reduced stochastic simulation is different than the distribution obtained for the full stochastic model. The volume of the reaction mixture is 10^{−17} l.

Comparison of the cumulative distribution of probability of a certain concentration of the substrate in the reaction volume after the corresponding deterministic simulation has reached steady state. The cumulative distribution for the numerical solution of the chemical master equation for the complete reaction system has same mean (full: 29.92, reduced (Eq. (18) ): 30.66, reduced (Eq. (19) ): 29.88, reduced analytical (using Eqs. (13) and (19) ): 30.12) as that of the reduced system but a different standard deviation (full: 2.60, reduced (Eq. (18) ): 8.36, reduced (Eq. (19) ): 3.21, reduced analytical (using Eqs. (13) and (19) ): 3.16). The Kolmogorov Smirnov test shows that the three distributions are significantly different. The analytical and the numerical solution of the reduced chemical master equation (Eq. (17) ) yield identical results.

Comparison of the cumulative distribution of probability of a certain concentration of the substrate in the reaction volume after the corresponding deterministic simulation has reached steady state. The cumulative distribution for the numerical solution of the chemical master equation for the complete reaction system has same mean (full: 29.92, reduced (Eq. (18) ): 30.66, reduced (Eq. (19) ): 29.88, reduced analytical (using Eqs. (13) and (19) ): 30.12) as that of the reduced system but a different standard deviation (full: 2.60, reduced (Eq. (18) ): 8.36, reduced (Eq. (19) ): 3.21, reduced analytical (using Eqs. (13) and (19) ): 3.16). The Kolmogorov Smirnov test shows that the three distributions are significantly different. The analytical and the numerical solution of the reduced chemical master equation (Eq. (17) ) yield identical results.

Relationship between the accuracy of the quasi-steady state approximation in predicting the deterministic transients and the accuracy of the reduced stochastic model to predict steady state distributions. (a) The accuracy of the stochastic simulations quantified by the KL divergence between the full and reduced models as a function of the parameter λ which quantifies the QSSA approximation in the deterministic case. The points denoted by (*) are consistent with the data points in the paper by Barik et al.. ^{17} The circle is around the data point from the example in the previous figures. (b) The accuracy of stochastic simulations quantified as in (a) by the KL distance, as a function of the actual accuracy of the deterministic simulation. Here the accuracy of the deterministic simulations is quantified by the average normalized mean distance (NMSE) between the full and reduced models (see Appendix).

Relationship between the accuracy of the quasi-steady state approximation in predicting the deterministic transients and the accuracy of the reduced stochastic model to predict steady state distributions. (a) The accuracy of the stochastic simulations quantified by the KL divergence between the full and reduced models as a function of the parameter λ which quantifies the QSSA approximation in the deterministic case. The points denoted by (*) are consistent with the data points in the paper by Barik et al.. ^{17} The circle is around the data point from the example in the previous figures. (b) The accuracy of stochastic simulations quantified as in (a) by the KL distance, as a function of the actual accuracy of the deterministic simulation. Here the accuracy of the deterministic simulations is quantified by the average normalized mean distance (NMSE) between the full and reduced models (see Appendix).

The dependence of the relative error between the variance of the full system and the variance of the system reduced using the QSSA approach. Solid lines are the analytical results using the LNA approach (Eq. (26) ), and symbols with error bars are obtained directly from simulations. Different colors represent different parameters sets. In each parameter set all parameters are fixed, except for k 1 which is varied, to obtain different values of ρ.

The dependence of the relative error between the variance of the full system and the variance of the system reduced using the QSSA approach. Solid lines are the analytical results using the LNA approach (Eq. (26) ), and symbols with error bars are obtained directly from simulations. Different colors represent different parameters sets. In each parameter set all parameters are fixed, except for k 1 which is varied, to obtain different values of ρ.

(a) Deterministic simulation for the kinase-phosphatase switch showing bistability. The solid lines are obtained from the simulation of the full system. The dashed lines are from simulations with QSSA. Blue lines are for initial condition K P (0) = 25 μM and black lines are for initial condition K P (0) = 4 μM. (b) The source and sink terms for the quasi-steady state reduced model as a function of phosphorylated kinase concentration. The points of intersection of the solid and dashed lines are the fixed points of the system. (c) Potential energy wells for the fixed points in the reduced model. Note that the fixed points of the reduced model are identical to the fixed points of the full model. (Parameters: r 1 = 0.001 [1/(sec* μM)], r −1 = 0.002, r 2 = 0.02 [1/sec], r 3 = 0.08 [1/(sec* μM)], r −3 = 0.001, r 4 = 0.0539, r 5 = 0.00212 [1/sec], K T = 60 μM, P T = 5 μM.)

(a) Deterministic simulation for the kinase-phosphatase switch showing bistability. The solid lines are obtained from the simulation of the full system. The dashed lines are from simulations with QSSA. Blue lines are for initial condition K P (0) = 25 μM and black lines are for initial condition K P (0) = 4 μM. (b) The source and sink terms for the quasi-steady state reduced model as a function of phosphorylated kinase concentration. The points of intersection of the solid and dashed lines are the fixed points of the system. (c) Potential energy wells for the fixed points in the reduced model. Note that the fixed points of the reduced model are identical to the fixed points of the full model. (Parameters: r 1 = 0.001 [1/(sec* μM)], r −1 = 0.002, r 2 = 0.02 [1/sec], r 3 = 0.08 [1/(sec* μM)], r −3 = 0.001, r 4 = 0.0539, r 5 = 0.00212 [1/sec], K T = 60 μM, P T = 5 μM.)

Stochastic switching between bistable states. (a) Green line is the stochastic simulation of the complete reaction system using the Gillespie algorithm. The initial conditions are set to the upper steady state value. The black horizontal line indicates the two stable steady state concentration levels for each species. Note that the mean of the stochastic fluctuations within an equilibrium state is close to the steady state concentration. (b) The red line is the stochastic simulation using the reduced master equation and a Gillespie-like algorithm. Note the difference in the scale of x axis in the plots on the left and right.

Stochastic switching between bistable states. (a) Green line is the stochastic simulation of the complete reaction system using the Gillespie algorithm. The initial conditions are set to the upper steady state value. The black horizontal line indicates the two stable steady state concentration levels for each species. Note that the mean of the stochastic fluctuations within an equilibrium state is close to the steady state concentration. (b) The red line is the stochastic simulation using the reduced master equation and a Gillespie-like algorithm. Note the difference in the scale of x axis in the plots on the left and right.

Histogram showing difference in upstate and downstate distribution for the full model and quasi-steady state model for stochastic simulations. Kp > 5.7 μM is classified as upstate and Kp < 5.7 μM is classified as downstate.

Histogram showing difference in upstate and downstate distribution for the full model and quasi-steady state model for stochastic simulations. Kp > 5.7 μM is classified as upstate and Kp < 5.7 μM is classified as downstate.

Distribution of residence times in the upstate and downstate for the full model and QSSA model. The total simulation time over which this histogram is calculated is 5 × 10^{8} s. The mean upstate residence time for the full model is 2.70 × 10^{4} s while that for reduced model is 4.93 × 10^{3} s. The mean downstate residence time for full model is 3.25 × 10^{4} s while that for reduced model is 4.70 × 10^{3} s.

Distribution of residence times in the upstate and downstate for the full model and QSSA model. The total simulation time over which this histogram is calculated is 5 × 10^{8} s. The mean upstate residence time for the full model is 2.70 × 10^{4} s while that for reduced model is 4.93 × 10^{3} s. The mean downstate residence time for full model is 3.25 × 10^{4} s while that for reduced model is 4.70 × 10^{3} s.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content