Misuse of the Michaelis–Menten rate law for protein interaction networks and its remedy

For over a century, the Michaelis–Menten (MM) rate law has been used to describe the rates of enzyme-catalyzed reactions and gene expression. Despite the ubiquity of the MM rate law, it accurately captures the dynamics of underlying biochemical reactions only so long as it is applied under the right condition, namely, that the substrate is in large excess over the enzyme-substrate complex. Unfortunately, in circumstances where its validity condition is not satisfied, especially so in protein interaction networks, the MM rate law has frequently been misused. In this review, we illustrate how inappropriate use of the MM rate law distorts the dynamics of the system, provides mistaken estimates of parameter values, and makes false predictions of dynamical features such as ultrasensitivity, bistability, and oscillations. We describe how these problems can be resolved with a slightly modified form of the MM rate law, based on the total quasi-steady state approximation (tQSSA). Furthermore, we show that the tQSSA can be used for accurate stochastic simulations at a lower computational cost than using the full set of mass-action rate laws. This review describes how to use quasi-steady state approximations in the right context, to prevent drawing erroneous conclusions from in silico simulations.


Introduction
The Michaelis-Menten (MM) rate law has been the dominant paradigm for describing the rates of enzyme-catalyzed reactions for 100 years [1][2][3][4]. While the MM rate law was proposed in the pioneering works of Henri [1] and of Michaelis and Menten [2], it was Briggs and Haldane [3] who derived the MM rate law from the underlying molecular mechanism by an approach known as the steady-state approximation. Since then, steady-state approximations have been used to describe other bimolecular interactions, such as reversible binding between a gene and a transcription factor (TF) [5,6] and between a receptor and a ligand [7,8]. In this paper, we refer to the Briggs-Haldane approach as the "standard" quasi-steady state approximation (sQSSA).
Much later, Segel and Slemrod [9] presented a thorough mathematical analysis of the timescale separation that underlies the Briggs-Haldane derivation of the MM rate law. They concluded that the MM rate law is accurate only when the enzyme concentration is low enough so that the concentration of enzyme-substrate complex is negligible compared to substrate concentration [9]. This constraint is acceptable for most metabolic reactions in living cells, where substrate concentrations are typically much larger than the concentrations of enzymes that catalyze the reactions of intermediary metabolism. However, for protein interaction networks that govern many aspects of cell physiology, the "enzymes" and "substrates" are proteins (kinases, phosphatases, TFs, stoichiometric inhibitors, etc.) that are typically present in cells at comparable concentrations [10][11][12][13]. Therefore, it is at least suspicious and at worst misleading to use the MM rate law to model enzyme-catalyzed reactions in protein interaction networks. Unfortunately, even in such situations outside its range of validity, the MM rate law has frequently been (mis)used to model enzyme-catalyzed reactions.
In this review, we illustrate how such unjustified use of the MM rate law for protein interaction networks leads to erroneous conclusions. We describe how such errors can often be prevented by replacing the MM rate law with a slightly modified form, based on the "total" quasisteady state approximation (tQSSA), which is generally accurate for any combination of substrate and enzyme concentrations [14][15][16][17][18][19][20][21]. Specifically, when enzyme concentration is not low enough for enzyme-substrate complex to be negligible compared to the substrate, the tQSSA (but not the MM rate law) can lead to accurate estimation of enzyme kinetic parameters. When the MM rate law is used, zero-order ultrasensitivity and bistability are predicted in cases when they are impossible, but the tQSSA gets it right. Furthermore, a model of a transcriptional negative feedback loop, based on the sQSSA, predicts "no oscillations," but the same model, based on the tQSSA, is consistent with oscillatory dynamics. Finally, we illustrate that the tQSSA can accurately capture stochastic dynamics of enzyme kinetics, zero-order ultrasensitivity, and oscillations.
This review is written specifically for computational biologists who build, analyze, and simulate mathematical models (deterministic and/or stochastic) of protein interaction networks, to warn them of the dangers of misusing the MM rate law and to provide a simple fix, based on the tQSSA. Our comparisons of sQSSA and tQSSA models are mainly based on simulations. For readers who are more interested in the experimental consequences or the mathematical/statistical properties of sQSSA versus tQSSA, we provide some references to the literature on these topics.

