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-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 tim<span class="Chemical">escale 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-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” quasi-steady state approximation (tQSSA), which is generally accurate for any combination of substrate and enzyme concentrations [14-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.

Results

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 the total enzyme concentration (ET ≡ C+E) is constant, as it should be. Furthermore, because the total concentration of substrate + product (ST ≡ S+C+P) is also constant, Eq 1 describes the dynamics of only two independent variables, say C and P. This two-variable 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 quasi-steady state (QSS) where , 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 : 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, ST ≡ S+C+P (see below).

Fig 1

Comparison of the sQSSA and tQSSA for an enzyme-catalyzed reaction.

(a) 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 (Eq 4). (b–c) Heat maps of the predicted maximal fraction of total enzyme that is in complex with substrate, i.e., C(ST)/ET, predicted with sQSSA (Eq 2) and tQSSA (Eq 4). In these panels, color represents the ratio C(ST)/ET, and gray dashed lines are contours of ET /(KM + ST) = 0.1 and 1.0. The sQSSA (Eq 2) predicts that this fraction is independent of total enzyme concentration, ET (e.g., C(ST)/ET = 0.5 when ST = KM regardless of ET) in panel b. On the other hand, the tQSSA (Eq 4) predicts that, as ET 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 ET /(KM + ST) ≪1), the predictions of the two models are nearly identical. (d–e) The simulated accumulation of P and the estimation of the Michaelis constant, KM, from the progress curve of P. For these calculations, we chose kf = 10, kb = 9, kcat = 1, KM = 1, ST = 9, and ET /(KM + ST) = 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) = ST, C(0) = 0, P(0) = 0. If ET is low (i.e., ET /(KM + ST) = 0.1), the simulation of P and the estimation of KM with the sQSSA model are accurate (d). On the other hand, if ET is high (i.e., ET /(KM + ST) = 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 kcat and KM are adjusted to fit the sQSSA model to the true progress curve, then the estimates of kcat and KM are far from the true values (inset; the estimates of kcat 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 ET (d and e). Inset: gray lines are the true value of KM, and the red and blue circles are the mean KM 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 KM that minimizes the fitting error of the sQSSA and tQSSA models. The true values of parameters are used as 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.

The QSSA (Eq 2) implies that and , i.e., that enzyme molecules are half-saturated with substrate when S = KM, regardless of total enzyme concentration, ET (Fig 1b). This seems unlikely because as ET 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 ET is small enough. Indeed, Segel and Slemrod showed that the QSSA is valid only when ET ≪ ST +KM [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, ST ≡ S+C+P, and S can be replaced with ST−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 (ET ≪ ST +KM), 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-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-21]. For the tQSSA, S is replaced with the total substrate concentration 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 .)

To make this change, we rewrite the mass-action rate expr<span class="Chemical">essions (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 : 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 <span class="Chemical">Tzafriri [16]:

Tzafriri also showed that Thus, if kcat ≪ kb, then εtQ ≪ 1 regardless of total enzyme concentration, ET. Furthermore, Tzafriri also showed that εtQ ≪ 1 if either ST or ET (or both) are much greater than KM. Importantly, because < 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-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 has the reasonable consequence that, as ET 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 : which is valid if C2 ≪ ETST [10,16,17]. According to the Padé approximant, enzyme molecules are half-saturated with substrate (i.e., ), when = KM + ET. That is, the total substrate concentration needed for half-saturation of enzyme increases linearly with total enzyme concentration, ET (Fig 1c). Furthermore, when the sQSSA is valid (i.e., ET ≪ ST + KM), the Padé approximant becomes nearly identical with C(S) derived with sQSSA (Eq 2). This explains nearly identical values of C(S) and when ET ≪ ST + KM (Fig 1b and 1c, below the gray dashed line where ).

By substituting (Eq 4) into the full model, we can d<span class="Chemical">escribe the rate of product accumulation in terms of rather than S, and replace the MM rate law (Eq 3) by where (Fig 1a). Note that the conservation of ST 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 ST + KM, 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., KM = 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, ET ≪ ST + KM, 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-31]. Importantly, in vivo enzyme concentrations can vary extensively inside a cell and are usually much higher than those used experimentally in vitro [11-13,32]. For instance, the enzyme cytochrome P450, which metabolizes many drugs in the liver, is present in concentrations greatly exceeding [drug] + KM [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-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 substrate–product pair (S and SP) 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 SP is dephosphorylated by phosphatase D, which can be described by the following system of ODEs: where concentrations of total substrate (ST ≡ S + SP + ES + DSP), total kinase (ET ≡ E + ES), and total phosphatase (DT ≡ D + DSP) 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).

Fig 2

sQSSA, but not tQSSA, makes false predictions of zero-order ultrasensitivity and bistability.

(a) In the GK mechanism, substrate (S) is phosphorylated by kinase (E), and then phosphorylated substrate (SP) 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 (DSP) with either the sQSSA (Eq 7) or the tQSSA (Eq 8). In the tQSSA model, and are total unphosphorylated substrate and phosphorylated substrate concentration, respectively. (b–d) As the total concentration of kinase (ET) increases, the steady-state concentration of phosphorylated substrate () increases. The sQSSA model predicts a steep sigmoidal response of regardless of the level of DT (in this case, , 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 DT increases. In these calculations, ST = 100, kfe = kfd = 10, kd = ke = 1.7, and KME = KMD = 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 -dependent synthesis of E (kpos) creates a positive feedback loop between and E. (i–j) In panel i, the steady states of the extended positive feedback models (circles) are intersections of the nullclines of ET defined by (solid black lines) and of the original models of the GK mechanism defined by (colored lines; for the sQSSA model as complexes are assumed to be negligible). As the nullcline of ET changes, depending on the basal synthesis rate of E (kbasal), 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 kpos = kdeg = = 0.01. GK, Goldbeter–Koshland; sQSSA, standard quasi-steady state approximation; tQSSA, total quasi-steady state approximation.

In particular, assuming low concentrations of total kinase and phosphatase (ET,DT ≪ ST), Goldbeter and Koshland used the sQSSA for ES and DSP (Eq 2) to describe the rates of phosphorylation and dephosphorylation in terms of S and SP: where and (Fig 2a). Because the sQSSA assumes that ES and DSP are negligible, they can be neglected in the conservation condition ST ≡ S + SP + ES + DSP, and thus S can be replaced with ST − SP 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 DSP (Eq 4) [10,19,25,43,44]. This leads to the tQSSA model of the GK mechanism (Fig 2a): where and . Using the conservation condition can be replaced with in the above equation, leading to a one-variable tQSSA model. Note that the complexes (ES and DSP) are not neglected in the conservation condition unlike the sQSSA model.

As the concentration of kinase (ET) increases, the phosphorylation reaction is more likely to occur, and so the steady-state fraction of phosphorylated substrate (SP/ST) increases. For the sQSSA model (Eq 7), Goldbeter and Koshland derived an explicit function for SP/ST in terms of ET [10,41]. They showed that for KME, KMD ≪ ST, the function is sigmoidally shaped with an abrupt jump from SP/ST ≈ 0 to SP/ST ≈ 1 at keET ≈ kdDT, 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 KME, KMD ≪ ST (i.e., tight binding). A similar sharp transition, close to keET ≈ kdDT, is also seen for the steady-state fraction of total phosphorylated substrate () in the full model (Eq 6) and tQSSA model (Eq 8) when the enzyme concentrations are small compared to total substrate concentration (ET, DT ≪ ST; 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 () act as TFs for E, the kinase that converts S to SP (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 (ET):

The extended model of a positive feedback loop can be simplified using either the sQSSA (Eq 7) or the tQSSA (Eq 8), where for the sQSSA model (as <span class="Gene">DSP is assumed to be negligible) or for the tQSSA model.

The steady states of the sQSSA and tQSSA positive feedback models are found at the intersections of the nullclines of ET (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 (DT/ST = 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 (Fig 2i and 2j, filled red circles). On the other hand, when kbasal 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-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-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 SP 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-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-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-68].

In other cas<span class="Chemical">es, an intermediate molecule is involved in the transcriptional negative feedback loop [54]. For instance, repression of transcription can occur via protein sequestration [69-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.

Fig 3

sQSSA, but not tQSSA, leads to the false suppression of oscillations in a protein sequestration–based transcriptional negative feedback loop.

(a) A repressor protein (R) inhibits its own transcriptional activator (A) by forming a 1:1 stoichiometric complex (RA). The QSS of the free activator, which is not sequestered by the repressor, can be approximated with either the sQSSA or the tQSSA. As the QSSA is derived from reversible binding without any conversion rate, the dissociation constant Kd = kb/kf is used instead of the Michaelis constant KM. (b) In the transcriptional negative feedback loop, where the repressor protein sequesters the activator protein, the transcription of the repressor’s mRNA is assumed to be proportional to the fraction of total activator that is free (i.e., A/AT), which can be approximated by either the sQSSA or the tQSSA (a). M, mRNA; P, encoded protein; R, repressor protein; A, activator; RA, repressor–activator complex; , repressor + repressor–activator complex, AT, activator + repressor–activator complex. (c–e) Heat maps of the relative amplitude of simulated trajectories of M. Oscillations are not possible for the sQSSA model (c). On the other hand, when the validity condition for the sQSSA does not hold (i.e., above the gray dashed line), both the tQSSA (d) and the full (e) models simulate oscillations. In these calculations, AT and αM are varied so that (where is the steady-state value of ) varies between 10 and 105. Furthermore, kf = kb = 10 (i.e., Kd = 1) and all other parameter values = 1. Insets: for a particular choice of parameter values (star mark), trajectories of M are simulated, and all trajectories are normalized by the maximum value of the trajectory simulated by the full model. QSS, quasi-steady state; sQSSA, standard quasi-steady state approximation; tQSSA, total quasi-steady state approximation.

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, , is to be given by either the sQSSA or the tQSSA, as displayed in Fig 3a. for the tQSSA model and for the sQSSA model as the complex (RA) is assumed to be negligible.

In the sQSSA model, free activator concentration is given by , and so the rate of transcription of mRNA is (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 (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 , 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 (AT) and the total repressor (), provided Kd ≪ AT, i.e., tight binding between them (Fig 3d). The full model, based on mass-action 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-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 pr<span class="Chemical">esence of tim<span class="Chemical">escale 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-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-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-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, NX, with the relationship NX = [X]∙Ω, where Ω is the volume of the system. Then, based on the nonelementary propensity functions, stochastic simulations are performed with Gillespie’s algorithm [92-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-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-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 DT/ST ≪ 1 (Fig 4c) and when DT /ST = 1 (Fig 4d). When DT /ST ≪ 1, the GK mechanism exhibits zero-order ultrasenstivity, i.e., the fraction of total phosphorylated substrate changes sharply around ET = DT (Fig 2b). Reflecting this sharp change, the fraction of the total phosphorylated substrate () fluctuates wildly in the stochastic simulations, and the stationary distribution of is broad (Fig 4c). On the other hand, when DT/ST = 1, the fraction of total phosphorylated substrate changes very little around ET = DT (Fig 2d), and thus its stationary distribution is narrow (Fig 4d).

Fig 4

tQSSA provides accurate stochastic simulations.

(a) Catalysis-completion time τ (see inset) of the single-substrate enzyme-catalyzed reaction (Fig 1a) varies due to stochastic fluctuations. The distributions of τ simulated with the full model (3.24 ± 1.4) and the tQSSA model (3.22 ± 1.37) are nearly identical. Here, the same parameter set and initial condition are used as in Fig 1e and Ω = 1. Gillespie’s stochastic simulation algorithm is used [107]. (b) The peak-to-peak period τ (see inset) of simulated oscillatory trajectories for the model of a negative feedback loop (Fig 3b) fluctuates due to stochasticity. The distributions of τ simulated by the full model (5.3 ± 0.8) and tQSSA model (5.3 ± 0.78) are nearly identical. Here, the same parameter set and initial condition are used as in the inset of Fig 3d and 3e and Ω = 5 × 10−4. (c–d) In both the presence (c) and absence (d) of zero-order ultrasensitivity, stochastic simulations of the full and tQSSA models lead to nearly identical stationary distributions of : 0.5 ± 0.23 and 0.5 ± 0.22 (c) and 0.5 ± 0.05 and 0.5 ± 0.05 (d). Here, ET = DT and the values of the other parameters are the same as used in Fig 2b and 2d for (c) and (d), respectively. Initial conditions are Sp(0) = S(0) = ST/2, D(0) = E(0) = DT, and Ω = 1. tQSSA, total quasi-steady state approximation.

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 (ET ≪ ST + KM) 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-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 kcat ≫ kb and <span class="Gene">ET + ST is not in excess over KM [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 (kf and kb) 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, , matters, and KM 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-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 SP were replaced with the slower variables and (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.