Generalized Michaelis–Menten rate law with time-varying molecular concentrations

The Michaelis–Menten (MM) rate law has been the dominant paradigm of modeling biochemical rate processes for over a century with applications in biochemistry, biophysics, cell biology, systems biology, and chemical engineering. The MM rate law and its remedied form stand on the assumption that the concentration of the complex of interacting molecules, at each moment, approaches an equilibrium (quasi-steady state) much faster than the molecular concentrations change. Yet, this assumption is not always justified. Here, we relax this quasi-steady state requirement and propose the generalized MM rate law for the interactions of molecules with active concentration changes over time. Our approach for time-varying molecular concentrations, termed the effective time-delay scheme (ETS), is based on rigorously estimated time-delay effects in molecular complex formation. With particularly marked improvements in protein–protein and protein–DNA interaction modeling, the ETS provides an analytical framework to interpret and predict rich transient or rhythmic dynamics (such as autogenously-regulated cellular adaptation and circadian protein turnover), which goes beyond the quasi-steady state assumption.


Introduction
Since proposed by Henri [1] and Michaelis and Menten [2], the Michaelis-Menten (MM) rate law has been the dominant framework for modeling the rates of enzyme-catalyzed reactions for over a century [1][2][3][4].The MM rate law has also been widely adopted for describing other bimolecular interactions, such as reversible binding between proteins [5][6][7], between a gene and a transcription factor [8,9], and between a receptor and a ligand [10,11].The MM rate law hence serves as a common mathematical tool in both basic and applied fields, including biochemistry, biophysics, pharmacology, systems biology, and many subfields of chemical engineering [12].The derivation of the MM rate law from the underlying biochemical mechanism is based on the steady-state approximation by Briggs and Haldane [3], referred to as the standard quasi-steady state approximation (sQSSA) [12][13][14][15][16][17].The sQSSA, however, is only valid when the enzyme concentration is low enough and thus the concentration of enzyme-substrate complex is negligible compared to substrate concentration [14].This condition may be acceptable for many metabolic reactions with substrate concentrations that are typically far higher than the enzyme concentrations.
Nevertheless, in the case of protein-protein interactions in various cellular activities, the interacting proteins as the "enzymes" and "substrates" often show the concentrations comparable with each other [18][19][20].Therefore, the use of the MM rate law for describing protein-protein interactions has been challenged in its rationale, with the modified alternative formula from the total quasi-steady state approximation (tQSSA) [12,13,[21][22][23][24][25][26][27].
The tQSSA-based form is generally more accurate than the MM rate law from the sQSSA, for a broad range of combined molecular concentrations and thus for protein-protein interactions as well [12,13,[21][22][23][24][25][26][27].The superiority of the tQSSA has not only been proven in the quantitative, but also in the qualitative outcomes of systems, which the sQSSA sometimes fails to predict [12,21].Later, we will provide the overview of the tQSSA and its relationship with the conventional MM rate law from the sQSSA.
Despite the correction of the MM rate law by the tQSSA, both the tQSSA and sQSSA still rely on the assumption that the concentration of the complex of interacting molecules, at each moment, approaches an equilibrium (quasi-steady state) much faster than the molecular concentrations change [12,14,24].Although this quasi-steady state assumption may work for a range of biochemical systems, the exact extent of such systems to follow that assumption is not clear.Numerous cellular processes do exhibit active molecular concentration changes over time, such as in signal responses, circadian oscillations, and cell cycles [6,7,21,[28][29][30][31], calling for a better approach to even cover the time-varying molecular concentrations that may not strictly adhere to the quasi-steady state assumption.
In this study, we report the generalization of the MM rate law, whereby the interaction of time-varying molecular components is more properly described than by the tQSSA and sQSSA.This generalization is the correction of the tQSSA with rigorously estimated, timedelay effects affected by free molecule availability.Our formulation, termed the effective time-delay scheme (ETS), well accounts for the transient or oscillatory dynamics and experimental data patterns of biochemical systems with the relevant analytical insights, which are not captured by the previous methods.Surprisingly, we reveal that the existing quasi-steady state assumption can even fail for extremely slow changes in protein concentrations under autogenous regulation, whereas the ETS does not.In addition, the ETS allows the natural explanation of rhythmic degradation of circadian proteins without requiring explicitly-rhythmic post-translational mechanisms; this is not straightforward within the quasi-steady state assumption.As an added feat, the ETS improves kinetic parameter estimation.As demonstrated in a number of contexts such as autogenouslyregulated cellular adaptation and circadian oscillations, our approach offers a useful theoretical framework to interpret and predict rich transient or rhythmic dynamics of biochemical systems with a wide range of applicability.

Theory development
First, we present the outline of the tQSSA, sQSSA, and our generalized MM rate law.
Consider two different molecules A and B that bind to each other and form complex AB, as illustrated in Fig. 1(a).For example, A and B may represent two participant proteins in heterodimer formation, a chemical substrate and an enzyme in a metabolic reaction, and a solute and a transporter in membrane transport.The concentration of the complex AB at time , denoted by   , changes over time as in the following equation from the massaction law: .
(1) Here,   and   denote the total concentrations of A and B, respectively, and hence

𝐴 𝑡
and     are the concentrations of free A and B. The temporal profiles of   and   are allowed to be very generic, e.g., even with their own feedback effects as addressed later. denotes the association rate of free A and B.  is the effective "decay" rate of AB with  ≡     where  ,  , and  stand for the dissociation, translocation, and dilution rates of AB, respectively, and  for the chemical conversion or translocation rate of A or B upon the formation of AB.In other words, for the sake of generality,  is not limited to a dissociation event but encompasses all rate events to lower the level of AB [Fig.1(a)].
In the tQSSA, the assumption is that   approaches the quasi-steady state fast enough each time, given the values of   and   [12,24].This assumption and the notation  ≡  / lead Eq. ( 1) to an estimate     with the following form (Text S1): (3) Although the tQSSA looks a little complex, it only involves a single parameter  and is easy to implement in a computer program.As mentioned earlier, the tQSSA is generally more accurate than the conventional MM rate law [12,13,[21][22][23][24][25][26][27].To obtain the MM rate law, consider a rather specific condition,   ≪    or   ≪    .(4) In this condition, the Padé approximant for   takes the following form:   .
In other words, Eq. ( 5) would be valid when the concentration of AB complex is negligible compared to either A or B's concentration.This condition is essentially identical to the assumption in the sQSSA resulting in the MM rate law [14].In the example of a typical metabolic reaction with   ≪   for substrate A and enzyme B, Eq. ( 4) is automatically satisfied and Eq. ( 5) further reduces to the familiar MM rate law       /    , i.e., the outcome of the sQSSA [1][2][3][4][12][13][14].To be precise, the sQSSA uses the concentration of free A instead of   , but we refer to this formula with   as the sQSSA because the complex is assumed to be negligible in that scheme.Clearly,  here is the Michaelis constant, commonly known as  .
The application of the MM rate law beyond the condition in Eq. ( 4) invites a risk of erroneous modeling results, whereas the tQSSA is relatively free of such errors and has wider applicability [12,13,[21][22][23][24][25][26][27].Still, both the tQSSA and sQSSA stand on the assumption that   approaches the quasi-steady state fast enough each time before the marked temporal change of   or   .We now relax this quasi-steady state assumption and generalize the approximation of   to the case of time-varying   and   , as the main objective of this study.and ( 6)- (8).Simulations in (a) and (b) are based on periodic oscillation models in Texts S7 and S8, respectively, with their parameters in Table S7.
Suppose that   may not necessarily approach the quasi-steady state each time but stays within some distance from it.As detailed in Text S1, we linearize the right-hand side of Eq.
(   and   is notably sustained in the complex formation, as shown in Text S1.We will refer to this formulation as the effective time-delay scheme (ETS), and its relationship with the tQSSA is depicted in Fig. 1(a).
We propose the ETS as the generalization of the MM rate law for time-varying molecular concentrations that may not strictly adhere to the quasi-steady state assumption.If the relaxation time in complex formation is so short that the effective time delay in Eq. ( 6) can be ignored, the ETS returns to the tQSSA in its form.Surprisingly, we proved that, unlike the ETS, any simpler new rate law without a time-delay term would not properly work for active concentration changes over time (Text S3).Nevertheless, one may question the analytical utility of the ETS, regarding the apparent complexity of its time-delay term.In the examples of autogenously-regulated cellular adaptation and rhythmic protein turnover below, we will use the ETS to deliver valuable analytical insights into the systems whose dynamics is otherwise ill-explained by the conventional approaches.
About the physical interpretation of the ETS, we notice that the effective time delay is inversely linked to free molecule availability, as  ∆   1      2  from Eq. ( 2).Here,     2          , which is the total free molecule concentration at the quasi-steady state each time.In other words, the less the free molecules, the more the time delay, which is at most  .One can understand this observation as follows:    in Eq. ( 1) gives the expectation that the decay time-scale ( ) of the complex may approximate the relaxation time.Yet, the relaxation time is shorter than  , because free A and B are getting depleted over time as a result of their complex formation and therefore the complex formation rate          in Eq. ( 1) continues to decline towards quicker relaxation of the complex level.This free-molecule depletion effect to shorten the relaxation time is roughly proportional to the free molecule concentration itself (Text S1).Hence, the relaxation time takes a decreasing function of the free molecule concentration, consistent with the above observation.Clearly, the free molecule concentration would be low for relatively few A and B molecules with comparable concentrations-i.e., small     and     in Eq. ( 3).In this case, the relaxation time would be relatively long and the ETS shall be deployed instead of the tQSSA or sQSSA.We thus expect that proteinprotein interactions would often be the cases in need of the ETS compared to metabolic reactions with much excess substrates not binding to enzymes, as will be shown later.
Thus far, we have implicitly assumed the continuous nature of molecular concentrations as in Eq. ( 1).However, there exist biomolecular events that fundamentally deviate from this assumption.For example, a transcription factor (TF) binds to a DNA molecule in the nucleus to regulate mRNA expression and the number of such a TF-DNA assembly would be either 1 or 0 for a DNA site that can afford at most one copy of the TF [Fig.1(b)].This inherently discrete and stochastic nature of the TF-DNA assembly is seemingly contrasted with the continuous and deterministic nature of the molecular complex in Eq. ( 1).To rigorously describe this TF-DNA binding dynamics, we harness the chemical master equation [32] and introduce quantities   and   , which are the total TF concentration and the TF-DNA assembly concentration averaged over the cell population, respectively (Text S4).
According to our calculation, the quasi-steady state assumption leads to the following approximant for   :   ≡ , (7) where  ≡  / with  and  as the TF-DNA binding and unbinding rates, respectively, and  is the nuclear volume (Text S4).In fact,   in Eq. ( 7) corresponds to a special case of the previously-studied, stochastic quasi-steady state approximation (stochastic QSSA) [33,34] for arbitrary molecular copy numbers such as for multiple DNA binding sites.Of note, the stochastic QSSA becomes close to the tQSSA as its deterministic version, if ∆  ≫ 1 [34].
in Eq. ( 7) looks very similar to the MM rate law, considering the "concentration" of the DNA site ( ).Nevertheless,   is not a mere continuum of Eq. ( 5), because the denominator in   includes    , but not     .In fact, the discrepancy between   and Eq. ( 5) comes from the inherent stochasticity in the TF-DNA assembly (Text S4).In this regard, directly relevant to   is the stochastic version of the MM rate law with denominator       proposed by Levine and Hwa [35], because the DNA concentration   is  in our case.  is a fundamentally more correct approximant for the DNA-binding TF level than both the tQSSA and sQSSA in Eqs. ( 2) and ( 5).Therefore, we will just refer to   as the QSSA for the TF-DNA interactions.
Still, the use of   stands on the quasi-steady state assumption.We relax this assumption and generalize the approximation of   to the case of time-varying TF concentration.In a similar way to obtain   in Eq. ( 6), we propose the following approximant for   (Text S4): ) This formula represents the TF-DNA version of the ETS, and its relationship with the QSSA is illustrated in Fig. 1(b).The time-delay term in Eq. ( 8) has a similar physical interpretation to that in Eq. ( 6).Besides, this term is directly proportional to the probability of the DNA unoccupancy at the quasi-steady state, according to Eq. (7).
Through numerical simulations of various systems with empirical data analyses, we found that the ETS provides the reasonably accurate description of the deviations of time-course molecular profiles from the quasi-steady states (Texts S6-S8).This result was particularly evident for the cases of protein-protein and TF-DNA interactions with time-varying protein concentrations.In these cases, the ETS unveils the importance of the relaxation time (effective time delay) in complex formation to the shaping of molecular profiles, otherwise difficult to clarify.Yet, the use of the sQSSA or tQSSA is practically enough for typical metabolic reaction and transport systems, without the need for the ETS.The strict mathematical conditions for the validity of the ETS as well as those of the quasi-steady state assumption are derived in Text S5.

Autogenous control
Adaptation to changing environments is a process of biological control.The ETS offers an analytical tool for understanding transient dynamics of such adaptation processes, exemplified by autogenously regulated systems where TFs regulate their own transcription.
This autogenous control underlies cellular responses to various internal and external stimuli [36,37].We here explore the case of positive autoregulation and show that the quasi-steady state assumption does not even work for extremely-slow protein changes near a tipping point.The case of negative autoregulation is covered in Text S10.
In the case of positive autoregulation, consider a scenario in Fig. 2(a) that proteins enhance their own transcription after homodimer formation and this dimer-promoter interaction is facilitated by inducer molecules.The inherent cooperativity from the dimerization is known to give a sigmoid TF-DNA binding curve, resulting in abrupt and history-dependent transition events [36,38].We here built the full kinetic model of the system without the ETS, tQSSA, or other approximations of the dimerization and dimer-promoter interaction (Text S9).As the simulated inducer level increases, Fig. 2(b) demonstrates that an initially low, steady-state protein level undergoes an abrupt leap at some point  , known as a transition or tipping point.This discontinuous transition with only a slight inducer increase signifies a qualitative change in the protein expression state.Reducing the inducer level just back to the transition point  does not reverse the protein state, which is sustained until more reduction in the inducer level [Fig.2(b)].This history-dependent behavior, hysteresis, indicates the coexistence of two different stable states of the protein level (bistability) between the forward and backward transitions [36,38].
Other than steady states, we examine how fast the system responds to signals.Upon acute induction from a zero to certain inducer level (  ), the protein level grows over time towards its new steady state and this response becomes rapider at stronger induction away from the transition point [Fig.2(c)].Conversely, as the inducer level decreases towards the transition point, the response time continues to increase and eventually becomes diverging (in this study, response time is defined as the time taken for a protein level to reach 90% of its steady state).This phenomenon has been called "critical slowing down" [39][40][41].
Regarding this near-transition much slow protein growth, one may expect that the quasi- where   and   represent the concentrations of the unmodified and -th modified proteins, respectively ( 1, 2, ⋯ ,  and  is the total number of the PTMs with  1),   is the protein production rate through mRNA expression and translation,  denotes the protein's ( 1)-th modification rate ( 0, 1, ⋯ ,  1), and  denotes the turnover rate of the -th modified protein.
Next, we apply the ETS to the PTM processes in the model for the analytical estimation of the protein degradation rate.We observed the mathematical equivalence of the PTM processes and the above-discussed TF-DNA interactions, despite their different biological contexts (Text S11).This observation leads to the estimate   of the protein degradation rate as   ≡ min   ,   1 ⋯ , (12) where   is the total protein concentration,  and  are the rates of the two slowest PTM and turnover steps in the protein degradation pathway (the step of  precedes that of  ; see Text S11), and the rightmost formula is to simplify   with the Taylor expansion.
The use of   may not satisfactorily work for the degradation depending on many preceding PTMs, but still helps to capture the core feature of the dynamics.
Strikingly, the quasi-steady state assumption does not predict a rhythmic degradation rate, as the QSSA version of Eq. ( 12) gives rise to a constant degradation rate,     ⁄ (Text S11).In contrast, the ETS naturally accounts for the degradation rhythmicity through the effective time delay in the degradation pathway.The rightmost formula in Eq. (12) indicates that the degradation rate would be an approximately increasing function of   /  and thus increase as time goes from the ascending to descending phase of the protein profile.This predicted tendency well matches the experimental data patterns in Figs.
3(a) and 3(b).Fundamentally, this degradation rhythmicity roots in the unsynchronized interplay between protein translation, modification, and turnover events [42].For example, in the case of protein ubiquitination, ubiquitin ligases with a finite binding affinity would not always capture all newly-translated substrates, and therefore a lower proportion of the substrates can be ubiquitinated during the ascending phase of the substrate profile than during the descending phase.The degradation rate partially follows this ubiquitination pattern.Additional PTMs like phosphorylation, if required for the ubiquitination, can further retard the full substrate modification and thereby increase the degradation rhythmicity for a given substrate profile.One may expect that these effects would be enhanced with more limited ubiquitin ligases or kinases, under the condition when the substrate level shows a strong oscillation.This expectation is supported by the relative amplitude of the degradation rate estimated by Eq. ( 12): where 〈•〉 denotes a time average.Here, the relative amplitude of the degradation rate is proportional to 1   ⁄ as well as to the amplitude of   /  .Therefore, limited ubiquitin ligases or kinases, and strong substrate oscillations increase the rhythmicity of the degradation rate.Given a substrate profile, multiple PTMs can further enhance this degradation rhythmicity because they invite the possibility of smaller  and  values than expected for the case of only a single PTM.Moreover, Eq. ( 12) predicts that the degradation rate would peak around the peak time of   /  .
In the example of Fig.   [44,45,59].Horizontal white and black segments correspond to light and dark intervals, respectively.(c) A simulated protein degradation rate from the full kinetic model and its ETS-and QSSA-based estimates, when the degradation depends on a single PTM.In addition, the protein abundance profile is presented here (gray solid line).A vertical dashed line corresponds to the peak time of   /  where   is a protein abundance.The parameters are provided in Table S9.(d) The probability distribution of the peak-time difference between a degradation rate and   /  for each number of PTMs () required for the degradation.The probability distribution was obtained with randomly-sampled parameter sets in Table S6.(e) The probability distribution of the relative amplitude of a simulated degradation rate (top) or its estimate in Eq. ( 13) (bottom) for each , when the relative amplitude of a protein abundance is 1. (f) The probability distribution of the ratio of the simulated to estimated relative amplitude of a degradation rate for each .For more details of (a)-(f), refer to Text S11.Together, the ETS provides a useful theoretical framework of rhythmic degradation of circadian proteins, which is hardly explained by the quasi-steady state assumption.

Parameter estimation
The use of an accurate function of variables and parameters is important for good parameter estimation by the fitting of the parameters [13,48,49].Parameter estimation is a crucial part of pharmacokinetic-pharmacodynamic (PK-PD) analysis for drug development and clinical study design [50,51].Yet, the MM rate law is widely deployed for PK-PD models integrated into popular simulation and statistical analysis tools.To raise a caution against the unconditional use of the quasi-steady state assumption in parameter estimation, we here compare the accuracies of the tQSSA-and ETS-based parameter estimates.Because the sQSSA-based parameter estimates have already been known as less accurate than the tQSSA-based ones [12,49] and our own analysis supports this claim (Fig. S8), we will henceforth skip the use of the sQSSA.
Specifically, we consider a protein-protein interaction model with time-varying protein concentrations (Text S12).To the "true" profile of the protein complex [i.e.,   in Eq. ( 1)], we fit the ETS [  in Eq. ( 6)] or the tQSSA [  in Eq. ( 2

Discussion
The quasi-steady state assumption involves the approximation by time-scale separation where the "fast" components of a system undergo instantaneous equilibrium and only the "slow" components govern the relevant dynamics.The time-scale separation has been a long practice in many different areas, such as the Monod-Wyman-Changeux model of allosteric effects, the Ackers-Johnson-Shea model of gene regulation by λ phage repressor, and the Born-Oppenheimer approximation in quantum chemistry [53][54][55].If some prediction from the time-scale separation deviates from empirical data, our study may provide a useful intuition about this deviation based on an overlooked time-delay effect in that system.
We here proposed the ETS as a theoretical framework of molecular interaction kinetics with time-varying molecular concentrations.The utility of the ETS for transient or oscillatory dynamics originates in the rigorous estimation of the relaxation time in complex formation, i.e., the effective time delay.In the cases of protein-protein and TF-DNA interactions, the ETS manifests the importance of the effective time delay for the time-course molecular profiles distinct from the quasi-steady states.Accordingly, the ETS provides valuable analytical insights into the signal response time under autogenous regulation and the spontaneous establishment of the rhythmic degradation rates of circadian proteins.In addition, the ETS improves kinetic parameter estimation with a caution against the unconditional use of the quasi-steady state assumption.Our approach enhances the mathematical understanding of the time-varying behaviors of complex-complete massaction models [38,42,56] beyond only their steady states.
Further elaboration and physical interpretation of our framework, in concert with extensive experimental profiling of molecular complexes in regulatory or signaling pathways [18,19], are warranted for the correct explanation of the interplay of cellular components and its functional consequences.Although the simulation and empirical data presented here are supportive of the ETS, direct experimental validation is clearly warranted.This validation could involve the measurements of the time-series of molecular complex concentrations by co-immunoprecipitation assays or other techniques.High temporal resolution data are preferred for their comparison with the ETS-based profiles.Lastly, comprehensive consideration of stochastic fluctuations in molecular binding events [32,57,58] beyond the TF-DNA interactions in this study would be a fruitful endeavor for more complete development of our theory, through possible extension of the existing stochastic QSSA [33,34].

Fig. 1 .
Fig. 1.Generalization of the MM rate law for time-varying molecular concentrations, referred to as the ETS.(a) Two different molecules A and B bind to each other and form their complex.(b) A TF binds to a DNA molecule to regulate mRNA expression (RNA polymerase and other molecules are omitted here).In (a) and (b), the graphs show the comparison among the exact time-course profile of the complex concentration, the tQSSA-based (a) or QSSA-based (b) profile, and the ETS-based profile.The relationship between the tQSSA (or QSSA) and the ETS is illustrated through the effective time delay in the ETS.Notations  ,  ,  , , ∆  , , and   are defined in the description of Eqs.(1)-(3) steady state assumption would work properly near that transition point.To test this possibility, we modified the full model by the tQSSA and QSSA of the dimerization and dimerpromoter interaction, respectively, and call this modified model the QSSA-based model.For comparison, we created another version of the model by the ETS of the dimerization and dimer-promoter interaction and call this version the ETS-based model (Text S9).Across physiologically-relevant parameter conditions, we compared the QSSA-and ETS-based model simulation results to the full model's (Text S12 and Table S5).Surprisingly, the QSSA-based model often severely underestimated the response time, particularly near a transition point, while the ETS-based response time was relatively close to that from the full model [ 10 and Text S12; e.g., 8.5-hour shorter and 0.5-hour longer response times in Fig. 2(c) (left) in the QSSA and ETS cases, respectively].This unexpected mismatch between the QSSA and full model results comes from the following factors: because the QSSA model discards the effective time delay in dimerization and dimer-promoter interaction, this model accelerates positive feedback, transcription, and protein production, and thus shortens the response time.Near the transition point, although the protein level grows very slowly, a little higher transcription activity in the QSSA model substantially advances the protein growth with near-transition ultrasensitivity that we indicated above.Therefore, the QSSA model shortens the response time even near the transition point.

Fig. 2 .
Fig. 2. Positive autoregulation and induction kinetics.(a) Protein production mechanism with positive autoregulation in the presence of inducers.(b) Bifurcation diagram of the simulated protein level as a function of  (proxy for an inducer level).The steady state is plotted as  increases (solid line) or decreases (dashed line).Acute induction can be simulated by a sudden change of  0 to   in the shaded area.(c) Time-series of protein levels from the full, ETS, and QSSA models upon acute induction at time 0 h with  2.42 (left) or  200 (right).(d) The full model-to-QSSA difference in response time as a function of 1   / ⁄ .Both the simulated and analytically-estimated differences are presented.The analytical estimation is based on Eq. (9).For more details of (b)-(d), refer to Text S9 and TableS8.
3(c)  for a single PTM case, the simulated degradation rate from the aforementioned full kinetic model exhibits the rhythmic profile in excellent agreement with the ETS-predicted profile.Notably, the peak time of the simulated degradation rate is very close to that of   /  as predicted by the ETS.Indeed, the peaks of the degradation rates show only < 1h time differences from the maximum   /  values across most (89-99%) of the simulated conditions of single to triple PTM cases [Fig.3(d);Text S12 and TableS6].In addition, for each substrate profile, the simulated degradation rate tends to become more rhythmic and have a larger relative amplitude as the number of the PTMs increases [Fig.3(e)],supporting the above argument that multiple PTMs can facilitate degradation rhythmicity.The estimated relative amplitude in Eq. (13) also shows this tendency for single to double PTMs, yet not clearly for triple PTMs unlike the simulated relative amplitude [Fig.3(e)].This inaccuracy with the triple PTMs comes from the accumulated errors over multiple PTMs in our estimation, as we indicated early.Still, the estimate in Eq. (13) accounts for at least the order of magnitude of the simulated relative amplitude, as the ratio of the simulated to estimated relative amplitude almost equals 1 for a single PTM case and remains to be  1 for double and triple PTM cases [Fig.3(f)].

Fig. 3 .
Fig. 3. Rhythmic degradation of circadian proteins.(a) The experimental abundance levels (solid line) and degradation rates (open circles) of the mouse PERIOD2 (PER2) protein [43].(b) The experimental )] and estimate the original parameters of the model[52]: the ETS-based fitting can estimate both parameters  and  , and the tQSSA-based fitting can estimate only .Likewise, we consider a TF-DNA interaction model with time-varying TF concentration (Text S12).The ETS-based fitting can estimate both  and  , and the QSSA-based fitting can estimate only .In the case of protein-protein interactions, Fig.4(a) reveals that the ETS tends to improve the parameter estimation over the tQSSA, with more accurately estimated : most of  values (89.4%) estimated by the ETS show smaller relative errors than the tQSSA-based estimates and their 69.3% even show relative errors less than half the tQSSA's.In the case of TF-DNA interactions, the ETS still offers an improvement in the estimation of  [Fig.4(b)]:most of the ETS-estimated  values (90.3%) show smaller relative errors than the QSSAestimated ones and their 51.8% even show relative errors less than half the QSSA's.

Fig. 4 .
Fig. 4. Parameter estimation for protein-protein and TF-DNA interaction models.(a) The scatter plot of the relative errors of the tQSSA-and ETS-estimated  values for a protein-protein interaction model.(b) The scatter plot of the relative errors of the QSSA-and ETS-estimated  values for a TF-DNA interaction model.In (a) and (b), a diagonal line corresponds to the cases where the two estimates have the same relative errors.(c) The probability distribution of the relative error of the ETS-estimated  for the protein-protein interaction model in (a).(d) The probability distribution of the relative error of the ETS-estimated  for the TF-DNA interaction model in (b).Regarding (a)-(d),a subset of simulated conditions gave relative errors outside the presented ranges here, but they did not alter the observed patterns.For more details, refer to Text S12.
1) around     and estimate   's solution as the time integral of   (where the inclusion of an effective time delay  ∆  in complex formation.This delay is the rigorous estimate of the molecular relaxation time during which the effect of instantaneous

Table S8 .
den47] an inducer level and its value at the transition point, respectively,  is the sum of protein degradation and dilution rates, and  and  are parameters inversely proportional to the effective time delays in dimerization and dimer-promoter interaction, respectively.The additional details and the definition of parameter  are provided in Text S9. degradation rates of proteins with circadian production can spontaneously emerge without any explicitly time-dependent regulatory mechanism of the degradation processes[42,47].If the rhythmic degradation rate peaks at the descending phase of the protein profile and stays relatively low elsewhere, it is supposed to save much of the Related to this point, the ETS allows the analytical calculation of response time and its QSSAbased estimate.In this calculation, we considered two different stages of protein growth-its early and late stages (Text S9) and found that the QSSA model underestimates response time mainly at the early stage.This calculation suggests that the exact response time would be difference indefinitely grows as  decreases towards  , as a linear function of 1/   / .This prediction can serve as a testbed for our theory and highlights far excessive elongation of near-transition response time (compared to the QSSA) as an amplified effect of the relaxation time in complex formation.This amplified effect is the result of the near-transition ultrasensitivity that we indicated above.Consistent with our prediction, the full model simulation always shows longer response time than the QSSA model simulation and the difference is linearly scaled to 1/   / as exemplified by Fig.2(d) ( 0.98 in simulated conditions; see Text S12).Moreover, its predicted slope against 1/   / [i.e., 2 (6.28 ⋯) multiplied by    ] is comparable with the simulation results [7.3 ± 0.3 (avg.± s.d. in simulated conditions) multiplied by    ; see Text S12].The agreement of these nontrivial predictions with the numerical simulation results proves the theoretical value of the ETS.Again, we raise a caution against the quasi-steady state assumption, which unexpectedly fails for very slow dynamics with severe underestimation of response time, e.g., by a few tens of hours in the case of Fig. 2(d). of pathophysiological conditions [7,9,21,28-30].According to previous reports, some circadian clock proteins are not only rhythmically produced but also decompose with rhythmic degradation rates [Figs.3(a) and 3(b)] [42-46].Recently, we have suggested that the rhythmic of the protein to a circadian mRNA expression or translation rate.Yet, a protein degradation rate in the model is not based on any explicitly time-dependent regulatory processes, but on constantly-maintained proteolytic mediators such as constant E3 ubiquitin ligases and kinases.In realistic situations, the protein turnover may require multiple preceding PTMs, like mono-or multisite phosphorylation and subsequent ubiquitination.Our model covers these cases, as well.The model comprises the following equations: