Quantitative analysis reveals how EGFR activation and downregulation are coupled in normal but not in cancer cells

Ubiquitination of the epidermal growth factor receptor (EGFR) that occurs when Cbl and Grb2 bind to three phosphotyrosine residues (pY1045, pY1068 and pY1086) on the receptor displays a sharp threshold effect as a function of EGF concentration. Here we use a simple modelling approach together with experiments to show that the establishment of the threshold requires both the multiplicity of binding sites and cooperative binding of Cbl and Grb2 to the EGFR. While the threshold is remarkably robust, a more sophisticated model predicted that it could be modulated as a function of EGFR levels on the cell surface. We confirmed experimentally that the system has evolved to perform optimally at physiological levels of EGFR. As a consequence, this system displays an intrinsic weakness that causes—at the supraphysiological levels of receptor and/or ligand associated with cancer—uncoupling of the mechanisms leading to signalling through phosphorylation and attenuation through ubiquitination.


Lysates of HeLa cells stably transfected with
To avoid using different concentrations for the same species, we computed concentrations using the volume of the medium occupied by each cell !"#$%! , as shown in 6 . Notice that for this reason, the concentrations are not directly comparable to the experimental one: for the conversion, it will be enough to multiply for the proper volume. For example, the concentration of EGFR in the medium was computed from the concentration in the membrane as !"#$%! = !"!#$%&" × !"!#$%&" / !"#$%! . In a typical experiment, we have in the order of 3-5x10 6 cells growing in 3-5 ml of medium, which gives the volume per cell !"#$%! ≈ 0.5×10 !! ml. To estimate !"# , we approximated HeLa cells as spheres having a radius of 15 . We thus obtained !"# ≈ 10 !! ml.
!"!#$%&" was estimated as the surface of the cell x the height of EGFR, which is of the order of 10 nm, resulting in !"!#$%&" ≈ 10 !!! ml. The rescaling affected also the constants with an explicit dependence on concentration. Also in this case, as shown above, rescaling took into account the ratio of the volumes. Thus, for all the parameters we provided values in typical units (to be directly compared to the one described in the literature) and the corresponding values scaled to the volume of the medium (used to build the model).
Relations linking equilibrium constants to the corresponding forward and backward rate The parameters not listed can be computed using the definitions given above, e.g. ! = ! × .
a. Assuming that HeLa cells contain an average of 2.5--3×10 ! EGFRs (measured by saturation binding assay, see Methods) and that each cell is surrounded by 0.5×10 !! ml of medium. b. We estimated the amount of Cbl available for EGF interaction by measuring the levels of Cbl phosphorylated by the EGFR at the maximal EGF concentration (see main text). This represents the maximal "active Cbl" fraction, i.e. able to recruit E2 ubiquitin-conjugating enzymes and to ubiquitinate the EGFR. Thus, we take this as a proxy for the maximal amount of Cbl actually available for binding to the EGFR (measured in Fig. 3b). c. Measured in Fig. 3b

Supplementary note 1 -EGFR models
In this section, we provide synopses of the EGFR.
The Multi-site Phosphorylation Model (MPM) describes EGFR phosphorylation and dephosphorylation using the ordinary differential equations (ODE) provided in Eq. (2). The mapping between the amount of EGF and the extent of kinase activity of EGFR is modeled with a Hill function, see Eq.
where ! denotes an EGFR with α phosphate groups, while !"# and !"! denote the rate constants for phosphorylation and dephosphorylation of ! , respectively. The chemical reactions of Eq.
Free-free and saturated-saturated regimes studied with the Gillespie method. In Fig. 2b, we show that two limiting regimes (free and saturated) can be identified for the reaction catalyzed by kinases, depending on the stability of the complex between the catalytic subunit of EGFR and its target Tyr. To simulate the reactions defined by Eq. (2), since we are dealing with a small number of molecules-one enzyme and a maximum of 9 substrates (Tyr) for EGFR-WT-we have to resort to stochastic simulations. We fixed !"# = 10 and !"" = 1 and varied !" . in units of !"# (the basic rate constant that governs the addition of a phosphoryl group in the absence of competing Tyr). For !" + !"" / !"# = 10 !! (blue dashed line), !"# controls the ratelimiting step, a condition that we call the saturated-enzyme regime. For !" + !"" / !"# = 10 ! (black line), !" controls the rate limiting step, a condition that we dubbed the free-enzyme regime. In summary, the figure show that the free-enzyme regime corresponds to !"# = − !"# for kinases and !"! = !"! for phosphatases, while the saturated-enzyme regime corresponds to !"# = !"# for kinases and !"! = !"! for phosphatases, where !"# and !"! denote the phosphorylation and dephosphorylation rate constants for a single Tyr. Since each regime may be applied independently to kinases and phosphatases, in principle we have four models of phosphorylation: the saturated-saturated, saturated-free, free-saturated and free-free where the first term always refers to kinases and the second to phosphatases. Only when the freefree regime is adopted (not shown), does the model reproduce faithfully the shape of the doseresponse curves, as well as the ratio of EGFR-WT and EGFR-3Y+ phosphorylation ( Fig. 2d-e).
Since EGFR follows the free-enzyme regime for both phosphorylation and dephosphorylation, the phosphorylation of each Tyr is independent of that of the other Tyr. As a consequence, modeling EGFR with = 9 or = 3 does not make any difference, because only Y1045, Y1068 and Y1086 are relevant for EGFR ubiquitination.
Mapping how EGFR activity depends on EGF concentration. We then needed to define how EGF modifies the extent of EGFR activity, i.e. to define !"# as a function of EGF. We disregarded altogether the mechanisms of EGFR activation and followed an empirical approach: we assumed that !"# is a function of EGF as described by a Hill function: In Eq. (3), L stands for the EGF concentration, the Hill coefficient represents how steeply !"# increases with , while is the EGF concentration where !"# is half-maximal. Eq. (3) describes a sigmoid varying between 0 and !"# !"# , which is the maximal kinetic activity of EGFRs.
Conversely, we assume that !"! stays constant for EGFR in the plasma membrane in the first two minutes after EGF stimulus. With this mapping, the equations describing the transitions among the phosphorylation states of EGFR are those described by Eq. (1) and (2), provided that !"# is substituted by !"# . In summary, the occupancy of each of the ! , which denotes an EGFR with α phosphate groups, depends on the concentration of EGF and on the parameters ! , , !"# !"# , and !"! , collectively represented by a vector . For a given EGF concentration L, ! , , represents the fraction of EGFR with phosphate groups at time t.
Measuring phosphorylation in MPM. Since the C-terminal tail of EGFR, where the relevant Tyr are located, is predicted to be unstructured and very flexible, we assumed that the pan-pY antibody (4G10, Millipore) recognizes phosphorylated Tyr independently of EGFR conformation. -! , , -which depends on time, EGF concentration L, and on the parameters ! , , !"# !"# , and !"! -has a number of pYs equal to , with = 0, … , . The expression for the total level of EGFR phosphorylation that can be directly compared to the experimental data is: EGFR binding to Cbl and Grb2. In this section, we discuss the mathematical formalization of EGFR ubiquitination, which depends on Cbl and Grb2 binding to the receptor. As summarized in  Having defined the new variables, we will now list all possible reactions including Cbl and Grb2. Note that despite the large number of reactions involved, they are all straightforward binding reactions. It is worth stressing however three peculiarities. First, if EGFR has both pY1068 and pY1086, the Grb2-binding rate doubles with respect to a form of EGFR having only pY1068 or pY1086 (doubling the binding sites doubles the chance to bind). Second, binding reactions are independent of one another. Third, we introduced a cooperative mechanism as explained in the main text, Fig. 3a.
The list of reactions begins with Cbl and Grb2 binding to free receptors and to each other. To shorten the list, we used an additional notation, without changing the number of species. To identify the Y1045 site, we use the symbol Z that can take two values: 'empty' and 'P'. Where 'empty' stands for Y1045 and 'P' for pY1045. For the Y1068/Y1086 site, X can take the values 'empty', 'P' and '2P' (the latter two collectively represented as ' ', with = 1, 2). The association and dissociation reactions are where: 45 ( 1) and (2), respectively.
The reactions regarding the binding of a second Cbl or Grb2 moiety to receptors are and those for a third moiety  (7) Finally, the three reactions that characterize the cooperative hypothesis for ubiquitination are which lead to the formation of only fully ubiquitinated EGFR complexes.
The binding reactions of EGFR, Cbl, and Grb2 can take place in different ways, for example the binding of Cbl to pY1045 of EGFR can follow a first-order kinetics as described in reactions (8), with rate kb45*, but also a second-order kinetics, as in (5) and (6), with rate kb45. A similar observation can be made for the other reactions described in (8). The binding is favored by a cooperative mechanism when the first-order reactions are favored compared to the second-order reactions, i.e. when the following condition is satisfied: (the latter also implies kbcg * ≫ kbcg × Cbl ). Equations (9) which is several orders of magnitude higher than the typical concentrations of Cbl that we simulate (≈10 -3 nM). The simulations shown in Fig. 4, Fig. 5, Fig. 6b, Fig. 7c, Fig. 8a-c of the main text satisfy Eqs. (9), a condition that we define as cooperativity. If we disregard the spatial localization, i.e., we consider the single molecule as distributed in the cytosol with volume V cyt ≈ 10 −8 ml and Fig. 5a-b of the main text), we have a non-cooperative mechanism.
Phosphorylation/dephosphorylation in the presence of Cbl/Grb2 binding. It is likely that the phosphate groups involved in Grb2 or Cbl binding to the EGFR cannot be dephosphorylated when occupied. Thus, the dephosphorylation kinetics of the EGFR needs to be modified slightly to take into account Cbl and Grb2 binding. For example, when Cbl is bound to pY1045, it seems natural to assume that pY1045 is inaccessible to phosphatases and cannot be dephosphorylated.
Unfortunately, this constraint renders the phosphorylation kinetics dependent on the EGFR configuration and we cannot fully adopt the phosphorylation model where we do not distinguish between different Tyr. Otherwise, using the same example as before, we would have the spurious effect that a pY1045-bound Cbl protects also pY1068 from dephosphorylation. To avoid this inconsistency, we single out Y1045, while we still consider collectively Y1068 and Y1086. Rather then using variables ! , ! , ! and ! , we thus use the already introduced variables , ! , ! , !" , The total EGFR Tyr phosphorylation is expressed as To compute total ubiquitination, the relevant species are the fully bound trimer ! ! : : 2 and the complexes ! !:!"# , ! !:!"#:!"#$ , !!:!"#$ !:!"# and !!:!"#$ !:!"#:!"#$ where Cbl is only singly bound to the receptor. For brevity, again, we introduce two intermediate species: The ubiquitination then equals Eqs (14) show, in fact, that !"## and !" are Cbl-bound EGFR species rather than true ubiquitinated EGFRs. To be precise, they should be regarded as ubiquitination probabilities. To avoid introducing ad-hoc sources of sigmoidicity, we modeled actual ubiquitination as a completely first-order process, dR Ub dt = k ub Ub − k dub R Ub , which at steady state becomes R Ub = k ub Ub k dub so that, eventually, the normalized Ub is computed as where !"# is the maximum Ub value in the simulations. Note that we did not allow for saturation of receptor ubiquitination, because we assumed that EGFR can be ubiquitinated at multiple sites.
We used Eqs. (12) and (14) [or (15) at steady state that, for normalized ubiquitination, is equivalent to (14)] to compute the curves shown in Fig. 4, Fig. 5, Supplementary Fig. 3, Supplementary Fig. 4b-c, Supplementary Fig. 5 and Supplementary Fig. 6. Possible dimer forms are ! ∶ ! , ! ∶ ! , and ! ∶ ! , where , = 0, … , . Here, and in the following model, we do not distinguish between the position of monomers in the dimer, i.e., there is no 'left' or 'right' distinction. Thus, ! ∶ ! is equivalent to ! ∶ ! , ! ∶ ! is equivalent to ! ∶ ! , and ! ∶ ! is equivalent to ! ∶ ! . While in general ! ∶ ! and ! ∶ ! are two distinct species because the EGF-bound receptor is distinguishable by the number of phosphate groups attached, when = they represent the same species. To avoid duplication, we use the convention that the left index is always smaller than or equal to the right index, and we put the ligand on the right. Thus, ! ∶ ! only exists when < , while the species ! ∶ ! does not exist. The reader should remember these conventions when reading the ODEs that refer to dimers.

Explicit EGFR activation models (EAM)
We consider that the total number of both EGFR and EGF molecules are conserved, because we are studying early events of EGFR activation, 2 min after EGF addition, when protein degradation and synthesis are insignificant. Conservation laws are: Given these conservation relations, we will now introduce the ODEs of the model. So far, we enumerated 4(N+1) monomer species, 2 2 (N+1) 2 dimer species and two conservation relations. The rate of change of the species between these states is determined by the kinetics of ligand binding, receptor dimerization and phosphorylation. Since these dynamics do not affect each other (e.g., phosphorylation state of EGFR does not affect its ability to bind EGF), we wrote independent ODEs for each of them. The final ODEs, for each variable used in the simulations, take into account all the dynamics simultaneously.
Ligand binding to EGFR. Classically, EGF binding to EGFR is represented by the chemical

which for the open EGFR species reads as
where ! and ! are the binding and unbinding rate constants of EGF-EGFR binding, respectively.
Since we always use mass-action kinetics, the ODEs corresponding to reactions (18) are: At steady-state, the ratio of bound to unbound EGFR is determined by the dissociation constant  (18) and (19) ! , ! , 1 ! and 1 ! for ! , ! , ! and ! , respectively.
In EGFR dimers, we assumed that each monomer binds to EGF independently. Thus, the EGF binding rate constant to ! ∶ ! is twice that of the monomers. The same is true for the unbinding rate constant of ! ∶ ! . Using, Eq. (19) and the independence assumption, we obtained: In Eqs. (20), the symbol !" is Kronecker's delta, which takes a value of 1 when = , and 0 otherwise. Rather than explicitly solving an ODE for L, we computed the free ligand from the conservation relation Eq. (17). The same was done for species Rc 0 , for which we took advantage of the conservation relation for ! . We verified that Eqs. (19) and (20) the ratio of the two rates ! = ! ! determines, at steady-state, the ratio of closed to open EGFRs.