Quasi-steady state approximation for enzyme-catalyzed reactions
A single-substrate enzyme-catalyzed reaction (Fig 1a) can be described by the following system of ordinary differential equations (ODEs) based on mass-action kinetics: where S, E, C, and P are the time-dependent concentrations of substrate (S), enzyme (E), enzyme-substrate complex (C), and product (P), respectively. Note that we use a different font style to distinguish the name of a substance from its concentration throughout this review. Because dE dt ¼ À dC dt ; the total enzyme concentration (E T � C+E) is constant, as it should be. The MM mechanism for a single-substrate enzyme-catalyzed reaction; S, substrate; E, enzyme; C, enzyme-substrate complex; P, product. Note that we use a different font style to distinguish a substance, S, from its concentration, S in Eq 1. The full model, based on mass-action kinetics (Eq 1), can be simplified by either the sQSSA (Eq 3) or the tQSSA (Eq 5). Both approximations assume that C rapidly reaches a QSSA, which is solely determined by either S (Eq 2) orŜ ¼ S þ C (Eq 4). (b-c) Heat maps of the predicted maximal fraction of total enzyme that is in complex with substrate, i.e., C(S T )/E T , predicted with sQSSA (Eq 2) and tQSSA (Eq 4). In these panels, color represents the ratio C(S T )/E T , and gray dashed lines are contours of E T /(K M + S T ) = 0.1 and 1.0. The sQSSA (Eq 2) predicts that this fraction is independent of total enzyme concentration, E T (e.g., C(S T )/E T = 0.5 when S T = K M regardless of E T ) in panel b. On the other hand, the tQSSA (Eq 4) predicts that, as E T increases, more substrate is needed to achieve that same level of substrate saturation of enzyme molecules in panel c. When the validity condition for the sQSSA holds (i.e., below the gray dashed line where E T /(K M + S T ) �1), the predictions of the two models are nearly identical. (d-e) The simulated accumulation of P and the estimation of the Michaelis constant, K M , from the progress curve of P. For these calculations, we chose k f = 10, k b = 9, k cat = 1, K M = 1, S T = 9, and E T /(K M + S T ) = 0.1 (panel d; triangle mark in (b) and (c)) or 1 (panel e; star mark in (b) and (c)). Initial conditions are S(0) = S T , C(0) = 0, P(0) = 0. If E T is low (i.e., E T /(K M + S T ) = 0.1), the simulation of P and the estimation of K M with the sQSSA model are accurate (d). On the other hand, if E T is high (i.e., E T /(K M + S T ) = 1), the simulated P with the sQSSA model using the true rate constants is wildly divergent from the true progress curve simulated with the full model (e). As a result, if k cat and K M are adjusted to fit the sQSSA model to the true progress curve, then the estimates of k cat and K M are far from the true values (inset; the estimates of k cat is not shown). On the other hand, the tQSSA model accurately reproduces the progress curve of P and estimates the rate constants regardless of the value of E T (d and e). Inset: gray lines are the true value of K M , and the red and blue circles are the mean K M estimates from the tQSSA and sQSSA models, respectively. The bars represent 99% confidence intervals. The Quasi-Newton method is applied to 21 sampled points of the progress curve of P simulated with the full model to identify the value of K M that minimizes the fitting error of the sQSSA and tQSSA models. The true values of parameters are used as Furthermore, because the total concentration of substrate + product (S T � S+C+P) is also constant, Eq 1 describes the dynamics of only two independent variables, say C and P. This twovariable model can be simplified to a single-variable model by a "quasi-steady state" approximation (QSSA), which eliminates C under the assumption that C rapidly reaches a quasisteady state (QSS) where dC dt � 0, and thereafter C slowly tracks along the QSS, often referred to as the "slow manifold." A good approximation to the QSS (i.e., the QSSA) can be obtained by solving dC dt ¼ 0: where is the "Michaelis constant" of the enzyme. By substituting the QSSA (Eq 2) for C in the full model (Eq 1), we can describe the rate of product accumulation as which is known as the MM rate law (Fig 1a). Of course, because S is continually changing with time as the reaction proceeds, S needs to be calculated from the conservation condition, S T � S +C+P (see below). The QSSA (Eq 2) implies that CðSÞ , that enzyme molecules are halfsaturated with substrate when S = K M , regardless of total enzyme concentration, E T (Fig 1b). This seems unlikely because as E T increases, more substrate is likely to be needed for half-saturation of enzyme molecules, suggesting that the QSSA is likely to be valid only if E T is small enough. Indeed, Segel and Slemrod showed that the QSSA is valid only when E T � S T +K M [9]. Under this condition, implying that the complex is negligible compared to the total amount of substrate. Thus, C can be neglected in the expression for total substrate concentration, S T � S+C+P, and S can be replaced with S T −P in Eq 3, leading to a single-variable model. We refer to this approximation (Eq 2) as the sQSSA [9], and we refer to Eq 3 (the MM rate law) as the sQSSA model of a single-substrate enzyme-catalyzed reaction. The constraint of low enzyme concentration (E T � S T +K M ), which is required to use the sQSSA, is acceptable for enzyme-catalyzed reactions in intermediary metabolism, where enzyme concentrations are typically 1,000-fold smaller than metabolite (substrate) concentrations.
On the other hand, in protein interaction networks, both enzyme and substrate are proteins whose intracellular concentrations are roughly comparable [10][11][12][13]. In this case, because C does not reach the sQSSA (Eq 2) on a short timescale compared to changes in S [9], we need to analyze the kinetics of the reaction by a different approximation, the tQSSA [14][15][16][17][18][19][20][21]. For the tQSSA, S is replaced with the total substrate concentrationŜ � S þ C: Now, the validity condition for the QSSA is changed from "C reaches QSS before S changes appreciably" to "C reaches QSS beforeŜ changes appreciably." This new condition is easier to satisfy than the old condition because, compared to S,Ŝ (the sum of S and C) is less likely to change substantially until the product (P) has accumulated substantially (see [19,22] for detailed timescale analysis). What a brilliant idea "to make C the winner of the race to QSS, change its competitor from a rabbit (S) to a turtle (Ŝ)"! (Of course, the rabbit (S) is as slow as the turtle (Ŝ) when the validity condition for the MM rate law holds, and thusŜ � S.) the initial condition for the fitting procedure in order to avoid parameter unidentifiability. MM, Michaelis-Menten; sQSSA, standard quasi-steady state approximation; tQSSA, total quasi-steady state approximation. https://doi.org/10.1371/journal.pcbi.1008258.g001 To make this change, we rewrite the mass-action rate expressions (Eq 1) in terms ofŜ instead of S: which has the same underlying molecular mechanism and dynamics as the original model (Eq 1). Then, the tQSSA is derived in terms ofŜ (rather than S) by solving the quadratic equation for dC dt ¼ 0: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi CðŜÞ is not as elegant as the C(S) derived with sQSSA (Eq 2), but it is just as easy to implement in a computer program. A validity condition for the tQSSA (Eq 4) was derived by Tzafriri [16]: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Tzafriri also showed that ε tQ S T ; E T ð Þ � ε tQ 0; K M ð Þ ¼ k cat 4ðk b þk cat Þ : Thus, if k cat � k b , then ε tQ � 1 regardless of total enzyme concentration, E T . Furthermore, Tzafriri also showed that ε tQ � 1 if either S T or E T (or both) are much greater than K M . Importantly, because ε tQ � k cat 4ðk b þk cat Þ < 0.25 always holds, the condition ε tQ � 1 is not so badly violated for any cases. This indicates that the tQSSA is generally a good approximation for any ratio of enzyme concentration to substrate concentration and for any ratio of timescales, thanks to the simple but amazing idea of "changing C's competitor from S toŜ." We refer readers to [22][23][24] for stricter validity conditions for sQSSA and tQSSA than those presented in this review, conditions that ensure for the long-time validity of the approximation with a rigorous error analysis.
Unlike CðSÞ=E T ; CðŜÞ=E T has the reasonable consequence that, as E T increases, more substrate is needed to half-saturate the enzyme (cf . Fig 1b and 1c). This can be seen more clearly with the Padé approximant for CðŜÞ: [10,16,17]. According to the Padé approximant, enzyme molecules are half-saturated with substrate (i.e., CðŜÞ E T ¼ 1 2 ), whenŜ = K M + E T . That is, the total substrate concentration needed for half-saturation of enzyme increases linearly with total enzyme concentration, E T (Fig 1c). Furthermore, when the sQSSA is valid (i.e., E T � S T + K M ), the Padé approximant becomes nearly identical with C(S) derived with sQSSA (Eq 2). This explains nearly identical values of C(S) and CðŜÞ when E T � S T + K M (Fig 1b and 1c, below the gray dashed line where E T S T þK M ¼ 0:1). By substituting CðŜÞ (Eq 4) into the full model, we can describe the rate of product accumulation in terms ofŜ rather than S, and replace the MM rate law (Eq 3) by ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi whereŜ ¼ S T À P (Fig 1a). Note that the conservation of S T is exact in this tQSSA model, unlike the MM rate law (Eq 3), where the complex is neglected in the conservation condition.
Reflecting the validity of tQSSA over a wide range of conditions, simulations of the temporal accumulation of product with the tQSSA model (Eq 5) are consistent with simulations of the full model (Eq 1) regardless of enzyme concentration (Fig 1d and 1e). On the other hand, the sQSSA model (Eq 3) is accurate only when enzyme concentration is low enough compared to S T + K M , i.e., when the sQSSA is valid (Fig 1d and 1e).

Estimating the rate constants of an enzyme-catalyzed reaction
An accurate model is required for accurate estimation of parameters based on model fitting [13,25,26]. Indeed, in the case of low enzyme concentration, both the sQSSA and the tQSSA models provide good estimates of the "true" rate constants when the models are fitted to the progress curve P(t) simulated with the full model (i.e., a "progress-curve assay") (Fig 1d inset). However, for the case of high enzyme concentration, only the tQSSA model provides accurate estimates of the rate constants (Fig 1e inset). Here, to avoid potential unidentifiability of the parameters and thus focus on the bias of estimation caused by model inaccuracy, we used the true parameter values as the initial condition of fitting, under the assumption that we have good prior knowledge. However, in the absence of prior knowledge, the validity of the model reduction (i.e., forward problem) does not guarantee the identifiability of the parameter (i.e., inverse problem) [13]. For instance, if model solutions give equally acceptable fits to the constraining experimental data for several different sets of parameter values, then the "true" parameter values are unidentifiable. Indeed, kinetic parameters only can be estimated with the MM rate law under a stricter condition than its validity condition [27]. Although the identifiability of kinetic parameters with the tQSSA has rarely been explored, a recent numerical study used a Bayesian approach to show that the identifiability of kinetic parameters is not hampered by using the seemingly more complex tQSSA model (Eq 5) instead of the MM rate law for a progress-curve assay [26]. Importantly, due to the unbiased estimation afforded by the tQSSA model regardless of enzyme concentrations, an optimal experiment to identify parameters can be easily designed [26]. However, before the tQSSA model can be used to interpret real experimental data, further experiments and mathematical analysis are needed to validate this in silico study. When the MM rate law is valid, it should be used for parameter estimation (rather than the tQSSA model) because the MM rate law has been extensively investigated and validated in myriad experimental circumstances [13,27,28]. Furthermore, with the "initial-velocity assay" (based on the initial rate of product accumulation), the estimation of kinetic parameters with the MM rate law is straightforward (e.g., K M = substrate concentration at which the initial velocity is half-maximal) [13]. Thus, when the MM validity condition holds, there is no need to consider the more complex tQSSA model instead of the MM rate law.
However, the constraint on using the MM rate law, E T � S T + K M , cannot always be satisfied even for in vitro experiments, e.g., a low concentration of substrate is required for sensitive QD-FRET-based probes [29][30][31]. Importantly, in vivo enzyme concentrations can vary extensively inside a cell and are usually much higher than those used experimentally in vitro [11][12][13]32]. For instance, the enzyme cytochrome P450, which metabolizes many drugs in the liver, is present in concentrations greatly exceeding [drug] + K M [33,34]. However, the MM rate law has been used to predict the rate of drug metabolism in liver and thus drug clearance in approximately 60,000 publications [35]. A recent study describes the errors caused by the MM rate law for predicting hepatic drug clearance and shows how the prediction can be improved by the tQSSA model [36]. However, further study is needed for the application of the tQSSA model to in vivo systems where enzyme concentrations might keep changing dynamically [19].

