Crosstalk and ultrasensitivity in protein degradation pathways

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.

We begin with a modification of the Goldbeter-Koshland loop that incorporates protein turnover. It is described by the following scheme of enzymatic reactions: Here M and D represent any modifying and demodifying enzyme, respectively. Any substrate in the S state is degraded at a first-order rate δ 1 , and any in the S * state at a rate δ 2 , with δ 2 > δ 1 . Substrate is also synthesized at a constant rate Q, and all the synthesized substrates are in the S (unmodified) state. Note that degradation does not consume M or D.
(We find that relaxing this assumption only changes Equations (5) and (6), the ODEs that correspond to the concentration of the free M and D enzyme. Interestingly, these two equations are not used anywhere in the derivation. As a result, relaxing this assumption has no effect on our findings.) It is straightforward to use the Law of Mass Action to formulate the corresponding set of ordinary differential equations (ODEs), as described below.

Model definitions
Using the formalism above, we can describe the Goldbeter-Koshland model as having no synthesis or degradation (i.e. with Q = δ 2 = δ 1 = 0). The Intermediate model has a positive synthesis rate and a uniform rate of degradation (i.e. with Q > 0 and δ 2 = δ 1 > 0.) Finally, the Full model also has a positive synthesis rate, but has different rates of degradation (i.e. Q > 0 and δ 2 > δ 1 > 0.)
For purposes of display, the case of equal K M s (i.e. K M,M = K M,D ) is analyzed below. Analyses of scenarios with substantially different K M s are left to future work. It follows that K M,1 ≈ K M,2 because k +,M , k +,D δ 1 , δ 2 for the majority of enzymes.

Parameter Values
The following figure illustrates the biological relevance of our choices of parameter values. The curve in each plot is a kernel density estimate of the experimental values obtained from the BRENDA enzyme database (1). Figure S1: Distributions of relevant steady-state kinetic parameters. (A) Logarithmicscale density plots of the K M s for protein kinases, protein phosphatases, E3 ligases, and DUB enzymes in the human genome. Experimental values were obtained from the BRENDA enzyme database (1) by using EC numbers for the tyrosine kinases/phosphatases and serine/threonine kinases/phosphatases. (B) Identical procedure as in the previous panel, implemented for k cat s.

Analytical expressions for r 50
Let r denote the ratio of maximum velocities of the two enzymes (i.e. r ≡ . This ratio represents the signal source for the system as defined previously (2). Let α ≡ [S * ]/[S] T , or the molar fraction of modified substrate at steady-state. The r 50 of the response is defined as the amount of r necessary to yield a 50% response in modified substrate, or α = 0.5.
We can rewrite eq.

Intermediate Model
For the intermediate model, since Q, δ 1 > 0, note that eq. (10) becomes: We can substitute α = 0.5 in the expression above to obtain r 50 . Upon simplification, we obtain

Full Model
For the full model, we directly solve for [S], [S * ] in terms of α. Observe that .
Plugging these expressions into eq. (10), we arrive at a closed-form expression for r as a function of α: Substituting α = 0.5 in eq. (13) and simplifying the resulting expression yields:

Analytical expressions for n ef f
The effective Hill coefficient n ef f is defined in (4) as log(81)/ log( EC 90 EC 10 ). From (5), we have the expression for n ef f in the Goldbeter-Koshland model: We can substitute α = 0.1, 0.9 in both eq. (11) and eq. (13) to yield both EC 10 and EC 90 for the intermediate and full models, respectively. Simplifying, we obtain: where x = 9δ 1 + δ 2 and y = 9δ 2 + δ 1 .

Analysis of r 50 in saturated regimes
We show that r 50 (GK) < r 50 (Intermediate) < r 50 (Full) when total substrate is at saturating concentrations (i.e. when [S] T K M ). First observe that r 50 (GK) = 1 in a saturated regime. This can be seen by setting δ 2 = 0 in eq. (10) and noting that [S], [S * ] K M : Furthermore, we can show that different modes of enzyme saturation (i.e. increasing Q vs. decreasing K M ) effect substrate responses in distinct ways. Specifically, varying Q has a stronger effect on r 50 than varying K M -this is true in both the Intermediate and Full models. We use the following chain of (unitless) inequalities to order the quantities of interest: The only condition on this chain is Q δ 1 > 2K M , which is trivially satisfied in the saturated regime (i.e. [S] T K M ).

Analysis of n ef f
When total substrate is at saturating concentrations, since n ef f (Intermediate) has an additional positive term in the denominator compared to n ef f (GK), n ef f (Intermediate) < n ef f (GK).
Although it would be desirable to compare n ef f (GK) with n ef f (Full), this is hard to do for two reasons. However, we can show that in the limit K M → 0, n ef f (Intermediate) is a strictly decreasing function in Q: has a negative partial derivative with respect to Q. Entering the expression into Mathematica (3) and simplifying yields: We can also show that n ef f (Intermediate) > 1 for all possible values of the parameters in the model. This is equivalent to showing is not constant with respect to r.
The half-maximal response for total substrate occurs at 1 2 (min [S] T + max [S] T ). Since we can also solve for [S] T in terms of [S]: From eq. (15), we see that We can plug in and solve for α in eq. (13) and simplifying the resulting expression yields:

Analytical expression for n ef f of [S] T in Full model
To obtain an expression for n ef f , we need to first derive EC 10 and EC 90 . Note that the 10% response for total substrate occurs at 9 10 min [S] T + 1 10 max [S] T and the 90% response for total substrate occurs at 1 10 min [S] T + 9 10 max [S] T .
To derive EC 10 (β), we can plug in and solve for α when β = 9 10 min [S] T + 1 10 max [S] T : in eq. (13) and simplifying the resulting expression (3) yields: Similar calculations can be done to derive EC 90 : Thus n ef f = log(81)/ log

Single Substrate, Multiple Modification States
We generalize the model in the previous section to various models with arbitrarily long chains of ubiquitin units. In these models, substrates with zero to three ubiquitin units attached are degraded at the first-order rate δ 1 and substrates with at least four units are degraded at the rate δ 2 > δ 1 . In what follows, the maximal length of the chain is denoted by . The indices i and j represent the number of the ubiquitin unit and range from 0 to 3 and 4 to respectively. The index x ranges from 1 to 3.

Parameter Values
The values for parameters in the table below were chosen to be consistent with those of the single substrate models. The parameters listed here correspond to the deterministic version of the models. However, [S] T is one order of magnitude smaller. Thus, the corresponding values for Q, k +,M and k +,D are also scaled accordingly (i.e. Q is scaled down by a factor of ten, while k +,M and k +,D are scaled up by the same factor). Our motivation for this choice of [S] T was to lower the computational cost associated with the stochastic agent-based simulations, effectively decreasing the total number of agents in these simulations (see Section 2.8).

Parameter Values
The same reasoning applies here as in Section 2.

Parameter Values
In this model version, the changes from the previous models are reflected in a lower association rate k +,D and higher dissociation rate k −,D for the D enzyme. This parameterization is consistent with both the mechanism of D as well as the experimentally observed values in (6). . . .

Parameter Values
This model variant has a parameterization similar to that of the model with Distributive E3 & Processive/Sequential DUB. However, the mechanisms of the E3 ligase and DUB enzyme are swapped here. Also, we have now introduced the parameters k cat,1,M and k cat,2,M in place of k cat,M . Since the chemistry of an E3 ligase binding a ubiquitin unit to a protein is fundamentally different from that of an E3 ligase binding a unit to another ubiquitin unit, our parameterization is consistent with both the mechanism of E3 as well as the experimentally observed values in (6).  (2) . . . (1) . . .  (2) . . . (1) . . .

Graphical results
In the following figures, in order to choose a reasonable value for , we first ran simulations such that changes in the r 50 and n ef f were negligible beyond a point. Using these results, we then chose = 50 by inspection. The full range of the transitions in the curves are visible for each model. The continuous lines indicate the numerically integrated deterministic solutions.
Note that the case corresponding to the Processive E3 & Distributive/Sequential DUB model has been presented in the Main Text.

Stochastic simulations
Note that the definition of a maximum length above is necessary in order to ensure a finite set of chemical reactions and ODEs. To determine if this truncation has any effect on the results, we compared our deterministic case to stochastic simulations in which we allow the ubiquitin chains to reach an arbitrary length. Our approach to developing these simulations is inspired by "agent-based" simulators developed for the stochastic simulation of rule-based models (7; 8). Due to technical considerations with the mechanism, however, we wrote our own dedicated code for these simulations, following closely the approach taken in our previous work on modeling length control in the bacterial Type III Secretion System (9).
Briefly, our simulations contain three types of agents: the M enzyme, the D enzyme, and the substrate S. These agents are represented independently; in other words, if there are 1000 S molecules in the simulation, this is represented by having 1000 distinct "S" agents in memory. Each S agent has associated with it a number that represents the length of its ubiquitin chain. These lengths can range from zero to the largest number that can be represented by the particular data structure. Since this number is much, much larger than the largest value ever practically observed in the simulation, this essentially corresponds to allowing for arbitrary chain lengths.
We wrote a separate simulation in C++ for all of the scenarios described above. All of the parameters from the deterministic simulations were converted to their corresponding stochastic values in a straightforward way (10). In particular, we specified a "compartment" volume of 100 fL for each model. Since these simulations are relatively expensive computationally, we performed simulations for a subset of the parameters considered in our deterministic simulations (see below). All simulations were performed until the system achieved a steady state. Simulation codes are available upon request.

No truncation effect
In the following set of figures, the dots correspond to averages from stochastic simulations. We conducted 100 simulations for each model.
We note that, within the parameter regimes we were able to simulate in a reasonable amount of time, there is excellent agreement between the deterministic and stochastic frameworks. In other words, there does not seem to be an effect of premature truncation in chain length for the deterministic model versions.

Robustness of results
In order to explore the robustness of our conclusions with respect to parameter variation, we randomly sampled a reasonable region of parameter space for each of our multiple modification state models and characterized how ODE models based on those parameters behaved at steady state. We focused on sampling different values of K M for the enzymes and different degradation rates. To sample K M values, we fixed the values of k cat and k − at 0.999 s −1 and 0.001 s −1 , respectively, and sampled values of k + . We chose to do this because all of our analysis is focused on steady-state responses, and so the specific values of the rate constant only influence our results through the ratio K M = (k cat + k − )/k + . Setting k cat and k − to these values allows us to focus our analysis on the case of catalytically efficient enzymes. To sample values of K M , we sampled k + from a log-uniform distribution ranging from 0.1× to 10× the K M we originally chose for our analysis. We did this independently for both the M and D enzymes in the model. Similarly, we sampled δ 1 from a log-uniform distribution between 0.1× and 10× the original δ 1 . We fixed δ 2 to equal 10 × δ 1 , since our model focuses on the case where the ubiquitin chain drives an increased rate of substrate degradation. Finally, two values of the Q parameter were chosen to correspond to the unsaturated and saturated regimes.
One of the key findings of our work is the fact that increasing saturation by increasing the production of substrate (our parameter Q) increases the ratio of E3 to DUB activity that is needed to see a transition in substrate concentration (i.e. it increases the r 50 of the transition). To see if this held for different parameter sets, we computed the ratio r 50,saturated /r 50,unsaturated for each of the parameter sets; values of this ratio greater than 1 indicate that saturation increases r 50 . A histogram of this ratio over all our randomly sampled parameters indicates that they all indeed have an r 50 ratio greater than 1.
Another major finding of our work is that the ultrasensitivity of the response (for both the molar fraction of modified substrate, [S * ]/[S] T and the total amount of substrate, [S] T ) is much less than predicted by the standard Goldbeter-Koshland model of PTM cycles. To characterize this, we calculated the ratio of the effective Hill coefficients, n ef f for the saturated vs. the unsaturated cases. We see that this ratio is generally around 1 to 2, and never much larger than 5, indicating that the large increase in ultrasensitivity observed for standard saturated PTM cycles is not observed for these parameter sets.
Thus, although we cannot exhaustively explore parameter space, these new results indicate that our main conclusions are likely robust to reasonable variations of the parameters.

Comments
Adding multiple substrates to the system results in additive effects. Similar conclusions hold for multiple substrates with just one modification state as for one substrate with one modification state (Section 1). The only changes in the relevant expressions are: (i) The substrate production rate Q becomes z Q z (ii) S becomes z S z (iii) S * becomes z S * By inspecting eq. (13) and its derivation, one can see that the expression for r in this case will be similar. Following the derivation in Section 1.8, we can obtain r 50 ([S 1 ] T ) as follows: