Skip to main content
Advertisement
  • Loading metrics

Crosstalk and ultrasensitivity in protein degradation pathways

  • Abhishek Mallela,

    Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

    Affiliation Department of Mathematics, University of California Davis, Davis, California, United States of America

  • Maulik K. Nariya,

    Roles Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – review & editing

    Affiliation Laboratory of Systems Pharmacology, Harvard Medical School, Boston, Massachusetts, United States of America

  • Eric J. Deeds

    Roles Conceptualization, Formal analysis, Investigation, Methodology, Supervision, Writing – original draft, Writing – review & editing

    deeds@ucla.edu

    Affiliations Department of Integrative Biology and Physiology, University of California, Los Angeles, Los Angeles, California, United States of America, Institute for Quantitative and Computational Biosciences, University of California, Los Angeles, Los Angeles, California, United States of America

Abstract

Protein turnover is vital to cellular homeostasis. Many proteins are degraded efficiently only after they have been post-translationally “tagged” with a polyubiquitin chain. Ubiquitylation is a form of Post-Translational Modification (PTM): addition of a ubiquitin to the chain is catalyzed by E3 ligases, and removal of ubiquitin is catalyzed by a De-UBiquitylating enzyme (DUB). Nearly four decades ago, Goldbeter and Koshland discovered that reversible PTM cycles function like on-off switches when the substrates are at saturating concentrations. Although this finding has had profound implications for the understanding of switch-like behavior in biochemical networks, the general behavior of PTM cycles subject to synthesis and degradation has not been studied. Using a mathematical modeling approach, we found that simply introducing protein turnover to a standard modification cycle has profound effects, including significantly reducing the switch-like nature of the response. Our findings suggest that many classic results on PTM cycles may not hold in vivo where protein turnover is ubiquitous. We also found that proteins sharing an E3 ligase can have closely related changes in their expression levels. These results imply that it may be difficult to interpret experimental results obtained from either overexpressing or knocking down protein levels, since changes in protein expression can be coupled via E3 ligase crosstalk. Understanding crosstalk and competition for E3 ligases will be key in ultimately developing a global picture of protein homeostasis.

Author summary

Previous work has shown that substrates of Post-Translational Modification (PTM) cycles can have coupled responses if those substrates share enzymes. This implies that modifications leading to substrate degradation (e.g. ubiquitylation by an E3 ligase) could introduce coupling in concentrations of substrates sharing a ligase. Using mathematical models, we found adding protein turnover to a PTM cycle diminishes both sensitivity and ultrasensitivity, particularly in models admitting long ubiquitin chains. We also found that proteins sharing an E3 ligase can indeed have coupled changes in both expression and sensitivity to signals. These results imply that accounting for crosstalk in protein degradation networks is crucial for the interpretation of results from a wide variety of common experimental perturbations to living systems.

Introduction

All proteins undergo some form of turnover. For instance, proteins can become damaged via deamidation or some other process and must be degraded in order to prevent unfolding and aggregation. Turnover is also important in signaling and the regulation of protein function. A classic example is the degradation of IκB proteins, which bind the NF-κB protein complex and sequester it in the cytoplasm. During response to different stimuli, IκBs are phosphorylated by IκB kinases, ubiquitylated, and tagged for degradation, which allows NF-κB to translocate to the nucleus [1]. In the cell, synthesis and degradation (i.e. protein turnover) act in concert to maintain an appropriate concentration of active protein (i.e. protein homeostasis). Given the centrality of protein turnover to all cellular processes, it is not surprising that dysregulation of protein homeostasis has been implicated in a vast array of neurodegenerative diseases and cancers [2,3]. In eukaryotes, degradation is often achieved through the ubiquitin-proteasome system, where proteins are tagged with polyubiquitin chains that are recognized by the proteasome, ultimately leading to protein degradation [4]. Polyubiquitylation represents a form of post-translational modification (PTM) cycle where ubiquitin subunits are covalently linked to substrates by E3 ligases and removed by deubiquitylating (DUB) enzymes [5].

Over 35 years ago, Goldbeter and Koshland studied the general properties of a PTM cycle comprised of a modifying and demodifying enzyme. They found that reversible cycles of protein modification, such as a kinase enzyme adding a phosphoryl group and a phosphatase enzyme removing it, work like on-off switches when the enzymes are saturated [6]. This phenomenon, known as “0th-order ultrasensitivity”, has had profound implications for understanding how biochemical networks can exhibit switch-like behavior. Despite decades of progress in understanding 0th-order ultrasensitivity and other aspects of PTM function [712], to date there have been few attempts to systematically characterize the general behavior of PTM cycles that drive protein degradation.

The only exception to this has been the study of ubiquitylation in the context of cell cycle oscillations and bistability [13,14]. While these studies have provided key insights about cell cycle control, they have not investigated how ubiquitylation controls the steady-state expression levels of proteins not involved in the cell cycle. It has also been shown that adding protein synthesis and degradation to models of gene expression and cell signaling can have dramatic effects on system dynamics, but the detailed impact of turnover on PTM cycles remains unclear [1517].

In addition to a general lack of understanding of the influence of protein homeostasis on PTM cycle behavior, we recently discovered that substrates in such cycles can have coupled steady-state responses if those substrates share modification/demodification enzymes. In particular, if one substrate is at saturating levels, or if the substrates collectively saturate the enzymes, then all substrates of that pair will respond in a coupled, switch-like manner [1820]. This implies that modification leading to substrate degradation (e.g. ubiquitylation by an E3 ligase) could introduce coupling in the concentrations of substrates sharing a ligase. Interestingly, Mather and colleagues have shown that substrate concentrations can be coupled through saturation of the downstream degradation machinery [21,22]. It is currently unclear, though, whether such coupling can arise due to “crosstalk” in the upstream mechanisms that tag proteins for degradation.

In this work, we used a set of mathematical models to show that perturbing a standard PTM cycle by simply adding synthesis and degradation has profound effects on the response of the system. Specifically, we found that the sensitivity of the system to incoming signals and the ultrasensitivity of the response are dramatically muted when the substrate is at saturating concentrations. When the modification in question drives protein degradation at a higher rate, these effects are even more pronounced. Furthermore, more realistic models allowing for long ubiquitin chains exhibit qualitatively similar behavior to the case with a single modification state, but with further decreases in sensitivity and ultrasensitivity. These findings are robust to changes in the specific mechanisms utilized by the E3 and DUB enzymes. Interestingly, we found that distinct modes of enzyme saturation (i.e. increasing substrate production rate vs. decreasing the Michaelis constant of the enzyme) also generate different substrate responses. This indicates that many classic results on PTM cycles, including the extremely ultrasensitive response they exhibit when the substrates are at saturating concentrations, may not hold in vivo where protein turnover is inevitable. We also found that proteins sharing an E3 ligase can indeed have closely related expression profiles. Moreover, the sensitivity of protein concentration to changes in E3 activity for any given protein is largely dependent upon the total expression level of the other proteins. This suggests that it may be difficult to interpret experimental results obtained from either overexpressing or reducing protein concentrations, since changes in protein expression can be coupled via E3 ligase crosstalk. Further experimental characterization of E3-ligase/DUB enzyme/substrate relationships will thus be vital to developing a global understanding of protein regulation within the cell.