Ultrasensitivity and bistability
Ultrasensitivity is a critical property of many signaling networks because it amplifies signals in a nonlinear manner that can trigger a switch-like response [37][38][39][40]. Here, we illustrate that, when the MM rate law is misused in models of protein interaction networks, the model may falsely predict ultrasensitivity and bistability.
In the Goldbeter-Koshland (GK) mechanism of "zero-order ultrasensitivity," a substrateproduct pair (S and S P ) are interconverted by two enzymes (E and D) (Fig 2a) [41]. This system is commonly observed in phosphorylation and dephosphorylation processes: S is phosphorylated by kinase E, and S P is dephosphorylated by phosphatase D, which can be described by the following system of ODEs: where concentrations of total substrate (S T � S + S P + ES + DS P ), total kinase (E T � E + ES), and total phosphatase (D T � D + DS P ) are all conserved quantities. As the GK mechanism is composed of two interlocked, single-substrate, enzyme-catalyzed reactions (Fig 1a), it can be simplified by using the QSSA as done in the previous section (Fig 2a).
In particular, assuming low concentrations of total kinase and phosphatase (E T ,D T � S T ), Goldbeter and Koshland used the sQSSA for ES and DS P (Eq 2) to describe the rates of phosphorylation and dephosphorylation in terms of S and S P : where K ME ¼ k be þk e k fe and K MD ¼ k bd þk d k fd (Fig 2a). Because the sQSSA assumes that ES and DS P are negligible, they can be neglected in the conservation condition S T � S + S P + ES + DS P , and thus S can be replaced with S T − S P in Eq 7. In this way, Goldbeter and Koshland derived a simple one-variable model, which we refer to as the sQSSA model of the GK mechanism.
Alternatively, the rates of phosphorylation and dephosphorylation can be described using the tQSSA for ES and DS P (Eq 4) [10,19,25,43,44]. This leads to the tQSSA model of the GK mechanism (Fig 2a): ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi In the GK mechanism, substrate (S) is phosphorylated by kinase (E), and then phosphorylated substrate (S P ) is dephosphorylated by phosphatase (D). The full model based on mass-action kinetics (Eq 6) can be simplified by replacing the concentrations of the substrate-kinase complex (ES) and phosphorylated substrate-phosphatase complex (DS P ) with either the sQSSA (Eq 7) or the tQSSA (Eq 8). In the tQSSA model,Ŝ � S þ ES andŜ P � S P þ DS P are total unphosphorylated substrate and phosphorylated substrate concentration, respectively. (b-d) As the total concentration of kinase (E T ) increases, the steady-state concentration of phosphorylated substrate (Ŝ P ) increases. The sQSSA model predicts a steep sigmoidal response ofŜ P regardless of the level of D T (in this case,Ŝ P ¼ S P , as complexes are assumed to be negligible). On the other hand, both the full and the tQSSA models predict that the steep response is lost as D T increases. In these calculations, S T = 100, k fe = k fd = 10, k d = k e = 1.7, and K ME = K MD = 1. (e-g) Heat maps of the predicted effective Hill exponent of the response function. When the validity of the sQSSA is not satisfied (i.e., above the dashed gray line), the tQSSA model (e) but not the sQSSA model (f) captures the results of the full model (g).
Here, the effective Hill exponent is estimated using log [81]/log[EC90/EC10] [42]. The circle, triangle, and star mark represent the parameter values used for (b), (c), and (d), respectively. (h) Synthesis and degradation of the kinase (E) are added to the sQSSA and tQSSA models for the GK mechanism (a). In particular, the addition ofŜ P -dependent synthesis of E (k pos ) creates a positive feedback loop betweenŜ P and E. (i-j) In panel i, the steady states of the extended positive feedback models (circles) are intersections of the nullclines of E T defined by dE T dt ¼ 0 (solid black lines) and of the original models of the GK mechanism defined by dŜ P dt ¼ 0 (colored lines;Ŝ P ¼ S P for the sQSSA model as complexes are assumed to be negligible). As the nullcline of E T changes, depending on the basal synthesis rate of E (k basal ), so do the steady states in panel i, which is also illustrated as a bifurcation diagram in panel j. The filled and open circles represent stable and unstable steady states, respectively. In the bifurcation diagram, the solid and dashed lines denote the stable and unstable steady states, respectively. Here, the same parameter set is used as in (d), with k pos = k deg = S À 1 T = 0.01. GK, Goldbeter-Koshland; sQSSA, standard quasi-steady state approximation; tQSSA, total quasisteady state approximation.
https://doi.org/10.1371/journal.pcbi.1008258.g002 whereŜ P � S P þ DS P andŜ � S þ ES. Using the conservation condition S T ¼ ðS þ ESÞ þ ðS P þ DS P Þ 1 S þŜ P ;Ŝ can be replaced with S T ÀŜ P in the above equation, leading to a one-variable tQSSA model. Note that the complexes (ES and DS P ) are not neglected in the conservation condition unlike the sQSSA model.
As the concentration of kinase (E T ) increases, the phosphorylation reaction is more likely to occur, and so the steady-state fraction of phosphorylated substrate (S P /S T ) increases. For the sQSSA model (Eq 7), Goldbeter and Koshland derived an explicit function for S P /S T in terms of E T [10,41]. They showed that for K ME , K MD � S T , the function is sigmoidally shaped with an abrupt jump from S P /S T � 0 to S P /S T � 1 at k e E T � k d D T , i.e., when the maximum rates of phosphorylation and dephosphorylation of substrate are equal (Fig 2b). This behavior is known as "zero-order" ultrasensitivity because the kinase and phosphatase enzymes are-for the most part-saturated with substrate molecules when K ME , K MD � S T (i.e., tight binding). A similar sharp transition, close to k e E T � k d D T , is also seen for the steady-state fraction of total phosphorylated substrate (Ŝ P =S T ) in the full model (Eq 6) and tQSSA model (Eq 8) when the enzyme concentrations are small compared to total substrate concentration (E T , D T � S T ;  Fig 2b and 2e-2g, green circle). On the other hand, they make different predictions as enzyme concentrations increase, and the validity condition for the sQSSA is violated (Fig 2c-2g, green triangle and star). In this case, while both the full and the tQSSA models predict less sensitive responses, the sQSSA model continues to predict ("phantom") ultrasensitivity.
Embedding an ultrasensitive response in a positive feedback loop is a common mechanism for creating a bistable switch [37,38]; consequently, modeling zero-order ultrasensitivity based on MM rate laws (Eq 7) in a protein interaction network can generate "phantom" bistability as well. To illustrate this, we embed the GK mechanism (Fig 2a) in a positive feedback loop, where both forms of phosphorylated substrate (Ŝ P ¼ S P þ DS P ) act as TFs for E, the kinase that converts S to S P (Fig 2h). Specifically, the model of the GK mechanism (Eq 6) is extended to model a positive feedback loop by including the following differential equation for total kinase (E T ): The extended model of a positive feedback loop can be simplified using either the sQSSA (Eq 7) or the tQSSA (Eq 8), whereŜ P ¼ S P for the sQSSA model (as DS P is assumed to be negligible) orŜ P ¼ S P þ DS P for the tQSSA model.
The steady states of the sQSSA and tQSSA positive feedback models are found at the intersections of the nullclines of E T (Fig 2i, black solid lines) and the models of GK mechanisms (Fig 2i, colored lines for either the full, the sQSSA, or the tQSSA model of the GK mechanism). To illustrate the problem of the sQSSA model, we consider the case of high enzyme concentrations (D T /S T = 1), where the sQSSA and tQSSA models predict different response curves (blue and red) for the GK mechanism (Fig 2d and 2i). Both the tQSSA and the full models have a single stable steady state regardless of k basal (Fig 2i and 2j, filled red circles). On the other hand, when k basal is between 0.19 and 0.82, the sQSSA model predicts two stable states and an unstable state (filled and open blue circles, respectively, in Fig 2i and 2j) due to the "phantom" ultrasensitivity of the sQSSA model (Fig 2d). In this case, the sQSSA model predicts "phantom" bistable behavior, whereas the full and the tQSSA models are monostable.
Such misuse of the MM rate law infected early models of the eukaryotic cell cycle control system proposed by Tyson, Novak, and colleagues [45][46][47][48]. They routinely employed sQSSA models of zero-order ultrasensitivity (Eq 7) in their protein interaction networks governing G1/S and G2/M transitions in the cell cycle, and they stressed the idea that bistability plays a crucial role in these transitions. Although the bistable behavior predicted by these models was confirmed in a number of independent experimental studies [49][50][51], the theoretical foundation of this behavior was shaky because bistability predicted by the sQSSA could not be trusted. Indeed, in later publications, it was shown that bistability is impossible in some of the Novak-Tyson models when the sQSSA is "unpacked" into full mass-action kinetics [52]. Ciliberto and colleagues [10] addressed this problem by comparing sQSSA, tQSSA, and full models of some representative protein interaction networks. First, they considered a simple model of two protein kinases, S and E, that mutually phosphorylate and inactivate each other. They showed that bistability predicted by the sQSSA model is spurious because the full model cannot exhibit bistability for any values of the kinetic rate constants. Bistability can easily be recovered by assuming that S P retains some catalytic activity for phosphorylating E, a fact that is correctly predicted by the tQSSA model. Then, they extended their first simple model to reflect the interactions among the major regulators (MPF, Wee1, and Cdc25) of the G2/M transition in the cell cycle. For biochemically realistic values of the rate constants, their updated models predict bistability, which is reassuring given that this transition is known to be bistable [50,51]; however, time-course curves of the transition from G2 into M simulated by the full model are correctly reproduced by the tQSSA model but not by the sQSSA model [10].
In addition to ultrasensitive responses and bistability, other dynamic properties of signaling networks are more accurately modeled by tQSSA than by sQSSA [19]. For instance, a careful study of phosphorylation and dephosphorylation cycles in the MAP kinase pathway identified various types of input-output characteristics and their noise filtering function with tQSSA, which were not possible with sQSSA [53].

Biochemical oscillations
To generate periodic oscillations, a negative feedback loop is required [54][55][56][57][58]. For instance, many synthetic oscillators in bacteria are based on transcriptional negative feedback, where a repressor protein directly binds with the promoter of its own gene to suppress its transcription [59][60][61][62][63]. In this case, we can expect that the number of repressor molecules (100s or 1,000s) is in large excess over gene promoter regions (1 or 2) [64,65]. Thus, the repressor-promoter complex is negligible compared to total repressor, and the validity condition for the sQSSA is satisfied. Consequently, MM-or Hill-type functions, derived with the sQSSA, can be used to describe the probability that gene promoter regions are not bound with repressors [66][67][68].
In other cases, an intermediate molecule is involved in the transcriptional negative feedback loop [54]. For instance, repression of transcription can occur via protein sequestration [69][70][71], when a repressor sequesters a transcriptional activator in a 1:1 stoichiometric complex rather than directly binding to its own gene promoter (Fig 3a). To describe the binding of repressor and activator proteins whose concentrations can be comparable, the sQSSA is not appropriate, and the tQSSA should be used.
To illustrate the problem, let's compare sQSSA and tQSSA models of a protein sequestration-based transcriptional negative feedback loop for mRNA (M), encoded protein (P), and repressor protein (R) (Fig 3b): where the ratio of free activator to total activator, AðRÞ=A T , is to be given by either the sQSSA or the tQSSA, as displayed in Fig 3a.R ¼ R þ RA for the tQSSA model andR ¼ R for the sQSSA model as the complex (RA) is assumed to be negligible. In the sQSSA model, free activator concentration is given by A R ð Þ ¼ A T K d K d þR , and so the rate of transcription of mRNA is a M K d K d þR (Fig 3a left). This model is a special case of Goodwin's classic model of periodic gene expression, where transcriptional suppression is described by a Hill function a M K n d K n d þR n (i.e., when n = 1) [72]. It is well known that Goodwin's model exhibits oscillations only if n > 8 [67,73,74]. That is to say, for the sQSSA model with a M K d K d þR , the feedback effect does not have a degree of cooperativity large enough to generate oscillations. Indeed, the sQSSA model cannot produce oscillations for any choice of parameter values (Fig 3c).
The tQSSA model, known as the Kim-Forger model, was developed to describe the circadian clock in mammalian cells [70,71,75,76]. Unlike the sQSSA model, it exhibits sustained oscillation for comparable concentrations of total activator (A T ) and the total repressor (R), provided K d � A T , i.e., tight binding between them (Fig 3d). The full model, based on massaction kinetics, exhibits oscillations that are nearly identical to the tQSSA model, provided fast, reversible, tight binding between repressor and activator (Fig 3e). Both the full and the tQSSA models exhibit oscillations in a region of parameter space where the sQSSA model is invalid, i.e., above the gray dashed line in Fig 3d and 3e. On the other hand, all three models agree that oscillations are impossible when activator concentration is low (below the gray dashed line in Fig 3c-3e). Once again, the sQSSA model is accurate when its validity condition holds but can cause errors outside of its range of validity. In this particular case, the sQSSA model destroys the ultrasensitive response of transcription rate to repressor accumulation that is required to generate oscillations [54,58,70,71]. Unlike the sQSSA model, the tQSSA model predicts oscillations when the amounts of activator and repressor are similar, and they bind tightly (Fig 3d). This is the case for the mammalian circadian clock: repression of transcription mainly occurs via protein sequestration [77][78][79] with a tight binding [80], and the amounts of the key activator (BMAL1) and repressor (PER1/2) are similar [81,82]. Importantly, as their ratio becomes closer to 1:1, the amplitude and stability of circadian rhythms increase [77,83], consistent with predictions of the tQSSA model (Fig 3d) [70,71]. Other properties of the mammalian circadian clock (e.g., robust network design) are also accurately captured by the tQSSA model [71]. Because the amounts of activator and repressor are similar in the mammalian circadian clock, using the sQSSA model would be inappropriate. However, mechanisms for which the sQSSA model would be appropriate, such as incorporation of an additional positive feedback loop or cooperative multisite phosphorylation, have been proposed as key oscillatory mechanisms for the mammalian circadian clock, but experimental evidence for these mechanisms are not presently available (see [71] for details).

Stochastic total quasi-steady state approximation
The chemical master equation (CME) is a popular way to capture the stochasticity of biochemical reactions. In the presence of timescale separation, the CME can be simplified by assuming that fast species are in QSS. The stochastic QSSA of a fast species is its conditional average number given the numbers of molecules of the slow species [84][85][86][87]. However, calculation of the conditional average has been possible so far in very limited circumstances, such as linear reaction kinetics, simple birth-death processes, feed-forward networks, and complex-balanced networks [88][89][90][91]. In the hope of making progress in the general case, modelers have approximated the stochastic QSSA with the more easily derived deterministic QSSA obtained from the corresponding deterministic ODE system [92][93][94][95]. Specifically, nonelementary reaction rates based on the deterministic QSSA (e.g., the MM function) are converted to stochastic propensity functions by converting the concentration of a species, [X], to its number of molecules, N X , with the relationship N X = [X]�O, where O is the volume of the system. Then, based on the nonelementary propensity functions, stochastic simulations are performed with Gillespie's algorithm [92][93][94][95]. In this way, a deterministic model simplified with QSSA can be used to perform stochastic simulations.
This heuristic approach has been widely used in the hope that simulations of the simplified stochastic model will be accurate; however, this is often not the case when the deterministic sQSSA is used to calculate propensity functions [96][97][98][99]. On the other hand, when the deterministic tQSSA (e.g., Eq 4) is used, stochastic simulations have been surprisingly accurate with much lower computation cost [26, [98][99][100][101][102][103][104][105]. To illustrate, consider first the case of a single-substrate enzyme-catalyzed reaction (Fig 1a). In Fig 4a inset, we simulate the catalysis-completion time τ, i.e., the time elapsed for all of the initial substrate to be converted to product. This time varies from one simulation to another, of course, because of stochastic fluctuations in chemical reactions involving small numbers of molecules. From many realizations of the stochastic process, the distributions of τ for the full and tQSSA models are calculated, and they are nearly identical (Fig 4a). Similarly, fluctuations of the oscillatory period (peak-to-peak duration) of a transcriptional negative feedback loop (Fig 3a) are also accurately captured by the stochastic tQSSA model (Fig 4b). As a third example, simulations performed with the stochastic tQSSA model of the GK mechanism (Fig 2a) are accurate both when D T /S T � 1 (Fig 4c) and when D T /S T = 1 (Fig 4d). When D T /S T � 1, the GK mechanism exhibits zero-order ultrasenstivity, i.e., the fraction of total phosphorylated substrate changes sharply around E T = D T (Fig  2b). Reflecting this sharp change, the fraction of the total phosphorylated substrate (Ŝ P =S T ) fluctuates wildly in the stochastic simulations, and the stationary distribution ofŜ P =S T is broad (Fig 4c). On the other hand, when D T /S T = 1, the fraction of total phosphorylated substrate changes very little around E T = D T (Fig 2d), and thus its stationary distribution is narrow (Fig 4d).
The accuracy of the stochastic simulation using the deterministic tQSSA stems from the accuracy of the deterministic tQSSA over a wide range of conditions [99]. In the presence of timescale separation, a deterministic solution evolves in two phases: an initial transient phase and a QSS phase. The two phases correspond, respectively, to times before and after the solution is in QSS (i.e., "on the slow manifold"). During an initial transient phase, the fast variable rapidly reaches the slow manifold. Then, the system evolves along the slow manifold during the QSS phase. By contrast, the trajectory of the stochastic system, unlike the trajectory of the deterministic system, keeps jumping off the slow manifold due to fluctuations. Thus, unlike the deterministic system, transient phases continually recur whenever the stochastic system leaves the slow manifold due to a fluctuation. Importantly, depending on exactly how the trajectory escapes the slow manifold, the transient phase will begin from different initial conditions. Thus, for accurate stochastic simulations, the deterministic QSSA must be accurate not only for specific initial conditions but also for a range of initial conditions that cover the most likely fluctuations away from the slow manifold. For this reason, the tQSSA is a better choice than the sQSSA for stochastic simulations, because the tQSSA is accurate over a much wider range of initial conditions than the sQSSA (see [99] for details).
However, in the presence of large fluctuations, the moment closure assumption underlying the replacement of the stochastic QSSA with the deterministic tQSSA [98,103,105] can cause considerable errors [90]. In this case, model reduction based on exact moment derivation without the moment closure assumption is needed [90], or a recently developed approximation based on curvature of the slow manifold and initial transient flow field can be used [106]. See [87] for the other types of approximations for CMEs with timescale separation. If all else fails, the full system of mass-action rate laws can be simulated by Gillespie's stochastic simulation algorithm.

Conclusion
The MM rate law has been a great asset in describing enzyme-catalyzed reactions for more than a century. It is valid for most metabolic reactions in living cells, where substrate concentrations are much larger than the concentrations of enzymes that catalyze the reactions of intermediary metabolism. However, when used outside of its domain of validity (low enzyme concentrations), the MM rate law can cause serious errors of dynamics, thereby reducing the reliability of MM-based models (Figs 1-3). Of most relevance in this regard are protein interaction networks, for which enzymes and their substrates are usually present in comparable concentrations; hence, the validity condition for the MM rate law is in doubt, and the conclusions of models based on the MM rate laws are seriously compromised.
It is often difficult to check the validity condition (E T � S T + K M ) of the MM rate law due to limited information about the underlying system. For instance, it is challenging to measure in vivo enzyme concentrations, which might be highly heterogeneous across cell types. For the case of gene regulation by TF binding to gene regulatory sequences, an MM-type rate law is likely to be valid because TFs are usually in large excess over their chromosome-binding sites [64,65]. However, when intermediate molecules are involved, such as cofactors that mediate the interactions between TFs and gene promoter regions, it is risky to use an MM rate law (a hyperbolic function for the probability of gene transcription) to describe such gene regulatory networks (Fig 3). Similarly, it is risky to use an MM-type function to describe the interaction between ligands and receptors, when their concentrations are comparable [8]. In all these situations, we need an alternative to MM-type rate laws.
Classically, the MM rate law was derived by Briggs and Haldane [3] by assuming that the enzyme-substrate complex rapidly reaches a quasi-steady state concentration determined by the substrate concentration. Their "standard" quasi-steady state approximation (sQSSA) is valid if the enzyme concentration is low enough so that the concentration of enzyme-substrate complex is negligible compared to substrate concentration. In cases where this proviso is not satisfied, there is a simple fix: to re-derive the rate law based on the "total" quasi-steady state approximation (tQSSA). The tQSSA is better than the sQSSA because it is accurate for a much wider range of conditions, as exemplified in Figs 1-3. Furthermore, the tQSSA is accurate even for stochastic simulations (Fig 4). Admittedly, the tQSSA rate law (Eq 4) is less intuitively appealing than the MM rate law (Eq 2) and less amenable to mathematical analysis, but for computer simulations, there is no drawback to using the tQSSA in place of an MM rate law.
In light of the uncertainties surrounding sQSSA models, why has this approach been so popular among modelers, including (sorry to say) the senior author of this review? Well, the MM rate law is well known and respected by molecular cell biologists, so, of course, it seemed like a good place to start. Signaling pathways often exhibit highly sigmoidal responses [38], and "zero-order ultrasensitivity" (which is based on the MM rate law) is a simple and elegant way to model such behavior. When problems with the sQSSA approach were recognized, the suspicious models were justified as "phenomenological" rather than "mechanistic." As a last resort, one could always fall back on the "spherical cow" defense, i.e., it's OK to make unrealistic assumptions if they capture the "relevant details" of the problem.
Indeed, it is a common and respectable practice to simplify complex systems by neglecting irrelevant details. However, according to Kruskal's minimum simplification principle [108], "no term should be neglected without a good reason." When the sQSSA is used to simplify a complex reaction mechanism without good reason (e.g., when its validity condition is not satisfied), the resulting model can cause errors of both steady-state and transient behaviors [109][110][111][112]. In particular, when the MM rate law is used outside its domain of validity, the combination of faulty QSS assumptions and neglect of enzyme-substrate complexes causes serious errors of dynamics, thereby reducing the reliability of MM-based models, as illustrated in Figs 1-3.
One must bear in mind that the tQSSA also has validity constraints: it is invalid if k cat � k b and E T + S T is not in excess over K M [16,24]. These constraints are not likely to limit tQSSA's applicability to protein interaction networks, but they must be checked if we are to take Kruskal's principle seriously.
In light of the uncertainties involved in any QSSA, why not use the full mass-action kinetic equations (e.g., Eq 1)? Because the QSSA reduces the dimension of a dynamical system, the simplified system of equations can often be analyzed more readily than the full set of equations. For example, the phase plane portrait of the reduced model in Fig 2i provides real insight into the properties of the full, six-variable model of a bistable switch. Furthermore, both sQSSA and tQSSA models have fewer kinetic parameters than the full model, which allows us to sidestep what might be serious parameter identifiability problems [26,27]; in particular, the forward and reverse binding constants (k f and k b ) of the enzyme-substrate complex disappear in the reduced models. Hence, it is not necessary to measure experimentally the rate constants of these (presumably) fast reactions, only their ratio, k f , matters, and K M can often be estimated experimentally (e.g., [80]).
Because the simple systems that we have considered in this review are often used as modules in more complex reaction mechanisms, the tQSSA can be used to simplify these complex systems [10,19,25,[113][114][115][116]. When doing so, as multiple timescales may exist, the validity conditions for the tQSSA should be carefully checked so that the tQSSA is not misused. For instance, the rates of association or dissociation of complexes should be fast compared to other rate processes in the network, such as synthesis, degradation, phosphorylation, and dephosphorylation. Furthermore, unlike the examples considered in this review, the algebraic equations underlying the tQSSA of more complex systems may share variables in ways that preclude explicit solution [10], and the QSS equations need to be solved numerically. Excellent computer algorithms for solving differential-algebraic equations are available in [117].
In circumstances where a tQSSA model may be preferred over an sQSSA model, the problem lies in identifying the proper change from fast variables to slower ones (i.e., rabbits to turtles) that ensures a separation of timescales over a wide range of conditions [53,113]. For instance, in enzyme kinetics, the concentration of free substrate, S, was replaced with the slower variable,Ŝ (Fig 1a), and in the GK mechanism, S and S P were replaced with the slower variablesŜ andŜ P (Fig 2a). For complicated reaction networks, the appropriate change of variables can be hard to see [113]. Thus, it would be nice to have an algorithm that identifies an optimal change of variables and performs symbolic or numerical calculation of QSS equations. Such facilities would greatly improve the reliability of dynamical models in computational cell biology.
Taking all these factors into consideration, we recommend to computational biologists who are modeling protein interaction networks with disparate timescales that they simplify their models using the tQSSA rather than the sQSSA. In situations where binding partners (e.g., enzyme and substrate, receptor and ligand, and transcriptional regulators) are in comparable concentrations, tQSSA models are more reliable in predicting the dynamical properties of deterministic interactions and in capturing the fluctuations of stochastic interactions.