Emergence and Enhancement of Ultrasensitivity through Posttranslational Modulation of Protein Stability

Signal amplification in biomolecular networks converts a linear input to a steeply sigmoid output and is central to a number of cellular functions including proliferation, differentiation, homeostasis, adaptation, and biological rhythms. One canonical signal amplifying motif is zero-order ultrasensitivity that is mediated through the posttranslational modification (PTM) cycle of signaling proteins. The functionality of this signaling motif has been examined conventionally by supposing that the total amount of the protein substrates remains constant, as by the classical Koshland–Goldbeter model. However, covalent modification of signaling proteins often results in changes in their stability, which affects the abundance of the protein substrates. Here, we use mathematical models to explore the signal amplification properties in such scenarios and report some novel aspects. Our analyses indicate that PTM-induced protein stabilization brings the enzymes closer to saturation. As a result, ultrasensitivity may emerge or is greatly enhanced, with a steeper sigmoidal response, higher magnitude, and generally longer response time. In cases where PTM destabilizes the protein, ultrasensitivity can be regained through changes in the activities of the involved enzymes or from increased protein synthesis. Importantly, ultrasensitivity is not limited to modified or unmodified protein substrates—when protein turnover is considered, the total free protein substrate can also exhibit ultrasensitivity under several conditions. When full enzymatic reactions are used instead of Michaelis–Menten kinetics for the modeling, the total free protein substrate can even exhibit nonmonotonic dose–response patterns. It is conceivable that cells use inducible protein stabilization as a strategy in the signaling network to boost signal amplification while saving energy by keeping the protein substrate levels low at basal conditions.


Regulation of Protein Stability through Posttranslational Modifications
It has been known for some time that posttranslational modifications (PTMs) are important mechanisms for regulating not only the activity of a protein, but also the abundance of a protein by means of changing its stability. A well-studied example is the DNA damage response. Once the tumor suppressor p53 is phosphorylated by upstream kinases, such as ATM (ataxia telangiectasia mutated), in response to DNA double-strand breaks, its half-life increases dramatically from less than 30 min to over 3 h ( Figure 1A), which causes the accumulation of p53 that can induce target gene expression [1,2]. A second example, in some sense of the opposite nature, occurs in the germinal center response of B lymphocytes. B cell receptor-activated MAPK phosphorylates BCL6 (B-cell lymphoma 6), resulting in accelerated degradation of BCL6 by the ubiquitin/proteasome pathway ( Figure 1B), which helps the B cells exit the germinal center response [3]. Many similar examples have pathway ( Figure 1B), which helps the B cells exit the germinal center response [3]. Many similar examples have been reported where protein stabilization or destabilization drives signaling, including IKK-mediated phosphorylation and degradation of IκB in the inflammatory response, Chk1-mediated phosphorylation and proteasomal degradation of Cdc25A during cell cycle arrest, and stabilization of ΔFosB by casein kinase 2-mediated phosphorylation, which might be responsible for long-term adaptation in the brain [4][5][6]. It is thus conceivable-and even likely-that altering protein stability and/or activity through the same PTM event may be an important controllable mode of dual regulation of cellular signaling networks in general. Expressed differently, if the abundance of a protein substrate can be fine-tuned through changes in protein stability, then these changes can in turn be used by the cell as modulators of both the dynamic and steady-state inputoutput (I/O) behaviors of covalent modification cycles (CMCs), which may or may not alter the activity of the protein.

Ultrasensitivity
Cell signaling networks display "ultrasensitivity" if small changes in input are amplified into much larger percentage changes in output [7,8]. An ultrasensitive I/O relationship is generally sigmoidal in shape and often approximated by a Hill function; the terminology suggests that an ultrasensitive response is steeper than the well-known hyperbolic trend of a Michaelian function [9,10]. Embedded in complex network structures such as feedback and feedforward loops, signal amplification is required for cells and organisms to achieve higher-order functions, including differentiation, proliferation, homeostasis, adaptation, and biological rhythms [11,12]. At least six major ultrasensitive response motifs (URM) have been identified in intracellular molecular networks, namely: (i) positive cooperative binding; (ii) homo-multimerization; (iii) multistep signaling; (iv) molecular titration; (v) zero-order CMCs; and (vi) positive feedback [12][13][14]. Each of these URMs has its own unique mechanism achieving signal amplification.

Ultrasensitivity through Zero-Order Covalent Modification Cycle
The ubiquitous zero-order CMC is particularly interesting, as it can generate nearly switch-like responses. A typical implementation is a modifying/demodifying cycle that is driven by PTMs involving phosphorylation, acetylation, oxidation, methylation, or sumoylation [15]. Specifically, protein activities can be regulated through covalent bonding of moieties to certain amino acid residues, such as phosphate to serine, threonine, and tyrosine in the case of phosphorylation, and an acetyl group to lysine in the case of acetylation. The local electrical charge, possibly accompanied by steric changes introduced by these moieties, can greatly affect the protein molecule's interaction with other large or small molecules, thereby turning on or off the activity of the protein as an enzyme, transcription factor, or signaling molecule. Covalent modifications of proteins often require specific enzymes, such as kinases, acetyltransferases, methyltransferases, and oxidases, as well as counteracting (demodification) enzymes catalyzing the reverse reactions, such as phosphatases, deacetylases, demethylases, and reductases.
Signal amplification through CMCs was first predicted and analyzed with a mathematical model by Goldbeter and Koshland Jr. in the early 1980s [16,17]. It occurs when the two opposing enzymes driving the modification cycle of a protein are operating near saturation. In a phosphorylation-dephosphorylation cycle, for example, zero-order ultrasensitivity arises when the amount of protein substrate is at a concentration high enough to saturate the available kinase and phosphatase. Here, the terminology "protein substrate" is used to distinguish this protein from the involved enzymes. Under these conditions, small changes in the amount or activity of either the kinase or phosphatase can dramatically change the steady-state fraction of the amounts of phosphorylated or dephosphorylated substrates. Since the theoretical predictions by the Goldbeter-Koshland model, zero-order ultrasensitivity via covalent modification has been reported in numerous biological settings, in both prokaryotes and eukaryotes [18][19][20][21][22][23].

Caveat of the Goldbeter-Koshland Model Suggests a Mechanism of Signaling Control
One important conceptual simplification of the original Goldbeter-Koshland model is that the total abundance of the protein substrate in the CMC is regarded as constant, which ignores turnover via de novo protein synthesis and degradation. This omission is possibly critical in the context of protein signaling, as proteins are constantly synthesized and degraded. The assumption of constancy may largely be valid when the signaling events driven by PTM occur rapidly in comparison to the protein substrate turnover. However, even if signaling is fast, it is possible-and indeed a frequent observation as mentioned before-that the PTM alters the stability of the protein substrate, which secondarily affects the total amount of the protein substrate. We first reported that, due to the "leakiness" caused by protein turnover, zero-order ultrasensitivity is compromised when turnover is present, and that the steepness of the sigmoidal response deteriorates as the overall protein turnover rate increases [24]. More recently, Mallela et al. further elaborated on the importance of protein synthesis and turnover in affecting zero-order ultrasensitivity of CMCs, especially in the context of multiple PTM cascades sharing the same E3 ligase responsible for protein degradation [25]. Thus, the formerly simple results described by Goldbeter and Koshland are in truth more complicated, as they depend on the kinetic features of the involved enzymes, their saturation, and the degree of protein synthesis and turnover.
Here, we pursue the question how cells may use alternate PTM-induced changes in protein stability as an additional layer of control to modulate the zero-order ultrasensitive response of a CMC. In particular, we ask whether such modulations are sufficient to render or enhance ultrasensitivity by stabilizing the protein substrate, or diminish or destroy it by destabilizing the protein substrate. To answer these questions, we systematically study the governing kinetic features of the protein cycle one by one, with mathematical modeling, which allows us to modify many aspects or combinations of aspects of a protein signaling cycle with full knowledge of the system features and behaviors. We demonstrate that ultrasensitivity can be gained, enhanced, or attenuated for the modified, unmodified, and total protein substrates depending on the conditions affecting stability.

Model Structure and Parameterization
Our goal is to explore how the behavior of a CMC is affected if protein turnover, protein stability, and kinetic features of the governing enzymes are explicitly taken into account. For this exploration, we consider the generic signaling motif of a protein phosphorylationdephosphorylation cycle ( Figure 1C) as an "order-of-magnitude" model, i.e., a numerical model without absolutely precise determination of parameter values and with an expectation of qualitative, rather than quantitative results.
We started with a simple model consisting of two ordinary differential equations (ODEs), formulated in the tradition of mass-action and Michaelis-Menten (MM) kinetics: In the following, this model is referred to as the MM model. R is the protein substrate that is newly synthesized with rate k 0 . It can either be phosphorylated into R p by a kinase X which, as a default, follows typical MM kinetics with Michaelis constant K m1 and maximal velocity V max1 = k 1 X, or it can be degraded with a first-order rate constant k 3 . Analogously, R p can be dephosphorylated by a phosphatase Y that follows MM kinetics with a Michaelis constant K m2 and maximal velocity V max2 = k 2 Y. R p can also be degraded, in this case with a first-order rate constant k 4 .
Default parameter values are presented in Table 1. Since covalent protein modifications such as phosphorylation and dephosphorylation occur rapidly, at the order of seconds to minutes, while protein degradation occurs at a much slower rate, often with halflives at the order of hours, the time scales between these two types of processes are often clearly separated by two or more orders of magnitude. Accordingly, we set default values for k 3 and k 4 to be 1/100 of k 1 /K m1 and k 2 /K m2 , respectively, because these two ratios approximate the apparent first-order time constants at which phosphorylation and dephosphorylation occur when the kinase and phosphatase are far from saturation. Unless otherwise specified, Y is kept constant with value 1. It is well known that the MM kinetics is a simplification that approximates the full mass-action enzymatic kinetics with certain constrains [26]. These constraints include that the substrate exists in excess of the enzyme and the substrate-enzyme complex quickly reaches a quasi-steady state. In many cellular circumstances, some of the assumptions may not be valid. An example is the case where the substrate and enzyme are proteins with comparable concentrations [27,28]. Under conditions where the assumptions are not satisfied, the full kinetic model should be utilized. Indeed, it has been shown that implementing the full model can generate nonlinear dynamics that may be missed by simple MM kinetics [26,29]. To assess the role of these differences, we used the full kinetic model ( Figure 1D and Table 2) to validate or refute the findings based on the approximating MM model. The ODEs for the full model are provided in the Supplemental Material. In this full model, we assume that R or R p , complexed with their respective enzymes X or Y, can still be independently degraded like the respective free forms and with the same half-lives, and when they are destroyed, X and Y are released and recycled. That a protein component in a multimeric protein complex can be specifically targeted and degraded by the ubiquitination/proteasomal pathway has been well established [30][31][32]. This matter is discussed further in the Discussion section. Table 2. Default parameter values for the full model in Figure 1D.

Metrics of Ultrasensitivity
Simulations start with the initial values of variables equal to the respective steady-state concentrations, which are achieved when X = 0 for the MM model or X tot = 0 for the full model, where the initial values of R p , RX, and R p Y are always zero. For fair comparisons, all dose-response (DR) curves are obtained once the model has achieved steady state. The increment size of the input level of X is 1% within the dose range indicated. The degree of ultrasensitivity of a steady-state DR curve can be evaluated with two related metrics that are commonly used in the field. First, the Hill coefficient, n H , may be approximated from the equation where X 0.9 and X 0.1 are the concentrations of X that produce 90% and 10%, respectively, of the maximal response (after subtracting the background response level when X = 0) [12]. n H represents the overall steepness or global degree of ultrasensitivity of the DR curve. Second, we evaluate the local response coefficient (LRC) of a DR curve by calculating all slopes of the curve on dual-log scales, which are equivalent to the ratios of the fractional change in response (R) to the fractional change in dose (D) [7]: The maximal |LRC| of a DR curve (|LRC| max ) represents the maximal amplification capacity of the signaling motif. The comparison between n H and LRC is important as these quantities are not necessarily equivalent and they vary depending on the basal response level and the shape of the DR curve. Typical ultrasensitive responses have |n H | and |LRC| max values substantially above 1. However, if a DR curve is very sigmoidal with high |n H |, but its |LRC| max < 1, then there is no signal amplification, and higherorder cellular functions such as bistability, robust adaptation, and oscillation cannot be achieved [12,33,34]. Thus, n H alone can misrepresent the actual degree of signal amplification and must be cross-examined with LRC.

Simulation Tools
The models were coded in MATLAB (MathWorks, Natick, MA, USA) and all simulations were run with differential equation solver ode15s. The model codes in MATLAB format as well as in R format are available as additional Supplemental files and can also be accessed at the GitHub repository https://github.com/pulsatility/Effects-of-Protein-Stability-on-Ultrasensitivity.git, accessed on: 17 November 2021.

Ultrasensitivity in the Absence of PTM-Induced Changes in Protein Stability
To create a baseline, we start with the default setting k 3 = k 4 = 0.01, which reflects that the phosphorylation status of R does not affect its stability. As a consequence, the total steady-state protein substrate concentration R tot (=R + R p ) remains constant even if the activity of the kinase X varies. Additionally, the k 3 and k 4 values are very small in comparison to k 1 and k 2 . Since R tot typically exceeds K m1 and K m2 by 10-fold or more, and as expected for the CMC motif, the steady-state DR curves of R vs. X and R p vs. X are sigmoidal on the linear scale ( Figure 2A) with n H at −3.51 and 3.51, respectively ( Figure 2B); the negative sign for n H indicates a decreasing or inhibitory response. On a log scale, the quasi-exponential rise in R p and decay of R flatten toward straight lines ( Figure 2C).
The degree of local ultrasensitivity, as measured by LRC, varies across the range of X and peaks in the center of the DR curves at approximately −3.0 and +3.1 for R and R p , respectively ( Figure 2B). Thus, |n H | in this case is a slight overestimate of the corresponding |LRC| max . The phosphorylation and dephosphorylation fluxes (with rates k 1 and k 2 , respectively) are dominant over the relatively small protein turnover fluxes with rates k 3 and k 4 at the steady state for large input values of X ( Figure 2D). In a logarithmic representation, these MM fluxes increase essentially linearly as X increases before approaching plateaus ( Figure 2D). When protein production and degradation are considered negligible, by setting k 0 , k 3 and k 4 to zero, the ultrasensitive responses are slightly enhanced, and the Hill coefficients and (|LRC| max rise in magnitude to 3.74 and 3.45 for R p and −3.74 and −3.34 for R (simulation results not shown). The degree of local ultrasensitivity, as measured by LRC, varies across the range of X and peaks in the center of the DR curves at approximately −3.0 and +3.1 for R and Rp, respectively ( Figure 2B). Thus, |nH| in this case is a slight overestimate of the corresponding |LRC|max. The phosphorylation and dephosphorylation fluxes (with rates k1 and k2, respectively) are dominant over the relatively small protein turnover fluxes with rates k3 and k4 at the steady state for large input values of X ( Figure 2D). In a logarithmic representation, these MM fluxes increase essentially linearly as X increases before approaching plateaus ( Figure 2D). When protein production and degradation are considered negligible, by setting k0, k3 and k4 to zero, the ultrasensitive responses are slightly enhanced, and the Hill coefficients and (|LRC|max rise in magnitude to 3.74 and 3.45 for Rp and −3.74 and −3.34 for R (simulation results not shown).

Effects of Protein Stability on Ultrasensitivity
In this section, we suppose that changes in the stability of Rp can be introduced by the PTM, and thus by means of the kinase X, as it has been observed numerous times [1,[3][4][5][6]. When the protein substrate R is phosphorylated to Rp, the stability of Rp is affected, which translates into an increasing or decreasing rate of degradation, k4. In particular, if the PTM stabilizes Rp, i.e., k4 decreases, the amount of Rp increases, and Rtot is expected to increase accordingly. This rise in Rtot secondarily alters the degree of saturation of the phosphorylation and dephosphorylation reactions and, consequently, is expected to affect the degree of ultrasensitivity in the DR curves. These overall effects could theoretically also be caused by changes in k3, but we focus on k4 because the PTM directly affects the stability of Rp, whereas R is affected only in a secondary manner.

Effects of Protein Stability on Ultrasensitivity
In this section, we suppose that changes in the stability of R p can be introduced by the PTM, and thus by means of the kinase X, as it has been observed numerous times [1,[3][4][5][6]. When the protein substrate R is phosphorylated to R p , the stability of R p is affected, which translates into an increasing or decreasing rate of degradation, k 4 . In particular, if the PTM stabilizes R p , i.e., k 4 decreases, the amount of R p increases, and R tot is expected to increase accordingly. This rise in R tot secondarily alters the degree of saturation of the phosphorylation and dephosphorylation reactions and, consequently, is expected to affect the degree of ultrasensitivity in the DR curves. These overall effects could theoretically also be caused by changes in k 3 , but we focus on k 4 because the PTM directly affects the stability of R p , whereas R is affected only in a secondary manner.

Effects on Steady-State R
When the PTM increases the stability of R p , i.e., k 4 decreases, the steady−state DR curves of R vs. X ( Figure 3A) and R p vs. X ( Figure 3B) both become steeper; conversely, when the stability of R p decreases, i.e., k 4 increases, the two curves become shallower. The changes in the steepness of the DR curves can be quantified by n H as well as the maximal local ultrasensitivity, |LRC| max . Both increase as k 4 decreases ( Figure 3D,E). Interestingly, however, for k 4 values comparable to or below the default value, |LRC| max is generally lower than |n H | for the R vs. X response, which is an indication that the Hill coefficient overestimates the maximal degree of signal amplification in these situations ( Figure 3D). For k 4 values rising above the default value, |LRC| max starts to match up with |n H | and eventually exceeds it. For very large k 4 values, |n H | approaches a constant value of approximately 1.58 and |LRC| max approaches a constant value of approximately 1.72. Thus, there is still ultrasensitivity, but its degree is modest. The value of the kinase activity X at which |LRC| is maximal shifts to the left as k 4 increases. changes in the steepness of the DR curves can be quantified by nH as well as the maximal local ultrasensitivity, |LRC|max. Both increase as k4 decreases ( Figure 3D,E). Interestingly, however, for k4 values comparable to or below the default value, |LRC|max is generally lower than |nH| for the R vs. X response, which is an indication that the Hill coefficient overestimates the maximal degree of signal amplification in these situations ( Figure 3D). For k4 values rising above the default value, |LRC|max starts to match up with |nH| and eventually exceeds it. For very large k4 values, |nH| approaches a constant value of approximately 1.58 and |LRC|max approaches a constant value of approximately 1.72. Thus, there is still ultrasensitivity, but its degree is modest. The value of the kinase activity X at which |LRC| is maximal shifts to the left as k4 increases.

Effects on Steady-State Rp
The elevated steepness of the Rp vs. X response, with increased stability of Rp, is evidently due to the increasing maximal Rp level when k4 decreases ( Figure 3B). Interestingly, and contrary to the effect on the response of R vs. X, LRCmax is generally higher than nH for k4 values below the default value, indicating that the Hill coefficient is underestimating the maximal degree of signal amplification ( Figure 3E). For k4 values higher than the

Effects on Steady-State R p
The elevated steepness of the R p vs. X response, with increased stability of R p , is evidently due to the increasing maximal R p level when k 4 decreases ( Figure 3B). Interestingly, and contrary to the effect on the response of R vs. X, LRC max is generally higher than n H for k 4 values below the default value, indicating that the Hill coefficient is underestimating the maximal degree of signal amplification ( Figure 3E). For k 4 values higher than the default value, LRC max starts to match n H and eventually drops below its value. For very large k 4 values, n H approaches 1.58, whereas LRC max settles at approximately 1. The value of kinase activity X for which LRC is maximal shifts to the left as k 4 increases.

Effects on Steady-State R tot
In the Goldbeter-Koshland model of the CMC, either R or R p is regarded as the output, because the activities of either one may change by the phosphorylation status. However, in some situations, the covalent modification status of an amino acid residue may only affect protein stability without affecting protein activity [35,36]. In these cases, R tot should be viewed as the output. Depending on the values of k 4 , the response of R tot vs. X can be either stimulatory or inhibitory ( Figure 3C), because either more or less R p is removed from the system. At the default level of k 4 , which is equal to k 3 , R tot does not change with X. However, as the PTM stabilizes R p , i.e., k 4 decreases from the default value, the steady-state response of R tot vs. X increases monotonically to a higher plateau than before and also becomes increasingly steeper, with LRC max surpassing n H for very low k 4 values ( Figure 3F). Conversely, as the PTM destabilizes R p , the steady-state response of R tot vs. X decreases monotonically toward a lower plateau and also becomes increasingly more sigmoidal (despite that the response of R p itself is no longer ultrasensitive), with |LRC| max approaching 1.72 for very high k 4 values. Surprisingly, |n H | changes in the opposite direction to |LRC| max for k 4 values above the default value ( Figure 3F). A small increase in k 4 above the default value first results in a very high |n H |, but as k 4 increases further, |n H | drops back and approaches 1.58. This inverse relationship between |LRC| max and |n H | demonstrates again that these two metrics do not always conform to each other and that reliance on the Hill coefficient as an estimate of the degree of signal amplification can be misleading. In summary, both stabilization and destabilization of R p can lead to the enhancement of ultrasensitivity in the steady-state response curve of R tot vs. X.
While R p and R are expected to exhibit ultrasensitivity due to the zero-order covalent modification effect, as revealed by the Goldbeter-Koshland model, it is interesting to note that R tot also exhibits various degrees of ultrasensitivity depending on the value of k 4 , i.e., the stability of R p . To dissect this mechanism leading to ultrasensitivity for R tot , we use the following two steady-state flux and mass conservation equations to solve for R tot : By substituting either R or R p from Equation (5) in Equation (6), we obtain two equations that exhibit symmetry with respect to k 3 and k 4 , namely The equations say that except for cases where k 3 and k 4 are equal, the steady-state R tot scales linearly with both R p or R. When k 3 > k 4 , i.e., phosphorylation results in R p stabilization, R tot has a basal level determined by k 0 /k 3 and increases as R p increases (Equation (7)). For very small k 4 , R tot ≈ k 0 /k 3 + R p . Since the response curve R p vs. X is always monotonically increasing ( Figure 3B), its ultrasensitivity is passed to R tot with comparable n H values. By contrast, the LRC of the R tot response will be lower than that of the R p response due to the presence of the basal level at k 0 /k 3 ( Figure 3E vs. Figure 3F). Conversely, if phosphorylation results in R p destabilization, i.e., k 3 < k 4 , R tot has a minimal level determined by k 0 /k 4 (Equation (8)). For very large k 4 , R tot ≈ k 0 /k 4 + R.
Since the response curve of R vs. X is always monotonically decreasing ( Figure 3A), its ultrasensitivity is passed to R tot with comparable n H values, and again, the |LRC| of the R tot response is lower than that of the R response, due to the presence of the minimal level k 0 /k 4 ( Figure 3D vs. Figure 3F).

Effects on Timing of Signaling
PTMs can have an effect on the timing of signaling. When they induce changes in protein stability, the time it takes the signaling motif to reach steady state in response to X is no longer determined only by the covalent modification reactions, but also by the half-lives of the protein substrate. Not surprisingly, for k 4 lower than the default value, it takes much longer time for R, R p and R tot to reach their steady state ( Figure 3G-I). The trajectory of R is nonmonotonic-it first decreases quickly as a result of the phosphorylation of pre-existing R and then rises slowly (because R tot increases) to settle at a new steady state ( Figure 3G). In comparison, R p first shoots up quickly as a result of the phosphorylation of pre-existing R into R p , and then rises slowly toward its new steady state ( Figure 3H). R tot does not exhibit a biphasic trend and instead increases gradually toward its new steady state ( Figure 3I). For k 4 higher than the default value, the time it takes to reach the steady state does not appear to be monotonically correlated with k 4 (Figure 3G-I). For k 4 values slightly higher than k 3 , the differential stability of R and R p causes the system to approach the steady state slowly because the protein half-life, rather than the fast MM reactions, dominates the long-term kinetics ( Figure 3G,H, purple vs. orange lines). However, as k 4 increases further, the responses are overall faster since the overall protein half-life becomes shorter ( Figure 3G-I, green vs. purple lines). Generally, R first decreases quickly as a result of phosphorylation of pre-existing R and then continues to decrease till it settles to a new steady state ( Figure 3G). In comparison, R p exhibits a nonmonotonic trajectory-it first rises quickly as a result of phosphorylation of pre-existing R into R p , and then decreases (because R tot decreases) slowly to settle at a new steady state ( Figure 3H). R tot has a similar monotonically decreasing profile as R ( Figure 3I).

Behaviors of the Full Model
For the full model, the steady-state DR behavior of R, R p and free R tot (R + R p ) with respect to X tot (the input parameter of the full model) and their dynamical responses are nearly identical to the simple MM model (simulation results not shown), indicating that the enhancement of ultrasensitivity by protein stabilization is a robust feature that is insensitive to implementation details under the current parameter conditions where the protein substrate is in excess of the enzymes (i.e., for the default setting R tot = 100 vs. Y tot = 1). While such strong differences in molecular abundances may be encountered in cell signaling and metabolic pathways, this situation is not general. Instead, the covalent modification enzymes may be at comparable levels to their protein substrates including those involved in the MAPK signaling cascade [27,28], suggesting for simulation use of the full model rather than the simplifying MM model. In particular, the apparent Michaelis constants K m1 and K m2 are substituted with their original definitions as (k 1b + k 1c )/k 1f and (k 2b + k 2c )/k 2f , respectively. The consequence of increasing the enzyme level is shown in Figure 4: increasing Y tot to 100, up to a point where Y tot = R tot , markedly weakens the degree of ultrasensitivity. Absent PTM-induced changes in R p stability, both R and R p are now only marginally ultrasensitive in response to X tot ( Figure 4A,B,D,E). However, when stabilization of R p occurs with lower k 4 values, the ultrasensitivity of R p is enhanced to some extent while the ultrasensitivity of R remains basically unchanged. Compared with the MM model (Figure 3 B,E), the enhancement of the ultrasensitivity of R p in the full model is not as strong. The reason is that the increase in the total protein substrate abundance, which can be achieved due to stabilization of R p in the full model, is not as high as can be achieved in the MM model. In the full model, a significant amount of R is titrated in the complex with enzyme X and thus not stabilized. As a result, the total protein substrate, including free substrates and substrates complexed with the two enzymes, only approaches 2.3-fold of the basal level when k 4 is decreased by 10-fold from the default value (results not shown), in contrast to the 10-fold increase in the MM model ( Figure 3C). Interestingly, when free R tot (R + R p ) is considered as the model output, it exhibits a monotonically decreasing but slightly ultrasensitive DR relationship at the default k 4 value ( Figure 4C). This outcome occurs because of sequestration of R and R p by the enzymes X and Y, respectively, as X tot increases. More interestingly, with stabilization of R p where k 4 value is further lowered, a nonmonotonic DR curve for free R tot emerges, which results from a rising total abundance of the protein substrate due to stabilization at higher X tot levels. When destabilization of R p occurs, the DR profile of free R tot follows R, the dominant free form, which is monotonically decreasing. outcome occurs because of sequestration of R and Rp by the enzymes X and Y, respectively, as Xtot increases. More interestingly, with stabilization of Rp where k4 value is further lowered, a nonmonotonic DR curve for free Rtot emerges, which results from a rising total abundance of the protein substrate due to stabilization at higher Xtot levels. When destabilization of Rp occurs, the DR profile of free Rtot follows R, the dominant free form, which is monotonically decreasing.

Protein Stabilization Can Lead to Emergence of Ultrasensitivity
As we demonstrated for a CMC with pre-existing ultrasensitivity, stabilization of Rp can enhance the degree of ultrasensitivity of the responses. In this section, we explore the possibility that stabilization of Rp can render a formerly non-ultrasensitive CMC ultrasensitive. To demonstrate this possibility, we first eliminate ultrasensitivity by raising the default values of the Michaelis constants 10-fold, such that in the MM model Km1 = Km2 = 100, a value that is the same as the substrate Rtot level with the default value of k4 at 0.01. As a result, the cycle no longer exhibits ultrasensitivity ( Figure 5A-C), as evaluated by |LRC|max ( Figure 5D-F). Starting with this new baseline, we now let k4 decrease below 0.01, which causes Rp to be more stable than R. Indeed, the responses, especially the steady-state DR curves for Rp vs. X and Rtot vs. X, all begin to show a trend toward ultrasensitivity, as the total protein substrate level approaches and eventually surpasses the Michaelis constants Km1 and Km2, thereby pushing the phosphorylation and dephosphorylation cycle toward saturation ( Figure 5A-C). These results demonstrate that ultrasensitivity can emerge de novo with PTM-induced protein stabilization.

Protein Stabilization Can Lead to Emergence of Ultrasensitivity
As we demonstrated for a CMC with pre-existing ultrasensitivity, stabilization of R p can enhance the degree of ultrasensitivity of the responses. In this section, we explore the possibility that stabilization of R p can render a formerly non-ultrasensitive CMC ultrasensitive. To demonstrate this possibility, we first eliminate ultrasensitivity by raising the default values of the Michaelis constants 10-fold, such that in the MM model K m1 = K m2 = 100, a value that is the same as the substrate R tot level with the default value of k 4 at 0.01. As a result, the cycle no longer exhibits ultrasensitivity ( Figure 5A-C), as evaluated by |LRC| max (Figure 5D-F). Starting with this new baseline, we now let k 4 decrease below 0.01, which causes R p to be more stable than R. Indeed, the responses, especially the steady-state DR curves for R p vs. X and R tot vs. X, all begin to show a trend toward ultrasensitivity, as the total protein substrate level approaches and eventually surpasses the Michaelis constants K m1 and K m2 , thereby pushing the phosphorylation and dephosphorylation cycle toward saturation ( Figure 5A-C). These results demonstrate that ultrasensitivity can emerge de novo with PTM-induced protein stabilization.
For comparison, we used the full model to validate the results of the simple MM model. Again, we eliminated ultrasensitivity by setting apparent K m1 = K m2 = 100, which was accomplished by increasing the values of both k 1b and k 2b to 990. With these settings, simulations show results that are nearly identical to those of the simple MM model; in particular, ultrasensitivity emerges for both free R p and free R tot ( Figure S1). We also explored the situation where ultrasensitivity is eliminated from the full model by increasing the phosphatase Y tot level. When Y tot is gradually increased to 300, the ultrasensitivity of R p basically disappears ( Figure S2B). At this new baseline level, Y tot has a value 3 times the basal R tot level. If k 4 is now decreased to mimic stabilization of R p , ultrasensitivity re-emerges for R p , albeit only marginally, as indicated by |LRC| max (Figure S2B). The R response is also ultrasensitive but the degree is not affected by k 4 ( Figure S2A). Free R tot exhibits a similar DR profile and ultrasensitivity as R ( Figure S2C). For comparison, we used the full model to validate the results of the simple MM model. Again, we eliminated ultrasensitivity by setting apparent Km1 = Km2 = 100, which was accomplished by increasing the values of both k1b and k2b to 990. With these settings, simulations show results that are nearly identical to those of the simple MM model; in particular, ultrasensitivity emerges for both free Rp and free Rtot ( Figure S1). We also explored the situation where ultrasensitivity is eliminated from the full model by increasing the phosphatase Ytot level. When Ytot is gradually increased to 300, the ultrasensitivity of Rp basically disappears ( Figure S2B). At this new baseline level, Ytot has a value 3 times the basal Rtot level. If k4 is now decreased to mimic stabilization of Rp, ultrasensitivity reemerges for Rp, albeit only marginally, as indicated by |LRC|max ( Figure S2B). The R response is also ultrasensitive but the degree is not affected by k4 ( Figure S2A). Free Rtot exhibits a similar DR profile and ultrasensitivity as R ( Figure S2C).

Regulation of Protein Modification Cycles through Alterations in Enzyme Features
Given the important role of enzyme saturation by the substrate in CMC-mediated ultrasensitivity, we explore further in this section whether changes in the kinetic features of the modifying or demodifying enzymes can modulate the DR curves and their ultrasensitivity. Specifically, we investigate how changes in the Michaelis constants Km1 and Km2 modulate the steady-state DR curves and their ultrasensitivity. As a first example, we consider Km1 and examine the case where phosphorylation of R to Rp results in destabilization, using as the baseline k4 = 0.1, which is 10-fold greater than k3. As Km1 decreases, the DR curves for R and Rtot in the MM model become increasingly more sigmoidal ( Figure  6A,C), with limited changes in the Rp responses ( Figure 6B). For low Km1 values, |LRC|max can be much greater than |nH|, whereas for high Km1 values, |nH| approaches 1.12, and |LRC|max approaches 1, indicating loss of ultrasensitivity ( Figure 6A,D). For the Rp response, increasing Km1 reduces the steepness of the DR curve with |nH| approaching 1.25, and ultrasensitivity is lost for high Km1 values as indicated by |LRC| below 1 ( Figure 6B,E). Lastly, increasing Km1 reduces the steepness of the DR curve for Rtot with |nH| approaching 1.25, and ultrasensitivity is lost for high Km1 values as indicated by |LRC| below 1 ( Figure  6C,F). Varying the Michaelis constant Km2 of the phosphatase has a similar effect on

Regulation of Protein Modification Cycles through Alterations in Enzyme Features
Given the important role of enzyme saturation by the substrate in CMC-mediated ultrasensitivity, we explore further in this section whether changes in the kinetic features of the modifying or demodifying enzymes can modulate the DR curves and their ultrasensitivity. Specifically, we investigate how changes in the Michaelis constants K m1 and K m2 modulate the steady-state DR curves and their ultrasensitivity. As a first example, we consider K m1 and examine the case where phosphorylation of R to R p results in destabilization, using as the baseline k 4 = 0.1, which is 10-fold greater than k 3 . As K m1 decreases, the DR curves for R and R tot in the MM model become increasingly more sigmoidal ( Figure 6A,C), with limited changes in the R p responses ( Figure 6B). For low K m1 values, |LRC| max can be much greater than |n H |, whereas for high K m1 values, |n H | approaches 1.12, and |LRC| max approaches 1, indicating loss of ultrasensitivity ( Figure 6A,D). For the R p response, increasing K m1 reduces the steepness of the DR curve with |n H | approaching 1.25, and ultrasensitivity is lost for high K m1 values as indicated by |LRC| below 1 ( Figure 6B,E). Lastly, increasing K m1 reduces the steepness of the DR curve for R tot with |n H | approaching 1.25, and ultrasensitivity is lost for high K m1 values as indicated by |LRC| below 1 ( Figure 6C,F). Varying the Michaelis constant K m2 of the phosphatase has a similar effect on ultrasensitivity ( Figure S3). The full model produces very similar responses when the apparent K m1 is varied by either changing k 1f or k 1b , even for the situation where the basal substrate and enzyme levels are comparable (e.g., Y tot = 100; simulation results not shown).
The rationale for a second analysis is the situation where phosphorylation of R into R p results in strong protein stabilization (k 4 = 0.001, 10-fold lower than k 3 ). When K m1 decreases below its baseline value of 10 in this situation in the MM model, the DR curves for R, R p and R tot become increasingly sigmoidal. For the response of R, |n H | obviously overestimates the degree of ultrasensitivity as evaluated by |LRC| max (Figure 7A,D). By contrast, for high K m1 values, |n H | approaches 1.93, and |LRC| max approaches 1, indicating loss of true ultrasensitivity. For the R p response, increasing K m1 reduces the steepness of the DR curve with |n H | approaching 2.61, and |LRC| max is reduced to 4.97 with some, but not a complete loss of ultrasensitivity ( Figure 7B,E). Except for very high K m1 values, |LRC| max is generally higher than |n H |. The reason that large K m1 values do not result in complete loss of ultrasensitivity is that K m2 is still kept at default value of 10, thus keeping the dephosphorylation step close to saturable. Lastly, increasing K m1 reduces the steepness of the DR curve for R p with |n H | approaching 2.61, while |LRC| max is reduced to 2.34 with some loss of ultrasensitivity ( Figure 7C,F). Except for very low K m1 values, |LRC| max is generally higher than |n H |. Varying K m2 has a similar effect on ultrasensitivity ( Figure S4). The full model produces very similar responses when the apparent K m1 is varied by either changing k 1f or k 1b (simulation results not shown). For the condition when the basal substrate and enzyme levels are comparable by setting Y tot = 100, the degree of ultrasensitivity of the R and R p responses is generally lower, but still increases as the apparent K m1 decreases. However, similar to Figure 4C, nonmonotonic dose response emerges for free R tot (simulation results not shown). The response patterns of ultrasensitivity modulation by apparent K m2 are generally similar.
Biomolecules 2021, 11, x FOR PEER REVIEW 13 of 20 ultrasensitivity ( Figure S3). The full model produces very similar responses when the apparent Km1 is varied by either changing k1f or k1b, even for the situation where the basal substrate and enzyme levels are comparable (e.g., Ytot = 100; simulation results not shown). The rationale for a second analysis is the situation where phosphorylation of R into Rp results in strong protein stabilization (k4 = 0.001, 10-fold lower than k3). When Km1 decreases below its baseline value of 10 in this situation in the MM model, the DR curves for R, Rp and Rtot become increasingly sigmoidal. For the response of R, |nH| obviously overestimates the degree of ultrasensitivity as evaluated by |LRC|max ( Figure 7A,D). By contrast, for high Km1 values, |nH| approaches 1.93, and |LRC|max approaches 1, indicating loss of true ultrasensitivity. For the Rp response, increasing Km1 reduces the steepness of the DR curve with |nH| approaching 2.61, and |LRC|max is reduced to 4.97 with some, but not a complete loss of ultrasensitivity ( Figure 7B,E). Except for very high Km1 values, |LRC|max is generally higher than |nH|. The reason that large Km1 values do not result in complete loss of ultrasensitivity is that Km2 is still kept at default value of 10, thus keeping the dephosphorylation step close to saturable. Lastly, increasing Km1 reduces the steepness of the DR curve for Rp with |nH| approaching 2.61, while |LRC|max is reduced to 2.34 with some loss of ultrasensitivity ( Figure 7C,F). Except for very low Km1 values, |LRC|max is generally higher than |nH|. Varying Km2 has a similar effect on ultrasensitivity ( Figure S4). The full model produces very similar responses when the apparent Km1 is varied by either changing k1f or k1b (simulation results not shown). For the condition when the basal substrate and enzyme levels are comparable by setting Ytot = 100, the degree of ultrasensitivity of the R and Rp responses is generally lower, but still increases as the apparent Km1 de- Figure 6. Effects of K m1 on ultrasensitivity under phosphorylation-induced protein destabilization in the MM model (k 4 = 0.1). (A-C) Steady-state DR curves for R vs. X, R p vs. X, and R tot vs. X, respectively, for different values of K m1 , as indicated in A. The same color scheme for K m1 values holds for all panels. The degree of ultrasensitivity increases for decreasing values of K m1 . (D-F) LRC (solid lines) and n H (dashed horizontal lines) for R, R p , and R tot , respectively. * K m1 = 10 is the default value.
In addition, we studied the effects of changing the catalytic constant k 2 (for the MM model) or k 2c (for the full model) of the phosphatase reaction on ultrasensitivity. In a nutshell, changes in k 2 barely affect the degree of ultrasensitivity when k 3 < k 4 ( Figure S5), ultrasensitivity increases considerably as k 2 increases when k 3 > k 4 ( Figure S6). Varying k 1 (for the MM model) or k 1c (for the full model) merely shifts the DR curves horizontally without changing the degree of ultrasensitivity (simulation results not shown).

Ultrasensitivity in Response to Changes in Protein Synthesis Rate
Lastly, we examined whether changes in the synthesis of R can lead to ultrasensitivity when PTM induces changes in protein stability. Suppose the kinase X displays an intermediate activity level of 1 matching the Y level, and the rate of synthesis of R, k 0 , is varied. Interestingly, when R p is destabilized in the MM model, i.e., k 4 > k 3 , R and R tot at steady state exhibit ultrasensitive responses for a certain range of values of k 0 even though their responses never plateau ( Figure 8A,C). By contrast, if k 0 is gradually increased, R p initially increases linearly (in log space), then plateaus, not exhibiting ultrasensitivity for any value of k 0 ( Figure 8B). When k 3 = k 4 , R tot is proportional to k 0 , and R is slightly ultrasensitive. For stabilization of R p , and thus k 3 > k 4 , the response of R vs. k 0 is linear, while the response of R tot vs. k 0 exhibits slight subsensitivity, with LRC dipping below 1 for some range of k 0 ( Figure 8F). In addition, we studied the effects of changing the catalytic constant k2 (for the MM model) or k2c (for the full model) of the phosphatase reaction on ultrasensitivity. In a nutshell, changes in k2 barely affect the degree of ultrasensitivity when k3 < k4 ( Figure S5), ultrasensitivity increases considerably as k2 increases when k3 > k4 ( Figure S6). Varying k1 (for the MM model) or k1c (for the full model) merely shifts the DR curves horizontally without changing the degree of ultrasensitivity (simulation results not shown).

Ultrasensitivity in Response to Changes in Protein Synthesis Rate
Lastly, we examined whether changes in the synthesis of R can lead to ultrasensitivity when PTM induces changes in protein stability. Suppose the kinase X displays an intermediate activity level of 1 matching the Y level, and the rate of synthesis of R, k0, is varied. Interestingly, when Rp is destabilized in the MM model, i.e., k4 > k3, R and Rtot at steady state exhibit ultrasensitive responses for a certain range of values of k0 even though their responses never plateau ( Figure 8A,C). By contrast, if k0 is gradually increased, Rp initially increases linearly (in log space), then plateaus, not exhibiting ultrasensitivity for any value of k0 ( Figure 8B). When k3 = k4, Rtot is proportional to k0, and R is slightly ultrasensitive. For stabilization of Rp, and thus k3 > k4, the response of R vs. k0 is linear, while the response of Rtot vs. k0 exhibits slight subsensitivity, with LRC dipping below 1 for some range of k0 ( Figure 8F). The emergence of ultrasensitivity in the responses of R and Rtot for high k4 values may be counterintuitive, since destabilization of Rp is believed to drive the enzymes away from saturation. The reason for ultrasensitivity to occur is the saturation of the flux through the phosphorylation (k1) step: when k0 approaches a high value like 10, any further small increase only leads to an increase in R, but not Rp, and the result is ultrasensitivity. Actually, this mechanism of ultrasensitivity is a variant of zero-order degradation, which no longer The emergence of ultrasensitivity in the responses of R and R tot for high k 4 values may be counterintuitive, since destabilization of R p is believed to drive the enzymes away from saturation. The reason for ultrasensitivity to occur is the saturation of the flux through the phosphorylation (k 1 ) step: when k 0 approaches a high value like 10, any further small increase only leads to an increase in R, but not R p , and the result is ultrasensitivity. Actually, this mechanism of ultrasensitivity is a variant of zero-order degradation, which no longer requires the dephosphorylation reaction. By setting k 2 = 0, i.e., disabling dephosphorylation, ultrasensitivity in the R and R tot responses remains strong ( Figure S7).
With the default parameter values, the full model produces responses that are nearly identical to the results from the MM model (simulation results not shown). In comparison, when both X tot and Y tot are increased to 100, the ultrasensitivity of both R and free R tot (R + R p ) in response to k 0 is greatly enhanced, especially with high destabilization of R p when k 4 = 1, compared with the MM model ( Figure S8A,C). This enhancement is likely due to the molecular titration effect by X tot and Y tot . However, when there is no change in stability of R p , or when it is stabilized, mild ultrasensitivity also occurs for R and free R tot , which is missing in the MM model. Moreover, R p also exhibits mild ultrasensitivity in response to k 0 , a feature that is missing in the MM model ( Figure S8B). A summary of the main results is provided in Table 3.
Note: ↑, ↓, and -denote increase, decrease, and no/very small effect, respectively. Unless otherwise indicated, the effects apply to both MM and full models. # Change that occurs with MM model only. * Change that occurs with full model only.

Discussion
Cellular signal transduction pathways and gene regulatory networks regularly involve PTMs of protein components as a means of regulating their activities and abundances. Nearly all PTM reactions require participation of specific enzymes that add or remove particular functional groups to the appropriate protein substrates. When these enzymes operate near saturation with respect to their substrates, nonlinear signaling may occur, where input signals are amplified and switch output signals on or off [16,17]. When the protein substrates in a CMC are in excess relative to the modification or demodification enzymes, the degree of saturation of these enzymes depends on the Michaelis constants and the abundances of the contributing substrates.
The covalent modification status of a protein substrate not only modulates its activity, but may also alter its affinity as a substrate for the ubiquitination-proteasomal pathway that mediates the degradation of the majority of intracellular proteins [37]. Depending on whether the covalently modified protein molecule is a better or less suited substrate for ubiquitination, PTMs can either stabilize or destabilize a protein and thereby regulate its abundance. For instance, under normoxia, HIF-1α is oxidized by prolyl hydroxylase domain-containing proteins (PHD) in an oxygen-dependent manner and thereby targeted by the pVHL ubiquitination pathway for degradation, thus keeping the hypoxic transcriptional program under control [36,38]. As a different example, phosphorylation of p53 by ATM during the DNA damage response leads to its stabilization [1]. Therefore, the overall protein half-life and abundance do not remain constant in these situations, rather, they can change dynamically depending on the covalently modified fraction of the protein molecules. The altered protein substrate abundance in turn affects the degree of enzyme saturation, and hence creates an important nonlinearity in signaling.
A paradigm scenario of this type is PTM-induced protein stabilization on top of zero-order ultrasensitivity that pre-exists even for basal abundances of the protein substrates. In this scenario, as our simulations demonstrate, the degree of ultrasensitivity for the phosphorylated protein response (R p ) with respect to the kinase X is considerably elevated, with LRC and the Hill coefficient increasing sharply as the half-life of R p is prolonged ( Figure 3B,E). The enhancement of ultrasensitivity is due to the concurrently increased total protein substrate abundance as the input signal X increases, which pushes the kinase and phosphatase further into a saturated mode of operation. When the protein substrate is not high enough to enable zero-order ultrasensitivity at the basal condition, the increased protein substrate abundance induced by PTM can move the signaling motif toward saturation, thereby causing the emergence of ultrasensitivity, as demonstrated in Figures 5 and S1. During the process of PTM-induced protein stabilization, the unmodified protein response is also enhanced for ultrasensitivity ( Figure 3A,D) or rendered ultrasensitive (Figures 5A,D, and S1A,D) although the response of R vs. X follows an inhibitory profile where R decreases as the input signal X increases.
An unexpected finding is the total protein response to the input signal (R tot vs. X for the MM model and free R tot vs. X tot for the full model), which can also exhibit ultrasensitivity, for both cases of PTM-induced protein stabilization and destabilization ( Figures 3C, S1C and S2C). The original Goldbeter-Koshland model was intended to examine either the covalently modified or unmodified protein responses under the condition of zero-order ultrasensitivity, while the total protein abundance stayed constant. Here, our simulations show that ultrasensitivity can emerge when there is an imbalance in the stability of the modified and unmodified proteins. When the modified protein is more stable, the total protein mass is dominated by the modified protein and thus resembles the modified protein response albeit with a non-zero basal level. By contrast, when the modified protein is less stable, the total protein response is dominated by the unmodified protein and thus resembles the unmodified protein response. In both situations, the response of R tot vs. X can be ultrasensitive because the modified or unmodified protein response is ultrasensitive. As an example, in the drosophila embryo, MAPK can phosphorylate transcriptional repressor Yan in response to morphogen gradients and thereby induce its degradation; this inducible degradation of Yan was proposed as part of a zero-order ultrasensitivity mechanism for the switch-like Yan response, which is responsible for the patterning of the embryonic ventral ectoderm [19]. Therefore, protein activity changes by PTM in a CMC are not mandatory for achieving zero-order ultrasensitivity if protein stability is also regulated by PTM. Using the full model, we also found that when there is a sufficient amount of demodifying enzyme to titrate the substrate, even under the condition that PTM results in protein stabilization, free total protein substrate R tot (R + R p ) can display a mildly ultrasensitive, decreasing DR with respect to X tot (Figures 4C and S2C). Moreover, if R p is further stabilized, free R tot exhibits nonmonotonic responses ( Figure 4C), which further demonstrates that PTM-induced changes in protein stability can bring about more complex dose-response patterns. We also demonstrate that if the input-driving signal is supposed to increase the production rate of the protein substrate, a saturable covalent modification reaction, coupled with decreased stability of the modified protein, can also lead to an ultrasensitive increase in either the unmodified or total protein levels ( Figures 8A,C and S8A,C); this result confirms a recent finding by Mallela et al. [25]. However, we additionally found that the ultrasensitivity is enhanced in the full model and also applies to the modified protein ( Figure S8). One uncertainty in our full model lies in the treatment of the stability of the protein substrates in complexes with their catalyzing enzymes. Since a protein in a multimeric protein complex can be selectively ubiquitinated and degraded, leaving the remaining components intact [30][31][32], it seems reasonable to assume that the protein substrate in complex with its enzyme is degraded with equal rate constant as the free substrate. It is of course possible that the stability of protein substrate in the complex can differ from that of its free form. The impact of this uncertainty is expected to be small if the substrate-enzyme complex is only a small fraction of the total protein substrate. However, if the fraction is significantly large, depending on the direction of the stability change, additional dynamics and DR response profiles may emerge, which will warrant further investigation.
In the absence of PTM-induced changes in protein stability, the CMC motif can launch a quick response amenable to the time scale associated with covalent modification reactions catalyzed by enzymes. However, when protein stability is altered by PTM with half-lives at the order of hours, it can take much longer for this signaling motif to reach a steady state ( Figure 3G-I). If the protein substrate or its downstream target is a transcription factor, such as p53, HIF-1, BCL-6 or Yan, a relatively slow rise or activation may not matter much as far as the timeliness of a response is concerned, because the ensuing transcriptional induction of downstream genes takes much more time to complete anyway. Importantly, we propose here that ultrasensitivity through protein stabilization can be a potential energy-saving strategy employed by cells, where maintaining a high, saturating level of the protein substrate at basal condition may no longer be necessary. In addition, the initial overshoot exhibited by the R or R p response, as shown in Figure 3G,I, can also be a signaling strategy utilized by cells to accelerate transcriptional induction for gene production with long half-lives [39].
Throughout the result section and the Supplemental Materials, we have compared the degree of steepness of the steady-state DR curve as quantified by n H with the degree of true ultrasensitivity quantified by LRC and confirmed their known differences in describing ultrasensitive DR curves [12,34]. While the two metrics in most situations move in the same direction in response to changes in a parameter value, the corresponding |n H | for a particular DR curve can be higher or lower than |LRC| max . A higher |n H | value means an overestimate of the degree of amplification of the DR curve, which often occurs when the DR curve has a significant basal level. There are also scenarios where the DR curve exhibits a profile comprising of an almost linear response followed immediately by a plateau (Figures 8B and S5B). Such a response profile may have an apparent n H = 2 despite the fact that its response is at most linear. We have also encountered DR curves having an |LRC| max value higher than |n H | ( Figure S6B,C); in these situations, n H underestimates the degree of amplification. Therefore, when examining the ultrasensitivity of a DR curve, it is advisable to quantify both LRC and n H .
Building upon Goldbeter and Koshland's concepts, Mallela et al. recently proposed mathematical models for protein modification cycles, focusing, in particular, on protein substrates that are ubiquitinated by the same E3 ligases, which mark both proteins for degradation [25]. Many E3 ligases are apparently promiscuous, thereby permitting competition between "similar" protein substrates. The authors observed that the sensitivity to incoming signals, as well as the ultrasensitivity of the response, is diminished or even destroyed when the protein substrate saturates the modifying enzyme. This ultrasensitivityweakening effect is more dramatic if the cycling proteins are degraded at a relatively high rate, consistent with our earlier findings [24]. However, even though the authors used the full model, their study only considered the situation when the total protein substrate is at least 1000-fold higher than the demodifying enzyme, while in the present study we systemically examined the condition when these two quantities are comparable or even when the enzyme is at a higher level. Mallela and colleagues also found that signaling cycles, in which the coupling of protein substrates collectively leads to saturation of the enzymes, can lead to a coupled, switch-like response in all protein substrates, likely due to the competition or "crosstalk" of the substrate proteins with respect to the same E3 ligases. The effects of protein turnover on ultrasensitivity does not seem to be limited to the CMC motif explored here. It also plays a modulatory role in the ultrasensitivity arising from molecular titration or protein sequestration [40].
The signaling motif of a CMC can exhibit complex dynamic behaviors and has been extensively studied with computational means. Wang et al. investigated and decomposed the tunability of the zero-order ultrasensitivity [41]. Xu and Gunawardena examined some more realistic intracellular situations where multiple enzyme intermediates exist due to co-substrate binding for both reversible and irreversible reactions and found that these complications modulate the zero-order switching behavior [42]. The operation of the CMC in the face of protein expression noise has been explored more recently [43,44]. It seems important to have correlated expression of the paired modification and demodification enzymes to prevent switch flipping, and bifunctional enzymes in a CMC may be an ideal solution in this regard [44]. Using linear reactions of the modification and demodification reactions, Soyer demonstrated that the CMC motif, like negative feedback or incoherent feedforward loops, can exhibit transient or persistent dynamic responses depending on the difference in protein stability [45]. In the present study, we have added a new aspect by demonstrating that PTM-associated changes in protein stability, enzyme features, or protein synthesis offer the cell another level of sophistication regarding the complex response behaviors of this long-studied signaling motif.