Definition of new dimeric species.
Despite the change of variables described so far, for N=9 we still have 400 dimer species, whose equations can be written simply by adding Eqs. (20) and (23). A second critical simplification comes by focusing on single receptors only, regardless of whether they are in a dimeric or monomeric configuration. Rather than following explicit dimers, we considered single molecules that are bound to a cognate EGFR. This class is denoted by a 'D', followed by an 'a' or an ' ' depending on whether the partner EGFR has an EGF bound or not, respectively ('a' stands for 'active', because phosphorylatable, ' ' for inactive because unphosphorylatable). If the receptor has an EGF bound, we also added a letter 'L', for ligand. With this simplification, for EGFR with N=9 Tyr, instead of following 400 dimeric species, we only needed to consider 40 EGFR dimer states, namely ! , ! , ! , and ! , where α = 0,..., 9 .
We have a factor of 2 increase from the 20 monomeric species because we distinguish whether their EGFR partner has an EGF moiety or not.
The remapping of the species follows the rules: In the following sections, we rewrite the rates for ligand binding and dimerization using these new variables.

EGF binding to collective species.
To obtain the ODEs of variables defined by Eqs. (26), we differentiated Eqs. (26) and substituted Eqs. (20). For example, for EGF binding to species Di α , we The ODEs for the three remaining D species were derived using the same procedure as outlined for !, giving the system of ODEs: Note that Eqs. (29) show no remaining dependence on the original dimeric species.
Dimerization for collective species. Differentiating Eqs. (26) and using Eqs. (23) we can write: In contrast to Eqs. (29), which only contain dimeric species, Eqs. (30) include also monomeric species. For example, the first equation tells us that ! becomes ! by binding to any monomer that has an EGF bound. Note that, each dimerization reaction appears twice on the right hand side of Eqs. (30), because it changes the concentration of two EGFRs in conformation 'D'. Nicely, these double contributions are automatically accounted for by the summations on the right hand side of the ODEs, as these ODEs satisfy the conservation law [Eq. (27)].

EGFR Phosphorylation for EAM species.
In ``EGFR phosphorylation as a multisite chain of events'', we explained how we reduced the kinase and phosphatase activity to a linear chain of increasingly phosphorylated EGFR species, i.e., to a multi-site phosphorylation model. In ``Freefree and saturated-saturated regimes studied with the Gillespie method'', we showed that since EGFR follows the free-enzyme regime for both phosphorylation and dephosphorylation, the phosphorylation of each Tyr is independent of that of the other Tyr, which allows us to model a 9Y

EGFR with only N=3 tyrosines. To obtain the phosphorylation/dephosphorylation ODEs for species
Da with N phosphorylatable Tyr, we employ the MPM in the free-free regime, which here translates into the following ODEs: To keep track of the three key Tyr, we solved numerically a model with N=3, which can be written explicitly as: , the other 'active' species, follows the same equations and can be obtained simply by replacing with . The four monomeric species , , , and , and the dimeric species and are inactive and thus have a negligible !"# , such that the equations reduce to dephosphorylation only: We obtained the ODEs for the dephosphorylation kinetics of , , , , and by replacing with these species in Eqs. (33) .

EAM phosphorylation in the absence of Cbl and Grb2
. This section extends the corresponding definitions given for the MPM-B model to the new model with explicit EGFR activation that contains a large number of species. We assumed that the anti-pY antibody recognizes pYs independently of EGFR conformation; thus, ! will have pYs as well as ! , ! , etc. Note that each of these terms depends upon EGF concentration, time, and all the parameters of the model (which we collectively represent as k ). It is convenient to define EGFR α (L, k,t) = R α + RL α + Rc α + RcL α + Di α + Da α + DiL α + DaL α ( ) (34) where we explicitly kept the ligand, time, and parameter dependence only on the left hand side of the equation. The overall phosphorylation level of the EGFR is finally given by where summation over ranges over all phosphorylation levels. Given the EGF concentration (L), to compute the total phosphorylation at any time t, we ran a simulation of EAM until time t and then used Eq. (35) to compute the total phosphorylation of the EGFR. To compute the phosphorylation of each of the species we consider that Cbl and Grb2 protect tyrosines from. The equations for the phosphorylations are analogous to Eq. (10), with the only exception that only the active species can be phosphorylated. To obtain the equations for the active species, we need to substitute and for R in Eq. (10), while for the inactive forms, i.e., for R , Rc , RL , RcL , Di and DiL , we also have to set !"# = 0 or, equivalently, remove the phosphorylation reactions altogether.
To compute the phosphorylation, for brevity of notation, we define a collective variable ! as the set of all possible EGFR species: ! = , , , , , , , . In other words, ! is the set of all 8 EGFR species and ∈ ! means that is any of the EGFR species. With this notation we have, for example, !∈! ! = + + + + + + + . We also defined 1PTyr, 2PTyr and 3PTyr as the sum of all EGFR species having 1, 2 and 3 phosphate groups, respectively: 1PTyr = Ρ P + Ρ P:Cbl + Ρ P:Cbl:Grb2 + Ρ P + Ρ P:Grb2 + Ρ P:Grb2:Cbl Finally, we defined total EGFR Tyr phosphorylation as Letting again Ρ T = {R, Rc, R, RL, Di, DiL, Da, DaL} , we define two species !"## and !" , which represent the fraction of EGFR with a moiety of Cbl either doubly or singly bound to pTyr1045, respectively. The two species are: where with ∈ ! can take any value of the set ! in the summation. The ubiquitination then equals, again, = !"## + !" and the normalized ubiquitination is computed as The same considerations made after Eq. (15) on ubiquination probabilities apply here.
We used Eqs. (37) and (39) to compute the curves shown in Fig. 6b-d, Fig. 7a, Fig. 7c, Fig. 8a-c and Supplementary Fig. 11. In the case of EAM, the PLE identifiability analysis classified all the EAM parameters as structurally identifiable. However, the amount of data available for the fit has been not enough to ensure also practical identifiability: all model parameters, with the exception of ! , ! !"# and !"# , have been classified as practically non-identifiable. 10 !! s !! µμM !! , respectively. Since the enzymatic parameters measured in these settings are unlikely to reproduce the kinetics in vivo of the full EGFR embedded on the plasma membrane, we decided to leave ample margin to the optimization algorithm.

Parameter optimization tools. Parameters were optimized for the EAM model (Supplementary
Mass spectrometry measurements of single EGFR tyrosines in vivo 1 show that tyrosine phosphorylation is maximal at around 2 minutes after EGF stimulus. Biochemical assays following kinetics of total EGFR phosphorylation or single tyrosine sites show maximal phosphorylation between 2 to 10 minutes depending on cell types 3 . Inhibition of the EGFR kinase at 2 or 10 minutes after EGFR stimulation, allowed for the estimation of the EGFR phosphorylation decay, which leads to the dephosphorylation rate constant ( !"! ) 3 . These experimental evidences suggest a phosphorylation time scale (the inverse of !"# ) of the order of tens of seconds. We use these qualitative results to set a (wide) range of optimization for the phosphorylation rate constant. We also fixed the dephosphorylation rate constant to the experimentally measured value for !"! at two min after EGF addiction reported in Kleiman et al. !"! = 0.016 !! 3 . Finally, mass spectrometry data 46 suggest that the simplifying assumption of using a single parameter for the phosphorylation kinetics of all EGFR tyrosines is acceptable in first approximation. There are reports of the dissociation constants for the binding of full length Grb2 to phosphopeptide pY1068 SH2 in solution 7,8,9 , and of Cbl TKB to phosphopeptide pY1045 10  Cbl and Grb2 interaction. Grb2 SH3 domains bind the Cbl proline-rich region. Grb2 has two SH3 domains, one located in the N-ter and one in the C-ter of the protein. Cbl has three proline-rich motifs located in the middle of the protein after the RING domain 47 . There are no studies that report in vitro measurements of Kd between purified Grb2 and Cbl. Instead, many studies focused on the binding between Grb2 and Sos-1 (which also contains proline-rich domains). Measurements in solution of SH3 domains of Grb2 with proline-rich region of Sos-1 vary a lot, in the range 1 nM-120 µM. This variability may depend, at least in part, on whether full-length Grb2 7,8,11 or isolated SH3 domains 9, 13 were used. In any case, all these studies suggested that binding of both N-ter and C-ter SH3 domains of Grb2 to multiple proline-rich sequences should decrease the Kd to the nM range [see also 48,49 ]. To be noted, Cbl has three proline-rich motifs that may contribute to the interaction with Grb2. In addition, early in vivo studies 50,51 , found that Cbl binds constitutively to Grb2 suggesting that Cbl-Grb2 binding may be stable. We thus set a wider optimization range for the dissociation constant , which represents the affinity of Cbl for Grb2.
We found no results for the dissociation rate of Cbl and Grb2, and therefore we leave an ample range of variation to the optimization algorithm for the corresponding parameter. by assuming that dimerization is diffusion limited and, consequently, using the theoretical expression ! !"# = 2 /ln ( / ) 52, 53 . We obtained ! !"# ≅ 2.7×10 !! µμm ! s !! setting the diffusion coefficient D for EGFR particles to ≈ 10 !!" cm ! s !! 19 , their radius to ≈ 5-7 nm (measured using the pdb ID 3I2T of the extracellular domain), and computing b, which is half the mean distance between receptors on a membrane, using the surface area A (we estimate a HeLa radius in the order of 15 and assume cells to be spherical) and the number of surface EGFR Γ T (of the order of 3-5×10 ! ) with the formula = / Γ ! . It is useful to transform the two dimensional rate constant ! !"# into a three dimensional one, by assuming that the third dimension is defined by the height of the EGFR, i.e. 10 nm. In the volume defined by the medium, the association rate constant for dimerization corresponding to the diffusion coefficient reported by Kusumi et al. is ! !"# ≅ 10nM !! s !! . Using all the estimates, we obtain the following range of optimization ! !"# ≅ 1--200nM !! s !! for the model parameter that corresponds to a range for the diffusion coefficient D, which is the measured observable, of 0.01--0.2 µμm ! s !!
There are various measurements of the backward dimerization rate constant ( ! !"# ), which lies in the range 0.27-2 s !! depending on the number of EGF bound and on the type of assay used 20,21 .

Supplementary note 3 -Sensitivity analysis
Although many parameters used were either available or determined ad hoc for this study, at least in terms of order of their magnitude, their precise values were fixed by means of a best-fit procedure described in Supplementary note 2. We thus decided to assess the robustness of the results to parameter variation. To this end, we performed a robustness analysis of the model, imposing a 1% increase and decrease for each parameter, to detect how the changes affect the output of the model in terms of steepness and half-maximal [EGF] of the dose-response curve for both phosphorylation and ubiquitination.
For each experiment analyzed in the main text [wild type (WT, Fig. 4a and Fig. 6b), Cbl overexpression (CblOE) and down-regulation (Cbl70Z) (Fig. 4b)], we measured the sensitivity coefficients € σ y i (k j ) = ∂ ln y i /ln k j = ∂y i / y i /∂k j /k j , where: 'i' identifies the observable € y i (either steepness n H or half-maximal value for total phosphorylation pY 0.5 or total ubiquitination x T ) and 'j' the parameter ! . In practice, we compute the sensitivity coefficient using the formula where stands for a vector containing all the parameters of the model, and * stands for the same vector with ! substituted by ! * . X (referring to either WT, CblOE or Cbl70Z), determines the experimental conditions corresponding to the computed sensitivity coefficient. Since we wanted to assess the sensitivity of both the WT and mutated cellular conditions, we measured € σ y i (k j ) = σ y i X (k j ) In Fig. 6c-d, we show the largest sensitivity coefficients computed using Eq. (41) for parameters variation of 1%.
The steepness of both EGFR phosphorylation and ubiquitination were remarkably robust: all sensitivity coefficients were below 0.5, which implies that this property of the curve is mainly encoded in the network wiring and not in the specific values of the parameters, at least around the best fit parameters. Even for the position of the half-maximal value, all sensitivity coefficients were below 1, with the most sensitive parameter being the first-order rate constant of phosphoryl group addition by kinases !"# (σ=0.9). Finally, it should be noted that all the Cbl-and Grb2-binding parameters showed a σ below 0.3, which ensures that primarily the phosphorylation kinetics determines the position of the ubiquitination threshold.