Results

Competition among E3 ligases

As mentioned in the introduction, shared E3 ligases have the potential to induce coupling in substrate responses. It is currently unclear, however, how widespread such “crosstalk” among E3 ligases might be. We searched the E3Net database [23] for statistics of E3-substrate interactions in human cells. For sake of comparison, we also obtained E3-specific statistics from the hUbiquitome database [24]. As of this writing, the total number of E3 ligases documented in E3Net was 415 and the total number of their substrates was 873, making the average ‘substrate load’ (substrate-to-ligase ratio) 2.10. Similarly, there are a total of 138 ligases and 279 substrates annotated in the hUbiquitome database, yielding a comparable ratio of 2.02. Thus, on average, most E3 ligases will ubiquitylate around two substrates.

In addition to providing the numbers of ligases and substrates, E3Net also captures information on specific E3-substrate interactions. We found that 54% of the E3 ligases in the database have no substrates listed; however, of the remaining E3 ligases, 52% have at least 2 substrates and 11% have more than 10 substrates (Fig 1A). Also, the maximum number of substrates for any ligase is 92. Given that the database is incomplete, it is likely that these numbers represent significant underestimates of E3 ligase crosstalk. Regardless, the phenomenon of E3 ligases acting on multiple substrates is likely widespread, and little is presently known about what influence crosstalk might have on the responses of these substrates to changes in E3 ligase activity.

thumbnail
Fig 1. Crosstalk among E3 ligases & schematic diagrams of single-substrate models.

(A) Probability distribution (on log-log scale) of E3 ligase-substrate specificity as recorded in the E3Net database. The average “substrate load” on a given E3 ligase is 2.1. In Panels B-D, the general representation of the canonical Goldbeter-Koshland (GK) loop is shown, with protein turnover included. (B) The case of the traditional GK loop, with no synthesis and degradation; (C) The “Intermediate” model, which introduces synthesis (at a rate Q > 0) and degradation (at a rate δ1 > 0); (D) The “Full” model, which includes synthesis (Q > 0) and to different degradation rates: δ1 for the unmodified substrate and δ2 for the modified substrate. Since the modification in this model is meant to represent a tag that drives degradation, the rate of degradation for the modified substrate is higher than that of the unmofied substrate (i.e. δ2 > δ1 > 0). Here “M” denotes modifying enzyme and “D” denotes demodifying enzyme. Modified substrate is indicated by the dark, filled circle.

https://doi.org/10.1371/journal.pcbi.1008492.g001

Adding synthesis and degradation to a PTM cycle

Even though E3 ligases generally attach long ubiquitin chains to their substrates [4], in order to simplify the problem to an analytically tractable form, we first considered a case with just a single modification state. Because ubiquitylation actively effects protein degradation, any investigation of the interplay between substrates competing for a protein and post-translational modifications (PTMs) leading to degradation must account for protein turnover. We thus focused first on studying how synthesis and degradation influence the behavior of the standard Goldbeter-Koshland loop (Fig 1B).

The first model, which we call the ‘Intermediate’ model, involves one substrate that can exist in two forms: modified and unmodified, denoted by S* and S respectively (Fig 1C). In this model, the modification (e.g. phosphorylation) does not lead to higher rates of degradation, so the modified and unmodified substrates are both degraded at the same first order rate δ1. Unmodified substrate is also synthesized at a constant rate Q. There is a modification enzyme M that catalyzes the addition of the PTM in question, and a demodification enzyme D that catalyzes removal of the modification. Since these are enzymes, they form enzyme-substrate complexes: D forms a complex with S* and M with S. To simplify the model, neither M nor D are subject to synthesis and degradation. We assume that the enzyme-substrate complexes are subject to “degradation” at the same rate δ1, but degradation in this case simply removes the substrate, essentially freeing the enzyme (note that relaxation of this assumption has no impact on our results, see S1 Text, Sec. 1).

The enzymatic reaction scheme corresponding to this model can be used to obtain a system of ordinary differential equations (ODEs) with the binding, dissociation, and catalysis steps treated explicitly (full details are in S1 Text, Sec. 1.1). We have denoted the kinetic rates of complex formation, complex dissociation, and catalysis by kx,y, where x represents the reaction step and y represents the enzyme—modifying (M) or demodifying (D) (S1 Text, Sec. 1.1). For example, kcat,D denotes the catalytic rate of the reaction catalyzed by the demodifying enzyme.

Traditional analyses of post-translational modification cycles (i.e. the GK loop) have examined the response of molar fraction of modified protein at steady-state (S* ≡ [S*]/[S]T, where [S]T = [S] + [S*] + [MS] + [DS*]) to changes in the input parameter , which quantifies the activity of the M enzyme relative to that of the D enzyme [6,18]. We varied r by simply changing [M] and numerically integrated the system to extract the steady-state solutions for unmodified and modified substrate ([S] and [S*]) at each value of r. Note that, for the Intermediate model, [S]T = Q/δ1, regardless of the values of other parameters (S1 Text, Sec. 1.3). Here, we focus on the case where the Michaelis constants of the enzymes are equal (KM,M = KM,D). We should note that our models only consider the standard Michaelis-Menten case where enzyme concentrations are significantly lower than substrate concentrations. Relaxing this assumption could have an influence on the behavior of the models when the concentration of the enzymes is similar to that of the substrates. That being said, analysis of such a scenario is beyond the scope of this work, and we leave exploration of that case to the future [25].

One key feature of GK loops is their capacity to exhibit 0th-order ultrasensitivity, which manifests as a switch-like transition in [S*] vs. r when the modification and demodification enzymes are saturated (i.e. [S]TKM) [6,9,18,25]. To explore this phenomenon in the Intermediate model, we initially increased Q to change saturation levels, since [S]T = Q/δ1. In order to conduct these simulations, we chose a set of reasonable values for the kinetic rate constants in the model (Table 1). In particular, kcat’s and KM’s were taken from experimentally observed ranges (S1 Text, Sec. 1.4, [26]) and δ1 was set based on the average observed half-life for proteins in living human cells [27]. The value for δ1 is also very similar to the shortest observed protein half-life in mouse C2C12 cells [20,28]. As can be seen from Fig 2A, there are dramatic differences between a GK loop and the Intermediate model upon saturation. For instance, defining r50 as the r-value when [S*] is half-maximal (i.e. S* = 0.5), we see that the curve of S* vs. r shifts to the right, indicating a higher r50. Secondly, the ultrasensitivity (i.e. the effective Hill coefficient neff) of the system under the Intermediate model is reduced. Since the Intermediate model is simply a GK loop with synthesis and degradation (Q and δ1, compare Fig 1B and 1C) added, these results indicate that adding synthesis and degradation to a PTM can have a dramatic effect on 0th-order ultrasensitivity.

thumbnail
Table 1. Parameter values used for Fig 2A and 2B.

Fig 2A: chosen on the basis of , which is the average reported protein half-life in human cells; Fig 2B: KM = 101 × [S]T (unsaturated) and KM = 10−1 × [S]T (saturated).

https://doi.org/10.1371/journal.pcbi.1008492.t001

While increasing the expression level of S (e.g. increasing Q) is a natural way to achieve saturation, one can also saturate the enzymes by decreasing KM, keeping Q fixed. In a standard GK loop, varying [S]T and varying KM are mathematically equivalent; in the Intermediate model, however, the effects of decreasing KM (with a lower bound of 100 nM for experimentally observed KM’s, S1 Text, Sec. 1.4) are dramatically different from the effects of increasing Q. In particular, the changes in r50 and neff are negligible (Fig 2B).

thumbnail
Fig 2. Effects on various single-substrate models of varying Q or KM as the measure of enzyme saturation.

(A) Modulating the rate of protein synthesis (Q) results in dramatic reduction of both sensitivity of the system to incoming signals and ultrasensitivity of the response, in the regime where the enzymes are saturated. This is indicated by the rightward shift of the r50 and the reduction in the Hill coefficient (neff) from GK to Full. Note the logarithmic axis used for r. (B) Varying the Michaelis-Menten constant (KM) results in a smaller reduction of r50 and neff, as compared to varying Q, in the regime of saturated enzymes. In Panels A-B, dashed lines and solid lines correspond to the unsaturated and saturated cases, respectively. Also, for all of the cases in these panels, [D]T = 10−1 nM and [M]T is varied in order to vary r. In Panel A, [S]T = 103 nM for the unsaturated cases and [S]T = 105 nM for the saturated cases. In Panel B, [S]T = 103 nM for all cases. (C) Increasing [S]T does not change the r50 for the GK model, but the reduction in sensitivity is highly pronounced for the Intermediate model, and even more so for the Full model. Axes are on a log scale. (D) Increasing KM results in a dramatic reduction of neff for the GK model. For systems that incorporate protein turnover (i.e. the Intermediate and Full models), the effect of reduction occurs for sufficiently large KM.

https://doi.org/10.1371/journal.pcbi.1008492.g002

The findings above were obtained for just a single set of parameters, and this raises the question of whether or not our observations are robust to parameter variation. Since this model is relatively simple, we were able to characterize this parameter dependence analytically by solving the system of ODEs at steady state. We obtained the following equation relating r50 to Q and KM when corresponding kinetic rate parameters for the modifying and demodifying enzymes are identical (S1 Text, Sec. 1.6):

Considering the endpoints of the plots in Fig 2C and the equation above, it is clear that as the rate of substrate production is made arbitrarily large, the r50 grows without bound. Thus the system described by the Intermediate model always becomes less and less sensitive to incoming signals as Q increases, regardless of the value of the other parameters. However, when we make KM as small as possible with all other parameters fixed (i.e. KM → 0), we see that , which is a constant. Note that this constant value is nevertheless always larger than r50 = 1 for the GK model, independent of the choices of the other parameters (S1 Text, Sec. 1.6).

In a similar fashion we can analyze the effective Hill coefficient (neff) for the Intermediate model. Note that the S* vs. r curves do not precisely follow the form of a Hill function; as such, we use the standard definition [6,29]. Our analytical results establish a lower bound on neff for the Intermediate model (which we will refer to as neff(I)), indicating that the Intermediate model always exhibits positive cooperativity (i.e. neff(I) > 1, S1 Text, Sec. 1.7). As with r50, we also find that varying saturation levels by changing total substrate concentration (i.e. changing Q) vs. changing KM results in opposing effects on neff (Fig 2D). While neff(GK) grows without bound as saturation increases, neff(I) increases only modestly and is generally smaller by several orders of magnitude. For instance, when Q is increased, neff(I) evaluates to exactly 2, regardless of the values of the other parameters.

The above results clearly demonstrate that changing saturation by varying Q and varying KM have very different consequences for the steady-state response of the PTM cycle in the Intermediate model. Specifically, increasing saturation of the enzymes by increasing total substrate levels (in other words, by increasing Q) results in responses that are much less sensitive (i.e. r50 increases, requiring greater M activity in order to achieve an appreciable response) and much less ultrasensitive (i.e. lower neff) than would be predicted by a traditional GK loop (Fig 1B) [6,7]. The intuitive reason behind this has to do with the ultimate source of extreme ultrasensitivity in the original GK model. In that case, there is neither synthesis nor degradation, so steady state can only be achieved when the flux of S* production by M precisely matches the flux of S production by D. If substrate is at high concentration ([S]TKM), then both enzymes can theoretically operate at their Vmax. If Vmax,MVmax,D (i.e. r ≠ 1), this means that one of the enzymes can operate faster than the other. Say that Vmax,M > Vmax,D. In that scenario, the only way for the enzymatic fluxes to be balanced is for the M enzyme to operate below its Vmax. In other words, the M enzyme must convert so much of the substrate from the S to the S* state that it becomes unsaturated, meaning [S] ≪ KM. Since we are in a condition with [S]TKM, this directly implies that [S]/KM ≪ 1, and as a result S* will be close to 1. When the D enzyme has a higher Vmax (r < 1), then the situation is reversed, and S* will be close to 0 at steady state. This leads to a very switch-like behavior in S* as r is increased in the original GK loop (Fig 2).

In the Intermediate model, however, additional fluxes have been added, namely synthesis and degradation. In particular, steady state is not achieved when the flux of S* production by M matches the flux of de-modification by D, but rather when the flux of S* production matches the flux of de-modification plus the flux of degradation of S*. In this model, degradation is first-order, so it can never be saturated; as [S*] increases, so will the flux of degradation. As more and more substrate is added to the system, the velocity of substrate modification must thus increase in order to match this increasing flux of degradation, if [S*] is to be relatively high at steady-state. Under such conditions, the M enzyme will be operating at Vmax, and so the only way to increase the rate of modification of the substrate is to increase Vmax,M, thus increasing r. This leads to the increase in r50 as the steady-state levels of [S]T increase. Similarly, since degradation helps balance S modification at steady-state, there is no longer an extreme switch-like response in S* vs. r (Fig 2). Because varying KM changes saturation without altering the steady-state level of substrate, the effect of changing KM is very different from varying Q in the Intermediate model.

We should note that since the value of KM depends on the underlying rate constants for the enzyme-substrate interaction, it is unlikely to vary on short, cellular time scales. In other words, it is hard to imagine how a cell would increase the saturation of an enzyme by dynamically lowering the KM of one of its enzymes during the course of a response to an environmental fluctuation. On evolutionary time scales, however, an enzyme’s KM can change, subject to reasonable physical and biological constraints. So, it is more likely that saturation will change by changing the production rate Q in vivo. Certainly, experimental manipulation of saturation generally occurs through changes in protein expression (e.g. by “overexpressing” the protein, which would correspond to increasing Q in this model). Thus, even for cases like phosphorylation where the PTM does not directly influence degradation rate, the steady-state responses of PTM cycles in vivo may thus be quite different from the standard predictions that have been made in the absence of any consideration of protein turnover [712].

Driving protein degradation: The “Full” model

While the results described above hold for any PTM cycle subject to turnover, we are ultimately interested in PTMs like ubiquitylation that drive protein degradation. This corresponds in our case to δ2 > δ1, which we term the “Full” model. For the purposes of display, we kept δ1 close to the average degradation rate of human proteins and set δ2 close to the fastest degradation rate observed in human cells (i.e. δ1 = 2 × 10−5s−1 and δ2 = 2 × 10−4s−1) [27]. The point here is to focus our analysis on a case where the PTM drives rapid degradation, as is often thought to be the case when a substrate is ubiquitylated in vivo [4].

We first considered how changes in E3 ligase activity relative to DUB activity would influence the modification state of the substrate. To facilitate comparison with the Intermediate and GK models, we initially focused our analysis on the steady-state level of substrate modification, S*. We found that transitions in S* are even less sensitive to incoming signals in the full model, as compared to the Intermediate model (Fig 2A and 2B). Indeed, the r50 for the full model is always greater than that for the Intermediate model as Q is increased (Fig 2C), and we have shown analytically that this is true for any reasonable set of kinetic parameters (S1 Text, Sec. 1.6). Intuitively, in the Full model, the modified substrate has a higher degradation rate than the unmodified state, so that the system needs to be driven harder towards the modified state (i.e. by increasing r) in order to appreciably increase its concentration at steady state. Interestingly, although neff(I) is always less than neff(GK) as discussed above, we find that neff(Full) is less than neff(I) only for very small Q. For instance, when Q is increased without bound neff(Full) evaluates to exactly , or approximately 7 in our case, which is larger than limQ→∞ neff (I) = 2.

Since E3 ligase activity drives higher levels of protein degradation in the full model, changes in the r parameter will change not only [S*] but also the total concentration of substrate ([S]T). Perhaps not surprisingly, we found that [S]T also exhibits an ultrasensitive transition in r. As with the transitions in [S*] discussed above, there is a rightward shift in r50 for the [S]T vs. r curve as Q is increased in the full model (S1 Text, Sec. 1.8). This phenomenon can generate interesting behaviors, as shown in Fig 3A. Suppose that we systematically increase the expression level of the protein while keeping the concentrations of the modifying/demodifying enzymes constant, which corresponds to a constant r in this model. As Q increases, the r50 of the curve increases from a point less than this constant value of r to a point greater than r. This leads to nonlinear changes in total substrate as Q increases (see below). The [S]T vs. r curve has a number of other similarities to the S* vs. r curve; for instance, we see that the effective Hill coefficient for total substrate in the full model does not change significantly when Q is increased, and never exceeds a value of 2 (Fig 3B).

thumbnail
Fig 3. Effects of varying Q on the r50 and neff in the Full model and its representative analog for multiple modification states.

(A) The shift in r50 for the transition in total substrate at increasing values of Q is clearly evident here. Each dashed curve indicates a different value for Q. The maximum and minimum values of [S]T for each curve is Q/δ1 and Q/δ2, respectively. (B) The effective Hill coefficient neff is relatively unaffected by increase in protein synthesis rate in the Full model. (C) Schematic diagram for one substrate with multiple modification states. Shown here is the model corresponding to the Processive E3, Distributive/Sequential DUB case. Each of the first three units is degraded at the rate δ1, which is smaller than the rate δ2 for the remaining units. The maximum length of the polyubiquitin chain is denoted by . (D) In the Processive E3, Distributive/Sequential DUB model, much more E3 ligase activity is necessary for a maximal response in the saturated regime. (E) Compared to Panel (A), the rightward shift is more pronounced in the presence of polyubiquitin chains.

https://doi.org/10.1371/journal.pcbi.1008492.g003

So, in the Full model, as in the Intermediate model, increasing substrate levels by increasing Q results in a decrease in sensitivity (increase in r50) and limited ultrasensitivity in the response, regardless of whether we consider the response parameter to be S* or [S]T (which is likely the more relevant response parameter in vivo). The similarities in this case all arise from the fact that, in these models, steady state can only be achieved when the flux of the M enzyme matches the sum total of the fluxes due to the D enzyme and the degradation process (note Eq. 8 in S1 Text, Sec. 1.3). This gives rise to a fundamentally different phenomenology than is observed in the standard GK loop (Figs 2 and 3).

Adding multiple modification states to the full model

While the full model is suggestive, it abstracts a number of details of the biological systems that control protein homeostasis. For instance, E3 ligases, rather than adding just a single ubiquitin to their substrates, instead tend to attach polyubiquitin chains of varying lengths. Most of the available literature suggests that this ubiquitin chain needs to be at least four units long in order for the proteasome to efficiently recognize and degrade the substrate [13,30,31]. While the specific details of both the requisite length of the chain and the E3/DUB enzyme mechanisms are still being elucidated, it is fairly clear that a single ubiquitin is not sufficient to drive effective degradation by the proteasome.

To capture these effects in our models, we surveyed available literature and found that multiple enzymatic mechanisms have been proposed for both E3 ligases and DUB enzymes [13,3038]. E3 ligases may be “processive”, in the sense that the ligase adds an ubiquitin unit to the polyubiquitin chain at each catalytic step and stays attached to the substrate while multiple ubiquitins are added sequentially. Alternatively, they may be “distributive”, meaning that the ligase disassociates from the substrate at the end of each catalytic reaction. In the one case that has been extensively studied experimentally, a form of E3 called a RING ligase works with the E2 Cdc34 to build polyubiquitin chains on substrates in a processive manner [37]. Of course, this does not mean that other E3 ligases might not display distributive kinetics. Regarding the DUB enzyme counterpart, 3 such enzymes have been found in 26S proteasomes: Rpn11, Usp14, and Uch37 [3234,38]. Rpn11 functions by truncating at the base of the chain (in a distributive manner), whereas Usp14 and Uch37 serve primarily to trim the ubiquitin chains sequentially (in a processive manner). Interestingly, more than one DUB might act on a given chain [32].

Although there are experimentally characterized examples for several of these possible mechanisms, little is actually known about how widespread each mechanism may be in nature. We thus employed an exhaustive approach, examining all combinations of the enzyme mechanisms and creating models of those scenarios. In our initial analysis, the parameter values we used for the distributive cases correspond to the values in the previous section (i.e. Single Substrate, Single Modification State). However, we used parameter values directly from literature for the processive cases [37].

While we have analyzed all 6 possible combinations of enzyme mechanisms (see S1 Text, Sec. 2), given the available experimental data [37], we focus our discussion in the main text on a reasonable and representative case (Processive E3 and Distributive/Sequential DUB) from this set of models. We have depicted this scheme in Fig 3C. As mentioned above, current consensus indicates that a polyubiquitin chain typically requires at least four ubiquitin units to be efficiently degraded by the proteasome [4]. We thus assumed that each of the first three modification states (0–3 ubiquitins) is degraded at a uniform rate, δ1, which is smaller than the corresponding (higher) rate δ2 for each of the remaining states (4 or more ubiquitins). In theory, the ubiquitin chain could reach an infinite length, though of course in practice the action of DUBs and degradation will limit the largest chain typically observed in the system. Denoting this maximum length of the ubiquitin chain by , we enumerated the chemical reaction networks for all of the possible mechanistic scenarios described above (S1 Text, Sec. 2). Due to the inherent complexity of the models, we could not obtain closed-form analytical solutions, and thus focused on numerical simulations.

Recall that in the full model, which has a single modification state, we found a significant reduction in both sensitivity and ultrasensitivity of the transition in S* when compared to the Intermediate model. To compare our more complex model with the full model, we defined S* for the case with multiple modification states as follows: , where k indexes the substrate modification state. To choose a reasonable value for , we systematically increased this parameter and found a threshold value such that changes in r50 and neff were negligible beyond that threshold. Using this approach, we chose a value of 50 for heuristically by visual inspection. In Fig 3D and 3E, we see that the inclusion of ubiquitin chains magnifies the aforementioned effects in both the S* vs. r and [S]T vs. r curves. Specifically, much more E3 ligase activity is necessary to achieve a maximal response in saturated regimes.

As described above, numerical integration of the ODEs for all of the multiple modification state models necessitated the definition of a maximum possible length for the ubiquitin chain (). To determine if this truncation has any influence on the results, we developed an “agent-based” stochastic simulation (see S1 Text, Sec. 2.8). Rather than representing the system as a set of concentrations for each species, in these simulations we consider a finite population of explicit substrate “agents.” Each of these substrate molecules has associated with it the length of its ubiquitin chain, which ranges from 0 (no ubiquitin) to the largest number we can represent as a floating point number, allowing the length to be essentially arbitrary. Because of their discrete representation of the substrates, these simulations are stochastic, and we parameterized these simulations so that the stochastic rates correspond exactly to the parameters we used for our deterministic simulations (see S1 Text, Sec. 2.8) [3942]. We found excellent agreement between the deterministic results and the agent-based simulations for all of the enzyme mechanisms we considered, suggesting that truncating the system at = 50 yields a reasonable approximation.

Interestingly, all of the models that we examined, arising from the various mechanisms proposed for the E3 ligase and DUB enzymes, generated similar qualitative behavior (S1 Text, Sec. 2.8–2.9). Thus, while the quantitative details of the response (e.g. the value of r50 and neff) vary somewhat between the models, our general findings are largely invariant with respect to the catalytic mechanisms utilized by the E3 ligase and DUB enzymes [4]. The results presented above, however, all focus on a single set of parameters; while these parameters are reasonable, it is not clear how reasonable variation in the parameters might affect our findings. We focused on how variation in the KM and δ parameters might influence our findings, sampling these parameters from log-uniform distributions across two orders of magnitude centered on the set of parameters considered above. Again, while changes in these parameters influence the quantitative details of the response, the qualitative behaviors of these models are all consistent with our original observations. In particular, we found that saturation increases the value of r50 in the transition of both S* and [S]T in all cases (S1 Text, Sec. 2.8–2.9).

As in the case of the Full model, the fundamental reason we observe similar behavior in these multiple modification state models, regardless of the specific enzyme mechanism, has to do with the fact that all of those details simply modify the relationship between M enzyme activity and the flux of degradation of the substrate. In other words, regardless of the kinetic scheme, there are two populations of molecules: those degraded more slowly (δ1) and those degraded more rapidly (δ2). In order to reach a steady state, the flux of the M enzyme that produces these rapidly degraded species must match the flux of degradation plus the flux of de-modification by the D enzyme. As substrate concentration increases, more M activity is required to do this, leading to the increase in r50 and the overall lower level of ultransensitivity that we observe. Thus, while modification of both the parameters and the enzymatic scheme change some of the quantitative features of the response, the fundamental behavior is similar in all cases.

One interesting thing to note here is that PTM systems with multiple modifications can readily lead to multistability, with distinct stable steady-states corresponding to different distributions of the concentration of the various modification states [4344]. We did not observe this kind of phenomenon in any of our models, most likely because we assumed that the rates of the enzymatic reactions do not depend on the modification state itself (see S1 Text, Sec. 2.1–2.6). Relaxation of this assumption may lead to multistability in ubiquitylation states, which could have important consequences for the function of the ubiquitin-proteasome system. To our knowledge, there has been no experimental observation of multistability in ubiquitylation, likely due in part to the fact that no one has to date designed an experiment aimed at testing this idea explicitly. We leave the exploration of the potential for multistability in these systems to future work.

Adding multiple substrates to the full model

As mentioned above (Fig 3A), there is an increase in r50 for the [S]T vs. r curve as Q is increased in the full model. As a consequence, increasing Q while keeping r fixed at a value of 100 results in the curve seen in Fig 4A. For low values of Q, the transition in r50 occurs before this fixed r-value, so [S]TQ/δ2; for large Q, the transition in r50 occurs after this fixed r-value, so [S]TQ/δ1. For intermediate Q, however, there is a distinct transition between these two regimes. In other words, as Q increases, the system will naturally transition between a point “after” the transition in [S]T to one “before” the transition in [S]T, due to the impact that saturation has on r50. The result in Fig 4A implies that if two substrates share an E3/DUB enzyme pair, the r50’s of the transitions in total substrate concentrations for these two proteins will be coupled. In other words, overexpression of one substrate could shift the total saturation, and thus the r50, of all the substrates coupled to that E3/DUB pair.

thumbnail
Fig 4. Effects of protein overexpression on total protein in single-substrate and multiple-substrate models.

(A) There is a non-linear transition in the [S]T vs. Q curve when Q is comparable in magnitude with the measure of saturation (in this case, [M]). The slope of the curve approaches 1 on either side of the transition, which is consistent with the maximum and minimum values of [S]T. (B) Schematic diagram for multiple substrates with one modification state each. Shown here is the model corresponding to two substrates, for simplicity. Here “M” denotes modifying enzyme and “D” denotes demodifying enzyme. Modified substrate is indicated by the dark, filled circle. (C) Plot of total concentration of first substrate vs. synthesis rate of the second substrate, for the multiple-substrate analogs of the Full model and the representative multiple modification state model. Axes are on a log scale. (D) Sensitivity to signal for first substrate vs. synthesis rate of the second substrate, for the multiple-substrate analogue of the Full model. The semi-analytical curve was obtained by substituting values for obtained empirically from simulation into α2 in the analytical expression for r50 ([S1]T). Axes are on a log scale.

https://doi.org/10.1371/journal.pcbi.1008492.g004

To test this, we introduced more than one substrate in the context of the full model. For the sake of display, we have taken the total number of substrates N in our model to be 2. As shown in Fig 4B, each E3 ligase and DUB now acts on two substrates with one modification state per substrate. The set of ODEs describing the model is given in S1 Text, Sec. 3. In contrast to the case of just one substrate, here we are interested in capturing the response of to changes in the synthesis rate of the second substrate, denoted by Q2. In Fig 4C, we see that increasing Q2 yields an increase in [S1]T for the Full model, as expected. The general behavior is similar in the Processive E3 and Distributive/Sequential DUB model, with the transition in [S1]T occurring at lower values of Q2.

Interestingly, r50([S1]T) also depends on Q2 (Fig 4D). Specifically, when Q2 is large enough in the Full model, r50 increases linearly with respect to Q2. In a similar manner to Fig 4C, a lower value of Q2 is sufficient to obtain a linear increase in r50 in the presence of multiple modification states. Note the excellent agreement in the Full model between the r50 values extracted empirically from simulation output and the analytical expression for r50([S1]T) (Fig 4D and S1 Text, Sec. 3.3). In fact, when we make Q2 arbitrarily large, we obtain: where α2 is the molar fraction of modified S2 at steady-state when [S1]T is half-maximal. In other words, as the concentration of the second substrate is increased, more and more activity of the E3 ligase (i.e. M) is needed to drive the transition in S1 concentration. Interestingly, all of these results can be readily generalized for any number of substrates, independent of substrate identity (S1 Text, Sec. 3.2–3.3). Thus, crosstalk in PTMs can lead to coupling of not only modification states [18,19], but also of overall protein levels.

Discussion

It has been over 35 years since Goldbeter and Koshland discovered the phenomenon of 0th-order ultrasensitivity. Since then, there has been extensive characterization of PTM cycles with 0th-order ultrasensitivity, both experimentally [4547] and computationally [18,20,25,48,49]. Until now, however, the properties of PTM cycles that drive protein degradation have not been studied in a systematic way. Using a mathematical modeling framework, we found that adding synthesis and degradation to a PTM cycle suppresses both sensitivity to signal and ultrasensitivity of the response, even when the PTM in question does not serve as a signal for protein degradation (Fig 5). Thus switch-like behaviors in vivo may or may not be the consequence of 0th-order ultrasensitivity, depending on the stability of the protein substrate. Although there are exceptions [1517], most models of signaling networks ignore protein turnover [50,51]. Our findings indicate that incorporating turnover, especially turnover based on actual protein stabilities, is key to capturing the global PTM dynamics of signaling systems.

thumbnail
Fig 5. A summary of the key features and results for each model considered here.

https://doi.org/10.1371/journal.pcbi.1008492.g005

Interestingly, we found the general trend of decreasing sensitivity and ultrasensitivity holds for PTMs that drive protein degradation, even when accounting for many of the complicated mechanisms that describe polyubiquitylation by E3 ligases and deubiquitylation by DUB enzymes (Fig 5 and S1 Text, Sec. 2). By adding E3 ligase crosstalk, we demonstrated that overexpressing one protein can elevate the concentration of another, and can also reduce the sensitivity of other proteins to incoming signals that would drive their degradation (Fig 4C and 4D). In other words, if one protein is overexpressed, it becomes more difficult to degrade any of its counterparts sharing the same E3/DUB enzyme pair.

Although there is some data available about the specificity of E3 ligases [34,52,53], this information is very far from complete. Consider the highly common experimental scenario where a primary aim is to characterize the function of a protein by manipulating its expression level. Our findings indicate that the interpretation of overexpression data in eukaryotic cells may be very difficult because some of the observed phenotypic or molecular effects could be directly due to the higher concentration of the protein that was expressed, but other effects could be due to E3 ligase coupling (Fig 4C). Additional complications could appear due to the change in sensitivity to the shared E3 ligases for other substrates in the system (Fig 4D). For instance, if a protein is being actively regulated by its E3 ligase and a degradation signal appears, then a high concentration of other proteins in the system would potentially inhibit the signal. This could have unforeseen large-scale effects on the overall system.

A global picture of E3-ligase/DUB enzyme specificity will thus likely be essential to comprehending the regulation of protein levels within cells. This will allow us to begin determining how to isolate direct effects of changes in protein expression levels from indirect effects. Equally necessary are mathematical or computational models of signaling dynamics, gene regulatory networks, and other cellular processes that describe the interplay between PTMs that do not lead to degradation and those that drive degradation. Incorporating the coupled dynamics of protein levels into our understanding of cell signaling and cellular physiology thus represents a grand challenge for both experimental and computational systems biology.

Materials and methods

Experimental methods

Our model behaviors can be described deterministically by systems of ordinary differential equations (ODEs). Numerical integration of the systems was performed by the stiff solver ode15s in MATLAB. All analyses were performed at steady-state. In parallel, agent-based stochastic simulations of the systems [3941] were conducted using custom-built software implemented in C++. Parameter values were chosen to ensure equivalence between the deterministic and stochastic systems. See the supporting information for full details regarding all of the models considered here.

Supporting information

S1 Text. This appendix contains further details about the models studied here, as well as additional mathematical and numerical results.

https://doi.org/10.1371/journal.pcbi.1008492.s001

(PDF)

Acknowledgments

We thank William Mather and Christian Ray for helpful discussions regarding this work.

References

  1. 1. Gasparian AV, Yao YJ, Kowalczyk D, Lyakh LA, Karseladze A, Slaga TJ, et al. The role of IKK in constitutive activation of NF-kappaB transcription factor in prostate carcinoma cells. J Cell Sci. 2002 Jan 1;115(Pt 1):141–51. pmid:11801732
  2. 2. Calamini B, Morimoto RI. Protein homeostasis as a therapeutic target for diseases of protein conformation. Curr Top Med Chem. 2012;12(22):2623–40. pmid:23339312
  3. 3. Powers ET, Morimoto RI, Dillin A, Kelly JW, Balch WE. Biological and chemical approaches to diseases of proteostasis deficiency. Annu Rev Biochem. 2009;78:959–91. pmid:19298183
  4. 4. Finley D. Recognition and processing of ubiquitin-protein conjugates by the proteasome. Annu Rev Biochem. 2009;78:477–513. pmid:19489727
  5. 5. Komander D, Rape M. The ubiquitin code. Annu Rev Biochem. 2012;81:203–29. pmid:22524316
  6. 6. Goldbeter A, Koshland DE. An amplified sensitivity arising from covalent modification in biological systems. Proc Natl Acad Sci U S A. 1981 Nov;78(11):6840–4. pmid:6947258
  7. 7. Goldbeter A, Koshland DE. Ultrasensitivity in biochemical systems controlled by covalent modification. Interplay between zero-order and multistep effects. J Biol Chem. 1984 Dec 10;259(23):14441–7. pmid:6501300
  8. 8. Xing J, Chen J. The Goldbeter-Koshland switch in the first-order region and its response to dynamic disorder. PLoS One. 2008 May 14;3(5):e2140. pmid:18478088
  9. 9. Xu Y, Gunawardena J. Realistic enzymology for post-translational modification: zero-order ultrasensitivity revisited. J Theor Biol. 2012 Oct 21;311:139–52. pmid:22828569
  10. 10. Martins BM, Swain PS. Ultrasensitivity in phosphorylation-dephosphorylation cycles with little substrate. PLoS Comput Biol. 2013;9(8):e1003175. pmid:23950701
  11. 11. Zhang Q, Bhattacharya S, Andersen ME. Ultrasensitive response motifs: basic amplifiers in molecular signalling networks. Open Biol. 2013 Apr 24;3(4):130031. pmid:23615029
  12. 12. Ferrell JE, Ha SH. Ultrasensitivity part I: Michaelian responses and zero-order ultrasensitivity. Trends Biochem Sci. 2014 Oct;39(10):496–503. pmid:25240485
  13. 13. Nguyen LK, Dobrzyński M, Fey D, Kholodenko BN. Polyubiquitin chain assembly and organization determine the dynamics of protein activation and degradation. Front Physiol. 2014;5:4. pmid:24478717
  14. 14. Nguyen LK, Muñoz-García J, Maccario H, Ciechanover A, Kolch W, Kholodenko BN. Switches, excitable responses and oscillations in the Ring1B/Bmi1 ubiquitination system. PLoS Comput Biol. 2011 Dec;7(12):e1002317. pmid:22194680
  15. 15. Loriaux PM, Hoffmann A. A framework for modeling the relationship between cellular steady-state and stimulus-responsiveness. Methods Cell Biol. 2012;110:81–109. pmid:22482946
  16. 16. Loriaux PM, Hoffmann A. A protein turnover signaling motif controls the stimulus-sensitivity of stress response pathways. PLoS Comput Biol. 2013;9(2):e1002932. pmid:23468615
  17. 17. Loriaux PM, Tesler G, Hoffmann A. Characterizing the relationship between steady state and response using analytical expressions for the steady states of mass action models. PLoS Comput Biol. 2013;9(2):e1002901. pmid:23509437
  18. 18. Rowland MA, Fontana W, Deeds EJ. Crosstalk and competition in signaling networks. Biophys J. 2012 Dec 5;103(11):2389–98. pmid:23283238
  19. 19. Rowland MA, Deeds EJ. Crosstalk and the evolution of specificity in two-component signaling. Proc Natl Acad Sci U S A. 2014 Apr 15;111(15):5550–5. pmid:24706803
  20. 20. Rowland MA, Harrison B, Deeds EJ. Phosphatase specificity and pathway insulation in signaling networks. Biophys J. 2015 Feb 17;108(4):986–96. pmid:25692603
  21. 21. Ogle CT, Mather WH. Proteolytic crosstalk in multi-protease networks. Phys Biol. 2016 Apr 4;13(2):025002. pmid:27042892
  22. 22. Cookson NA, Mather WH, Danino T, Mondragón-Palomino O, Williams RJ, Tsimring LS, et al. Queueing up for enzymatic processing: correlated signaling through coupled degradation. Mol Syst Biol. 2011 Dec 20;7:561. pmid:22186735
  23. 23. Han Y, Lee H, Park JC, Yi GS. E3Net: a system for exploring E3-mediated regulatory networks of cellular functions. Mol Cell Proteomics. 2012 Apr;11(4):O111.014076. pmid:22199232
  24. 24. Du Y, Xu N, Lu M, Li T. hUbiquitome: a database of experimentally verified ubiquitination cascades in humans. Database (Oxford). 2011;2011:bar055. pmid:22134927
  25. 25. Gomez-Uribe C, Verghese GC, Mirny LA. Operating regimes of signaling cycles: statics, dynamics, and noise filtering. PLoS Comput Biol. 2007 Dec;3(12):e246. pmid:18159939
  26. 26. Scheer M, Grote A, Chang A, Schomburg I, Munaretto C, Rother M, et al. BRENDA, the enzyme information system in 2011. Nucleic Acids Res. 2011 Jan;39(Database issue):D670–6. pmid:21062828
  27. 27. Eden E, Geva-Zatorsky N, Issaeva I, Cohen A, Dekel E, Danon T, et al. Proteome half-life dynamics in living human cells. Science. 2011 Feb 11;331(6018):764–8. pmid:21233346
  28. 28. Cambridge SB, Gnad F, Nguyen C, Bermejo JL, Krüger M, Mann M. Systems-wide proteomic analysis in mammalian cells reveals conserved, functional protein turnover. J Proteome Res. 2011 Dec 2;10(12):5275–84. pmid:22050367
  29. 29. Dubitzky W, Wolkenhauer O, Cho K, Yokota H. Encyclopedia of Systems Biology. Springer Publishing Company, Incorporated; 2013.
  30. 30. Proctor CJ, Tsirigotis M, Gray DA. An in silico model of the ubiquitin-proteasome system that incorporates normal homeostasis and age-related decline. BMC Syst Biol. 2007 Mar 21;1:17. pmid:17408507
  31. 31. Kleiger G, Saha A, Lewis S, Kuhlman B, Deshaies RJ. Rapid E2-E3 assembly and disassembly enable processive ubiquitylation of cullin-RING ubiquitin ligase substrates. Cell. 2009 Nov 25;139(5):957–68. pmid:19945379
  32. 32. Lee MJ, Lee BH, Hanna J, King RW, Finley D. Trimming of ubiquitin chains by proteasome-associated deubiquitinating enzymes. Mol Cell Proteomics. 2011 May;10(5):R110.003871. pmid:20823120
  33. 33. Livneh I, Cohen-Kaplan V, Cohen-Rosenzweig C, Avni N, Ciechanover A. The life cycle of the 26S proteasome: from birth, through regulation and function, and onto its death. Cell Res. 2016 8;26(8):869–85. pmid:27444871
  34. 34. Lu Y, Lee BH, King RW, Finley D, Kirschner MW. Substrate degradation by the proteasome: a single-molecule kinetic analysis. Science. 2015 Apr 10;348(6231):1250834. pmid:25859050
  35. 35. Metzger MB, Hristova VA, Weissman AM. HECT and RING finger families of E3 ubiquitin ligases at a glance. J Cell Sci. 2012 Feb 1;125(Pt 3):531–7. pmid:22389392
  36. 36. Metzger MB, Pruneda JN, Klevit RE, Weissman AM. RING-type E3 ligases: master manipulators of E2 ubiquitin-conjugating enzymes and ubiquitination. Biochim Biophys Acta. 2014 Jan;1843(1):47–60. pmid:23747565
  37. 37. Pierce NW, Kleiger G, Shan SO, Deshaies RJ. Detection of sequential polyubiquitylation on a millisecond timescale. Nature. 2009 Dec 3;462(7273):615–9. pmid:19956254
  38. 38. Ventii KH, Wilkinson KD. Protein partners of deubiquitinating enzymes. Biochem J. 2008 Sep 1;414(2):161–75. pmid:18687060
  39. 39. Feret J, Danos V, Krivine J, Harmer R, Fontana W. Internal coarse-graining of molecular systems. Proc Natl Acad Sci U S A. 2009 Apr 21;106(16):6453–8. pmid:19346467
  40. 40. Deeds EJ, Krivine J, Feret J, Danos V, Fontana W. Combinatorial complexity and compositional drift in protein interaction networks. PLoS One. 2012;7(3):e32032. pmid:22412851
  41. 41. Chylek LA, Harris LA, Tung CS, Faeder JR, Lopez CF, Hlavacek WS. Rule-based modeling: a computational approach for studying biomolecular site dynamics in cell signaling systems. Wiley Interdiscip Rev Syst Biol Med. 2014 Jan-Feb;6(1):13–36. pmid:24123887
  42. 42. Nariya MK, Israeli J, Shi JJ, Deeds EJ. Mathematical Model for Length Control by the Timing of Substrate Switching in the Type III Secretion System. PLoS Comput Biol. 2016 Apr;12(4):e1004851. pmid:27078235
  43. 43. Markevich NI, Hoek JB, Kholodenko BN. Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J Cell Biol. 2004 Feb 2;164(3):353–9. pmid:14744999
  44. 44. Takahashi K, Tanase-Nicola S, ten Wolde PR. Spatio-temporal correlations can drastically change the response of a MAPK pathway. Proc Natl Acad Sci U S A. 2010 Feb 9;107(6):2473–8. pmid:20133748
  45. 45. Hardie DG, Salt IP, Hawley SA, Davies SP. AMP-activated protein kinase: an ultrasensitive system for monitoring cellular energy charge. Biochem J. 1999 Mar 15;338 (Pt 3):717–22. pmid:10051444
  46. 46. LaPorte DC, Koshland DE. Phosphorylation of isocitrate dehydrogenase as a demonstration of enhanced sensitivity in covalent regulation. Nature. 1983 Sep 22–28;305(5932):286–90. pmid:6312317
  47. 47. Meinke MH, Bishop JS, Edstrom RD. Zero-order ultrasensitivity in the regulation of glycogen phosphorylase. Proc Natl Acad Sci U S A. 1986 May;83(9):2865–8. pmid:3458247
  48. 48. Gunawardena J. Multisite protein phosphorylation makes a good threshold but can be a poor switch. Proc Natl Acad Sci U S A. 2005 Oct 11;102(41):14617–22. pmid:16195377
  49. 49. Thomson M, Gunawardena J. Unlimited multistability in multisite phosphorylation systems. Nature. 2009 Jul 9;460(7252):274–7. pmid:19536158
  50. 50. Albeck JG, Burke JM, Spencer SL, Lauffenburger DA, Sorger PK. Modeling a snap-action, variable-delay switch controlling extrinsic cell death. PLoS Biol. 2008 Dec 2;6(12):2831–52. pmid:19053173
  51. 51. Suderman R, Deeds EJ. Machines vs. ensembles: effective MAPK signaling through heterogeneous sets of protein complexes. PLoS Comput Biol. 2013;9(10):e1003278. pmid:24130475
  52. 52. Merbl Y, Refour P, Patel H, Springer M, Kirschner MW. Profiling of ubiquitin-like modifications reveals features of mitotic control. Cell. 2013 Feb 28;152(5):1160–72. pmid:23452859
  53. 53. Merbl Y, Kirschner MW. Large-scale detection of ubiquitination substrates using cell extracts and protein microarrays. Proc Natl Acad Sci U S A. 2009 Feb 24;106(8):2543–8. pmid:19181856