Abstract
Background
The Epidermal Growth Factor (EGF) receptor has been shown to internalize via clathrinindependent endocytosis (CIE) in a ligand concentration dependent manner. From a modeling point of view, this resembles an ultrasensitive response, which is the ability of signaling networks to suppress a response for low input values and to increase to a predefined level for inputs exceeding a certain threshold. Several mechanisms to generate this behaviour have been described theoretically, the underlying assumptions of which, however, have not been experimentally demonstrated for the EGF receptor internalization network.
Results
Here, we present a mathematical model of receptor sorting into alternative pathways that explains the EGFconcentration dependent response of CIE. The described mechanism involves a saturation effect of the dominant clathrindependent endocytosis pathway and implies distinct steadystates into which the system is forced for low vs high EGF stimulations. The model is minimal since no experimentally unjustified reactions or parameter assumptions are imposed. We demonstrate the robustness of the sorting effect for large parameter variations and give an analytic derivation for alternative steadystates that are reached. Further, we describe extensibility of the model to more than two pathways which might play a role in contexts other than receptor internalization.
Conclusion
Our main result is that a scenario where different endocytosis routes consume the same form of receptor corroborates the observation of a clearcut, stimulus dependent sorting. This is especially important since a receptor modification discriminating between the pathways has not been found experimentally. The model is not restricted to EGF receptor internalization and might account for ultrasensitivity in other cellular contexts.
Background
Endocytosis is the process by which activated transmembrane receptors are directed into the endosomal system from the plasma membrane [14]. In the past years, it has emerged as a powerful mechanism for the cell to temporally and spatially control its signaling response [5]. Ligand induced phosphorylation of EGF receptor creates docking sites for adaptor proteins, such as EPS15, epsin and AP2 [6,7]. Via direct or indirect binding, adaptors recruit the receptor to special membrane regions which are characterized by a particular composition of cageproteins and/or lipids [8,9]. The forming vesicles pinch off the membrane and carry their cargo to distinct intracellular locations, which might account for the specificity of the invoked signal [1,10]. Endocytosis may direct the receptors for lysosomal degradation or recycle them back to the membrane [1012]. Proper sorting of the EGF receptor into the correct endocytosis route is crucial for cell functioning as indicated by the fact that corruption of the sorting e.g. by viral proteins [13,14] may result in impaired receptor downregulation and increased mitogenic activity [15].
Clathrindependent endocytosis (CDE) was the first receptor internalization mechanism to be discovered and is generally considered the major route for EGF receptor (reviewed in [1,5,6,9]). Nevertheless, receptor internalization mechanisms that do not employ the structural protein clathrin, but arise from lipid rafts and caveolinrich membrane regions exist (Clathrinindependent endocytosis, CIE) [8,9,16,17]. The important question which molecular events govern the sorting of the EGF receptor into the different endocytosis pathways remains unanswered [5,8,9,1820].
A study addressing the sorting between Clathrin vs lipid raft/Caveolaemediated Endocytosis in mammalian cells suggested an interesting mechanism for the sorting process [21]: the distribution of receptors into the two pathways was shown to be EGFconcentration dependent. In the presence of low concentrations of EGF, the receptor was exclusively internalized via CDE, whereas at high concentrations, receptors were equally distributed between CDE and CIE (Figure 1).
Figure 1. CDE and CIE pathways of EGF receptor. An illustration of CDE and CIE pathways of EGF receptor. High EGF concentrations induce CIE, whereas CDE is observed at low and high EGF concentrations. The adaptors for the respective endocytosis pathways are referred to as CDE or CIEadaptors, respectively. See list of abbreviations.
From a modeling point of view, the behaviour of the clathrinindependent pathway resembles an ultrasensitive response: activation of the pathway is suppressed for low input EGF values, to reach the same level as the clathrindependent pathway for high input levels. Theoretically, several different mechanisms can explain ultrasensitive behaviour. Multisite modifications lead to a sigmoidal response of the modified molecule [2224], an effect that can be enhanced by consecutive arrangement in the form of cascades [2529] which has also been validated experimentally [30].
Other models of ultrasensitivity have been derived for MichaelisMenten type enzyme reactions: the presence of a stoichiometric inhibitor of an enzyme can suppress a reaction up to a certain threshold [28]. In (de)modification cycles ultrasensitivity occurs when the opposing enzymes work in the zeroorder regime [31], a mechanism which has been shown to work during morphogen directed pattern formation [32], or if the abundance levels of unmodified substrate and enzyme are sufficiently high [33]. Mathematical modeling has previously played a significant role in elucidating the mechanisms of EGF receptor signaling and endocytosis [3442]. In a series of quantitative studies the interaction between receptors and endocytosis machinery was evaluated [34,35,38,43]. Here, the existence of at least two distinct internalization pathways with different affinities for the EGF receptor was discovered [35,43]. In [21] it was reported that monoubiquitination (monoUb) of the EGF Receptor could only be observed at high EGF concentrations, raising the question whether monoUb might serve as a discriminative feature, which, when appended to the receptor, selectively targets the receptor to CIE [19,44]. This, however, conflicts with reports on the involvement of ubiquitinbinding adaptor proteins such as epsin and EPS15 during CDE [19,20,4549].
To address this controversy, we built a mathematical model of the sorting process. We address the functional consequences of different affinities with which internalization pathways are entered and explain how a switchlike response of CIE may result simply from a saturation effect of the CDE pathway. Together with the observation of EGFconcentration dependence of CIE, this analysis invites attention to an ultrasensitive regulatory mechanism for endocytic sorting. We give an analytical derivation of the switcheffect and derive regimes of reaction parameters and initial values for which the switch is preserved. Further, we describe its extensibility to more than two pathways. Importantly, the mechanism imposes only weak assumptions on the underlying interaction structure and parameter values. In summary, we give evidence for the hypothesis that the main purpose of postligand binding modifications of the EGF receptor such as ubiquitination does not lie in the discrimination between alternative endocytosis pathways.
Results
We built a system of ordinary differential equations (ODEs) that models the sorting of EGF receptor into clathrindependent or independent endocytosis pathways. The equations read:
The model contains the binding reaction of EGF to receptor R, which leads to one form of activated receptor (R_EGF), capable of entering clathrindependent or independent endocytosis.
In order to simulate the entry of ligandbound receptor into an endocytosis pathway, we introduced variables CDE and CIE, representing adaptors for clathrindependent and independent endocytosis which R_EGF can enter with rates k_{cde }and k_{cie}, respectively. The variables CDE and CIE represent the amount of the limiting factor in each pathway, which could be adaptor or cageproteins. We assume that the affinity of activated receptors R_EGF is significantly higher for the CDEpathway (k_{cde }≫ k_{cie}). To quantify the fraction of receptor going either pathway, we introduced variables R_{i_cde }and R_{i_cie}. The steadystate values of these variables represent the amount of R_EGF internalized via CDE and CIE, respectively. The model equations were derived according to the law of massaction [50].
CIEinternalization depends on number of receptors and strength of EGFstimulation
We systematically scanned the space of initial values of the model (equations 1.1 – 1.7) to investigate the effect of EGF stimulation on receptor distribution into CDE or CIE (see Methods). Figure 2 shows the time trajectories of R_{i_cde }and R_{i_cie }for three different classes of initial conditions, each represented by four different sets of values (see Figure caption). The three classes of initial values have distinct effects on the internalizationbehavior of R_EGF, the activated receptor. They represent assumptions on the relative quantities of EGFmolecules, unbound receptors and endocytosis adaptors.
Figure 2. Temporal Evolution of internalized receptors. Time trajectories of R_{i_cde }(blue) and R_{i_cie }(green) for three different classes of initial conditions. In case A (B) at least one of EGF_{0 }or R_{0 }stays below CDE_{0} (CDE_{0 }+ CIE_{0}). In case C, both EGF_{0 }and R_{0 }exceed CDE_{0 }+ CIE_{0}. Initial values were chosen arbitrarily such that these conditions are satisfied. In all three cases CDE_{0 }= 2.0, CIE_{0 }= 2.0 and (A) (EGF_{0}, R_{0}) = (1.2, 1.5), (1.6, 1.3), (2.3, 1.3), (1.4, 2.5); (B) (2.5, 3.0), (3.1, 2.4), (3.1, 4.5), (5.0, 2.8); (C) (4.3, 4.7), (5.0, 5.2), (5.5, 5.0), (5.3, 6.0).
The first class of initial values (Figure 2A) represents the case that either the initial number of unbound receptors (R_{0}) or EGFmolecules (EGF_{0}) is lower than the capacity of the CDEpathway (CDE_{0}) (Generally, X_{0 }denotes the initial value of variable X). This corresponds to an experimental setting where cells are stimulated with low EGFconcentrations, i.e. EGF_{0 }< CDE_{0}. Initial values from the second class are such that either R_{0 }or EGF_{0 }are below the capacity of both internalization pathways (CDE_{0 }+ CIE_{0}) (Figure 2B), whereas the third class reflects the case that both R_{0 }and EGF_{0 }exceed the capacity of both internalization pathways (Figure 2C).
It can be seen that in case A, R_{i_cie}production stays close to zero. In case B, internalization via CIE does occur, albeit to a lesser degree than CDE. For case C, receptors are equally partitioned between CDE and CIE.
Conditions on receptor number for switcheffect of CIEinternalization
The simulations shown in Figure 2 suggest conditions on the receptor number under which an EGFdependent switch of CIEinternalization will occur. If a cell possesses less receptors than CDEadaptors, CIEinternalization will be low independent of EGFstimulation (cf. Figure 2A). If the cell exhibits more receptors than adaptors for CDE, but less than for both pathways, then, for EGFstimulations exceeding CDE_{0}, a moderate fraction of receptor will internalize via CIE (cf. Figure 2B). Finally, if the amount of receptors is higher than the combined capacity of both pathways, CIEinternalization will be switched on equally strong as CDEinternalization for EGFstimulations that are higher than this combined capacity (cf. Figure 2C).
To test this hypothesis, we performed the following virtual experiment. We chose three sets of initial values for receptor R, CDE and CIEadaptors such that they fall within the three respective classes: R_{0 }< CDE_{0}, R_{0 }< CDE_{0 }+ CIE_{0 }or CDE_{0 }+ CIE_{0 }< R_{0}, and stimulated the system with increasing amounts of EGF. Figure 3 shows the steadystate amounts of R_{i_cde }(blue) and R_{i_cie }(green) as a function of EGF_{0}.
Figure 3. SteadyState behaviour of internalized receptors. Plotted are steadystate values of R_{i_cde }(blue) and R_{i_cie }(green) as a function of EGFstimulation (EGF_{0}). An ultrasensitive response of CIEinternalization with respect to EGF occurs if R_{0 }> CDE_{0 }(B) or R_{0} > CDE_{0}+CIE_{0 }(C), but not if R_{0 }< CDE_{0 }(A). In (B) and (C), when EGF_{0} exceeds the amount of CDE_{0}, CIEinternalization switches on abruptly, with a maximal response if (CDE_{0}+CIE_{0}) < min{CIE_{0}, R_{0}} (C). In all three cases CDE_{0 }= 1.0, CIE_{0 }= 1.0 and (A) R_{0 }= 0.8, (B) R_{0 }= 1.3 and (C) R_{0 }= 2.4.
As predicted, for receptor levels lower than CDE_{0 }(Figure 3A), CIEinternalization stays close to zero, independently of EGFstimulation. If the initial amount of receptors is greater than the capacity of CDE, CIEinternalization sets in abruptly, albeit to a moderate degree compared to CDEinternalization, for EGFstimulations greater than CDE_{0 }(Figure 3B). Finally, if the initial number of receptors is greater than the capacity of both pathways, the CIE pathway switches on to an equal extent as CDEinternalization (Figure 3C).
We have thus derived an ultrasensitive response of CIEinternalization with respect to EGFstimulation, without assuming any discriminative receptor modifications. Rather, it is necessary and sufficient that the initial amount of receptors is higher than the capacity of the CDEpathway (CDE_{0 }< R_{0}, moderate switch) or both pathways (CDE_{0 }+ CIE_{0 }< R_{0}, maximal switch).
Correspondence to distinct classes of steadystate
For dynamical systems with multiple steadystates, a certain steadystate will be reached depending on whether the system starts in the corresponding basin of attraction [23,51]. Thus, a switch between steadystates occurs for different vectors of initial values, provided that the separatrix, i.e. the hypersurface between neighboring basins of attraction, is crossed.
We investigated, whether the switcheffect of CIEinternalization corresponds to such a transition between distinct steadystates of the system. Analytically, one derives two classes of steadystates (see Methods for derivation of conditions):
where X* denotes the steadystate concentration of the respective component.
Note that classes of steadystates are used since not all variables are assigned specific values. For example, in both cases and are not uniquely determined and depend on the corresponding initial values. Steadystate class I reflects the case that either all available EGF (EGF* = 0) or all free receptors R (R* = 0) have been absorbed in the binding reaction and all activated receptors R_EGF have been internalized. In steadystate class II neither receptors nor ligand are limiting for the internalization process and have come to an equilibrium with R_EGF. Instead, the capacity of both internalization pathways has been depleted (CIE* = CDE* = 0).
The systematic scan of initial values and subsequent solving of the system until steadystate revealed initial conditions under which each steadystate class is reached. We found that if both EGFstimulation and initial receptor level are higher than the capacity of both internalization pathways (CDE_{0 }+ CIE_{0 }< min{R_{0}, EGF_{0}}, initial value class C), the system tends towards steadystate class II. Otherwise, steadystate class I will be reached. In this case, if EGFstimulation is below the amount of receptors, all EGF will be depleted in the binding reaction (EGF* = 0), whereas if it is above, all receptors will be consumed (R* = 0).
This is exemplified in Fig. 4, where the steadystate value of ligandbound receptor (R_EGF*) is plotted as a function of EGFstimulation for different initial receptor levels R_{0}. Here, CIE_{0 }= CDE_{0 }= 1. For R_{0 }= 1.7 (orange), i.e. R_{0 }< CDE_{0 }+ CIE_{0}, the system reaches steadystate class I independently of EGFstimulation, as seen from R_EGF* = 0. If R_{0 }= 3 (green) or R_{0 }= 5(black), i.e. CDE_{0 }+ CIE_{0 }< R_{0}, the steadystate value of R_EGF becomes positive for EGFstimulations higher than 2 (steadystate class II).
Figure 4. Dependence on initial values for steadystate. Plotted are steadystate values of ligandbound receptor (R_EGF*) as a function of EGFstimulation for different initial receptor levels R_{0}. For R_{0 }= 1.7 (orange), i.e. R_{0 }< CDE_{0 }+ CIE_{0}, the system reaches steadystate class I independently of EGFstimulation, as seen from R_EGF* = 0. If R_{0 }= 3 (green) or R_{0 }= 5 (black), i.e. CDE_{0 }+ CIE_{0 }< R_{0}, the steadystate value of R_EGF becomes positive for EGFstimulations higher than 2 (steadystate class II).
Applying these derived conditions on the initial values, we can also show that the steadystates are stable. A steadystate is stable, if, for small perturbations, the system returns to this steadystate. Consider steadystate class I with R* = 0 and R_EGF* = 0. CIE* and CDE* are not clearly defined in this case, but according to the conditions on the initial values we derived, we know that R_{0 }< CDE_{0 }+ CIE_{0}. This means that in steadystate, at least one of the two adaptor variables must be greater than zero, i.e. 0 = R* < CDE* + CIE*. If we apply a sufficiently small perturbation to the system, and set the obtained value as the new start vector, this last inequality will still hold due to the continuity of the functions. According to the conditions on initial values we derived in the previous paragraph, the system will tend back to R* = 0 and R_EGF* = 0. Hence we showed stability of steadystate class I. An analogous argument can be used to show stability of steadystate class II.
In [21] it was reported that for high ligand concentrations the activated receptor is equally partitioned between CDE and CIE. Assuming similar initial abundance levels of adaptors, our simulations show that this is the case for initial receptor levels higher than the sum of both initial adaptor values (Fig. 2C, Fig. 3C). We thus hypothesize that in cells, where the steadystate levels of internalized receptors via CIE and CDE are similar, the amount of receptors exceeds the capacity of both pathways. In this case, treatment of the cells with low vs high EGFstimulations, corresponds to a transition of the system between steadystate classes I and II (see Table 1).
Table 1. Conditions on initial values for steadystates. Assuming that the initial number of receptors is higher than the combined capacity of both pathways (R_{0 }> CDE_{0 }+ CIE_{0}), low and high EGFstimulations lead to two different steadystates, respectively. In steadystate I, the receptor internalizes primarily via CDE, whereas in steadystate II it is equally partitioned between CDE and CIE.
Steepness of switch effect
An ultrasensitive response of a signaling system is characterized by a low, or damped response up to a certain threshold of stimulus, followed by an abrupt increase towards maximal response when this threshold is exceeded [23,26,27,50]. It has been derived to result from positive feedback or multisitemodification [24,25,50,52,53]. Ultrasensitivity has also been shown to arise in (de)modification cycles if the enzymes operate near saturation [31], which makes the mechanism very sensitive to small parameter changes [26], if the abundance levels of unmodified substrate and enzyme are sufficiently high (ultrasensitization, [33]) or if the enzyme is inhibited [28].
To characterize the steepness of the here discussed mechanism, we compared its response to a Hilltype reaction (see Methods). Figure 5 shows the reaction velocity V of the Hillformula, compared to production in our model (stimulusresponse curve) as a function of EGFstimulation. To generate the stimulusresponse curve, we chose the same parameter set as for Fig. 3C as a reference. From this curve we extracted the Hillcoefficienth, V_{max }and K_{m }to compute the corresponding Hillcurve, which will be used as a reference curve later on. The Hillcoeffcient is a measure of how much the input has to be increased in order to raise the response from 10% to 90% of its maximal value [28]. Stimulusresponse curves with Hillcoefficients of 5 or higher are generally considered ultrasensitive [25,26,28]. The Hillcoefficient obtained for the stimulusresponse curve shown in Figure 5 is 7.5.
Figure 5. Approximability of switcheffect by Hillcurve. Comparison between Hilltype response and the described switch effect. Plotted are steadystate values of R_{i_cie }or reaction velocity V (Hillkinetics). Parameter and initial values used were the same as for Fig. 3C. Hillparameters extracted from the stimulusresponse curve were: h = 7.5, V_{max }= 1.0, K_{m }= 1.4.
Robustness of solution
It has been argued that biologically functional modules or pathways need to be robust against variations of reaction parameters and protein concentrations in order to ensure proper functioning [54,55]. The concept of robustness refers to the 'purpose' of a certain module or pathway: it is expected that intracellular network structures have undergone an evolution that guarantees their proper functioning independently of precise parameter values [56].
To transfer this concept to the question of receptor sorting into alternative pathways, we asked to what extent the functioning of the here described module depends on exact parameter or initial values. As functioning we defined the clearcut sorting of the receptor into distinct routes, namely CDE at low, respectively CDE and CIE at high ligand concentrations.
The key parameter and initial concentrations that affect the strength of the switch effect are the initial concentrations of CDE and CIEadaptors as well as . Obviously, the ultrasensitive response will be steeper the higher the difference in binding kinetics for the respective pathways is, i.e. the greater r_{k}. We systematically varied r_{k }and from each thus obtained stimulus response curve of values extracted the Hillcoefficient h as a measure of steepness (see Methods). The Hillcoefficients varied from 2.8 (r_{k }= 3.3) to 7.7 (r_{k }= 200) as shown in Figure 6. We also computed the mean deviation between the obtained stimulus response curves and the reference Hillcurve (Figure 7), showing that for considerable variations of r_{k }the stimulusresponse curve remains approximable by a Hillcurve.
Figure 6. Robustness of switcheffect for parameter variations. Comparison between stimulusresponse curve (with Hillcoefficient h) and corresponding Hillcurve for selected r_{k }values. Here, k_{cde }= 1.0 and k_{cie }was varied. Plotted are steadystate values of R_{i_cie }(blue) or reaction velocity V (Hillkinetics, magenta).
Figure 7. Deviation from Hillcurve. Mean deviation from reference Hillcurve and Hillcoefficient of stimulus response curves for varying . Here, k_{cde }= 1.0 and k_{cie }was varied between 0.01 (r_{k }= 100) and 0.6 (r_{k }= 1.67). All other parameter and initial values were as in Fig. 3C.
In Figure 8A we plotted the steadystate values of R_{i_cde }and R_{i_cie }(V_{max }values of the stimulusresponse curves) as a function of initial adaptor values CDE_{0 }and CIE_{0}. Here, R_{0 }and EGF_{0 }were chosen 1.5.
Figure 8. Dependence of switch effect on abundance levels of endocytosis adaptors. Plotted are (blue) and (orange) as a function of CDE_{0 }and CIE_{0}. EGF_{0 }= R_{0 }= 1.5. See text for interpretation.
Consider the curve for (blue). For initial adaptor values such that CDE_{0 }+ CIE_{0 }< R_{0}, EGF_{0 }(see arrow), the curve is largely independent of CDE_{0 }and increases linearly as a function of CIE_{0 }up to the threshold of 1.5. The independence of CDE_{0 }reflects the fact that if neither receptors nor ligand are limiting for the internalization reaction (steadystate class II), the steadystate amount of receptor internalized via CIE is solely dependend on the initially available number of CIEadaptors. Outside of this range, decreases with increasing CDE_{0 }and becomes zero for CDE_{0 }> 1.5(cf. Figure 3A). (orange curve) is largely independent of CIE_{0 }and increases linearly with CDE_{0 }up to the threshold of 1.5 when ligand or receptor number become limiting.
The threshold of CIEinternalization (K_{m }of the stimulusresponse curves) is independent of CIE_{0 }(for CIE_{0 }> 0) and is equal to CDE_{0 }as shown in Figure 8B.
Role of Receptor modifications
It is wellknown that ligandinduced receptor modifications in the form of phosphorylation and/or ubiquitination play a functional role in signaling and contribute to the specificity of adaptorbinding. However, our analysis focused on the question, whether for a precise sorting of receptors into the two alternative endocytosis pathways discussed here a discriminative modification is necessary. In this light, R_EGF, which in our model indicates the activated receptor species capable of interacting with the endocytosis adaptors, could also represent an already modified form of the receptor. To illustrate this point, we extended the model as given in equations (1.1 – 1.7) to include the binding of the ubiquitinligase Cbl followed by ubiquitination of the receptor. The equations read
Here, k_{onCbl }and k_{offCbl }are the rate constants for the association and dissociation of Cbl to ligandbound receptor R_EGF, and k_{catCbl }is the rate of the ubiquitination step.
Again, we tested whether the assumption that both pathways consume the thus, i.e. equally, modified form of receptor would comply with the observation of an ultrasensitive response of CIEinternalization to increasing EGFstimulation. In Figure 9 we plotted the steadystate values of receptor internalized via CDE(R_{i_cde}, blue) and CIE (R_{i_cie}, green), respectively. Clearly, the response is comparable to the results of the simpler model (Figure 3C), proving that the existence of receptor modifications prior to internalization does not affect our results. The rate constants and initial values : k_{f }= 1.0, k_{r }= 0.01, k_{onCbl }= 1.0, k_{offCbl }= 0.01, k_{catCbl }= 1.0, k_{cde }= 1.0, k_{cie }= 0.01, R_{0 }= 2.0, Cbl_{0 }= 2.0, CDE_{0 }= 1.0, CIE_{0 }= 1.6. All other initial values are zero.
Figure 9. Inclusion of Receptor Modifications. Plotted are steadystate values of R_{i_cde }(blue) and R_{i_cie }(green) as a function of EGFstimulation (EGF_{0}) for the extended model, including receptor ubiquitination. The switch effect is preserved under the assumption of receptor modifications.
Model extension for more than two Pathways
For EGF receptor, evidence for the existence of more than just two independent endocytosis pathways has been given [19,57]. To our knowledge, EGFconcentration dependence has not been shown. Here, we describe an extension of the above model to more than two pathways which might play a role in other contexts.
Suppose that n > 2 pathways branch off from one activated signaling molecule R_{L }(playing the role of ligandbound receptor R_EGF). Assume that R_{L }can bind to n different types of adaptor molecules C_{i }with reaction rates k_{i}, i = 1,2, ..., n where k_{i }>> k_{i + 1}, i = 1,2, ... n  1. Then the sorting effect based on ultrasensitivity is extended to m cases if , m ≥ n, with X_{0 }denoting the initial value of molecule class X. The effect is illustrated in Fig. 10B for n = 4. The differential equations read:
Figure 10. Model extension to n > 2 pathways. (A) Extension of the model for n > 2 cases. (B) Simulation for n = 4 cases. Plotted are steadystate amounts of R_{i}, i = 1, 2, 3, 4. Rate constants and initial values used: k_{f }= 1.0, k_{r }= 0.01, k_{1 }= 2.0, k_{2 }= 0.05, k_{3 }= 0.005, k_{4 }= 0.001, R_{0 }= 9.0, C_{1 }= 1, C_{2 }= 1, C_{3 }= 1, C_{4 }= 1. All other initial values were zero.
Discussion
Our analysis addresses the experimentally observed dependence of endocytosis on EGFconcentration [5,21]. We propose an ultrasensitive sorting mechanism for EGF receptor internalization which does not require a discriminative receptor modification and give a systematic description of the parameter requirements to achieve proper sorting. We derive analytically the existence of alternative steadystates as well as conditions on the abundance levels of receptors, ligand and endocytic adaptors to reach these states.
Referring to previous models of sigmoidal responses based on cooperativity, for the EGF receptor in particular a cooperative binding effect of the ubiquitinligase Cbl during the ubiquitination reaction has been proposed to be necessary for the observed switcheffect of CIE internalization [44]. Our analysis explains how imposing weaker assumptions on the internalization machinery, namely that the two pathways are entered with distinct affinities, is sufficient to explain the observed switcheffect. Importantly, it is pointed out how by varying the abundance levels of active receptors or endocytic adaptors, cells may modulate their response to incoming EGFstimulations: depending on the initially available receptors or adaptors, the distribution of internalized receptors can be different for one and the same EGF concentration (cf. Figure 3, Figure 8A).
The lack of knowledge about the true parameter/initial values was accounted for by systematic variations over broad numerical ranges. The robustness of the switcheffect to exact parameter values argues for the plausibility of the introduced mechanism.
Mathematical models addressing the problem of receptor sorting into alternative endocytosis pathways do not currently exist. Previously proposed hypotheses based on experimental data only have not been able to give a satisfying answer to this question [9,12,1821,58]. Generally, the problem is considered at the 'singlemoleculelevel' : a single receptor is envisioned, which is thought to enter either the CDE or the CIEroute (see Figure 1). This picture misleadingly implies the necessity of a discriminative receptor modification.
Instead of thinking in terms of individual entities, we propose to consider the dynamical properties of a system of interacting molecule populations. Applying methods from the theory of dynamical systems, we were able to conceive that an increase in the ligand concentration above the capacity of the CDEpathway qualitatively alters the system behaviour by enforcing an alternative steadystate (Fig. 3, Table 1). Our model states that an abrupt, switchlike start of CIE occurs if the extracellular EGF concentration exceeds the capacity of the CDE machinery. This proposes an interesting implication of the regulation of receptor sorting: the cell achieves the switcheffect 'for free' since no extra cost has to be invested into a discriminative receptor modification. It can be assumed that cells have evolved to optimize energy efficiency [59]. Utilizing the kind of dynamics introduced here, where just one form of receptor is consumed by both pathways, could thus constitute an evolutionary advantage.
A second major observation we draw from the model is that the described mechanism provides a means for an individual cell to sense its surrounding medium: clathrinindependent endocytosis is switched on precisely when the extracellular ligand concentration exceeds the number of CDEadaptors. One might interpret this mechanism as a protein module, i.e. a small interaction network acting as a computational element, whose purpose is to store and process information [51,60,61].
One might argue that the simplicity of the model impairs its abilitiy to uncover unanticipated results. While an initiation of clathrinindependent endocytosis upon saturation of the clathrin pathway might have been proposed without mathematical modeling, the unexpected steepness of this switchingbehaviour as well as its robustness could not have been revealed by intuition alone [9,18,19]. Furthermore, we showed that extending the model by allowing a modification of the receptor does not increase the steepness of the response. Thus, we conclude that a modification of receptor is not required to discriminate between the pathways. This does notably not exclude the possibility that a modification of the receptor might have been chosen by nature to ensure proper endocytic sorting.
Finally, we want to emphasize that the generality of our model makes it applicable to ultrasensitivity in signaling processes other than the here discussed problem of receptor sorting. This is highlighted in the potential of the model to be extended to more than two overlapping signaling pathways (cf. Figure 9)
Conclusion
We describe the dynamical consequences of an interaction motif, whose molecular basis has already been well established experimentally and discuss its applicability to endocytic sorting of the EGF receptor. We give a simple, yet mathematically sound explanation how cell variation of endocytic sorting results from modulating abundance levels of involved cellular molecules. In the light of the difficulty to experimentally identify a discriminative receptor modification, our analysis points to the possibility of a systemslevel regulation of endocytic sorting. The natural extensibility of the model to more than two cases may prove applicable in other signaling contexts.
Methods
Numerical Simulations and SteadyState Analysis
To numerically solve the systems of equations (1.1 – 1.7) we used the MATLAB ODE15s function. The existence of distinct classes of steadystates was derived as follows: The system of ODEs exhibits four independent equations, namely equations (1.1), (1.3), (1.4) and (1.5). To obtain the values of the variables that represent a steadystate, we simultaneously set these equations equal to zero and solved for the variables.
Setting equations (1.4) and (1.5) equal to zero yields
Since all kinetic constants are assumed to be positive, one derives that either R_EGF* or both CDE* and CIE* have to be zero. Consider the case that R_EGF* = 0. Then, from setting equations (1.1) and (1.3) equal to zero, i.e.
it follows that either R* or EGF* have to be zero. We have thus derived steadystate class I (see Results). Consider now the second case derived from equations (1.4a) and (1.5a), namely that CDE* = 0 and CIE* = 0. Substituting these values into equation (1.3), one obtains
which also leaves equation (1.1) equal to zero. This second set of values represents steadystate class II (see Results).
Approximability by Hillcurve
In order to assess the steepness of the derived switcheffect, we compared the obtained stimulusresponse curve to the switcheffect described by the Hillformula. This formula was originally introduced as a phenomenological model to describe the binding of oxygen to hemoglobin [62] and later often used to for cooperative enzyme reactions:
Here, V denotes reaction velocity, V_{max }the maximal velocity, x the substrate concentration, K_{m }the substrate concentration where halfmaximal velocity is reached and h, the Hillcoefficient, represents the steepness of the response.
In order to investigate the sensitivity of the switcheffect of CIEinternalization on specific parameter values, we adapted a procedure applied in [26]. We systematically varied reaction parameters and initial values of the model (equations 1.1 – 1.7) using equidistant sampling points. For each thus obtained set we computed the stimulusresponse curve and extracted the Hillcoefficient as [25]. Here, x_{0.1 }and x_{0.9 }are the EGFconcentrations where 10% or 90% of the maximal responseis reached, espectively.
To evaluate the approximability of the stimulusresponse curve by the Hillequation, we computed the mean deviation between the obtained stimulusresponse curves and a reference Hillcurve. In this way we formally determined a broad range of parameter values for which an ultrasensitive response of CIEinternalization occurs (Fig. 5 and Fig. 6).
List of Abbreviations
The following abbreviations are used: CDE – Clathrin dependent endocytosis, CIE – Clathrin independent endocytosis, EGF – Epidermal Growth Factor, EGFR – Epidermal Growth Factor Receptor, monoUb – monoubiquitination
Authors' contributions
HSG designed the study, performed the simulations and mathematical analysis and wrote the manuscript. IV designed the study, helped with the mathematical analysis and wrote the manuscript. DH proposed the problem of endocytic sorting for mathematical modeling and helped design the study. ID proposed the problem of endocytic sorting for mathematical modeling and helped design the study. RE proposed the problem of endocytic sorting for mathematical modeling, designed the study and wrote the manuscript.
Acknowledgements
We would like to thank Frauke Melchior, Angela Stevens and Angel Alonso for helpful comments on the manuscript. We further acknowledge fruitful discussions on modelling aspects with Ivano Primi, Anna Marciniak, Sven Mesecke and Hauke Busch. This project was supported by the SMP Bioinformatics within the National Genome Research Network (NGFN) funded by the Federal Ministry of Education and Research (BMBF) and Endocyte (EU Foerderung).
References

Ceresa BP, Schmid SL: Regulation of signal transduction by endocytosis.
Curr Opin Cell Biol 2000, 12(2):204210. PubMed Abstract  Publisher Full Text

Miaczynska M, Pelkmans L, Zerial M: Not just a sink: endosomes in control of signal transduction.
Curr Opin Cell Biol 2004, 16:400406. PubMed Abstract  Publisher Full Text

Wiley HS: Trafficking of the ErbB receptors and its influence on signaling.
Exp Cell Res 2003, 284:7888. PubMed Abstract  Publisher Full Text

Perret E, Lakkaraju A, Deborde S, Schreiner R, RodriguezBoulan E: Evolving endosomes: how many varieties and why?
Curr Opin Cell Biol 2005, 17(4):423434. PubMed Abstract  Publisher Full Text

Polo S, Fiore PPD: Endocytosis conducts the cell signaling orchestra.
Cell 2006, 124(5):897900. PubMed Abstract  Publisher Full Text

Kirchhausen T: Clathrin Adaptors Really Adapt.
Cell 2002, 109(4):413416. PubMed Abstract  Publisher Full Text

Sorkin A: Cargo recognition during clathrinmediated endocytosis: a team effort.
Curr Opin Cell Biol 2004, 16(4):392399. PubMed Abstract  Publisher Full Text

Simons K, Toomre D: Lipid Rafts and Signal Transduction.
Nat Rev Mol Cell Biol 2000, 1:3139. PubMed Abstract  Publisher Full Text

Le Roy C, Wrana J: Clathrin and nonclathrin mediated endocytic regulation of cell signalling.
Nat Rev Mol Cell Biol 2005, 6:112126. PubMed Abstract  Publisher Full Text

Hoeller D, Volarevic S, Dikic I: Compartmentalization of growth factor receptor signalling.
Curr Opin Cell Biol 2005., 17(2) PubMed Abstract  Publisher Full Text

Dikic I: Mechanisms controlling EGF receptor endocytosis and degradation.
Biochem Soc Trans 2003, 31:11781181. PubMed Abstract  Publisher Full Text

Di Guglielmo GM, Le Roy C, Goodfellow AF, Wrana JL: Distinct endocytic pathways regulate TGF[beta] receptor signalling and turnover.
Nat Cell Biol 2003, 5(5):410421. PubMed Abstract  Publisher Full Text

Straight SW, Hinkle PM, Jewers RJ, McCance DJ: The E5 oncoprotein of human papillomavirus type 16 transforms fibroblasts and effects the downregulation of the epidermal growth factor receptor in keratinocytes.
J Virol 1993, 67(8):452132. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Straight SW, Herman B, McCance DJ: The E5 oncoprotein of human papillomavirus type 16 inhibits the acidification of endosomes in human keratinocytes.
J Virol 1995, 69(5):318592. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Waterman H, Yarden Y: Molecular mechanisms underlying endocytosis and sorting of ErbB receptor tyrosine kinases.
FEBS Lett 2001, 490(3):142152. PubMed Abstract  Publisher Full Text

Razani B, Woodman SE, Lisanti MP: Caveolae: From Cell Biology to Animal Physiology.
Pharmacol Rev 2002, 54(3):431467. PubMed Abstract  Publisher Full Text

Lucas PelkmansBTZM, Helenius A: CaveolinStabilized Membrane Domains as Multifunctional Transport and Sorting Devices in Endocytic Membrane Traffic.
Cell 2004, 118:767780. PubMed Abstract  Publisher Full Text

FelberbaumCorti M, Van Der Goot FG, Gruenberg J: Sliding doors: clathrincoated pits or caveolae?
Nat Cell Biol 2003, 5(5):382384. PubMed Abstract  Publisher Full Text

Aguilar RC, Wendland B: Endocytosis of membrane receptors: Two pathways are better than one.
PNAS 2005, 102(8):26792680. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hoeller D, Dikic I: Receptor endocytosis via ubiquitindependent and independent pathways.
Biochem Pharmacol 2004, 67:10131017. PubMed Abstract  Publisher Full Text

Sigismund S, Woelk T, Puri C, Maspero E, Tacchetti C, Transidico P, Di Fiore PP, Polo S: From the Cover: Clathrinindependent endocytosis of ubiquitinated cargos.
PNAS 2005, 102(8):27602765. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Segel LA: Modeling dynamic phenomena in molecular and cellular biology. Cambridge University Press; 1984.

EdelsteinKeshet L: Mathematical Models in Biology. SIAM, Philadelphia, USA; 2005.

Markevich NI, Hoek JB, Kholodenko BN: Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades.
J Cell Biol 2004, 164:353359. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Huang CYF, Ferrell J, James E: Ultrasensitivity in the mitogenactivated protein kinase cascade.
PNAS 1996, 93:1007810083. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bluthgen N, Herzel H: How robust are switches in intracellular signaling cascades?
J Theor Biol 2003, 225(3):293300. PubMed Abstract  Publisher Full Text

Ortega F, Acerenza L, Westerho3 HV, Mas F, Cascante M: Product dependence and bifunctionality compromise the ultrasensitivity of signal transduction cascades.
PNAS 2002, 99:11701175. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ferrell JE: Tripping the switch fantastic: how a protein kinase cascade can convert graded inputs into switchlike outputs.
Trends Biochem Sci 1996, 21(12):460466. PubMed Abstract  Publisher Full Text

Ferrell JE: How responses get more switchlike as you move down a protein kinase cascade.
Trends Biochem Sci 1997, 22(8):288289. PubMed Abstract  Publisher Full Text

Ferrell JE, Machleder EM: The biochemical basis of an allornone cell fate switch in Xenopus oocytes.
Science 1998, 280:8958. PubMed Abstract  Publisher Full Text

Goldbeter A, Koshland JDE: An amplified sensitivity arising from covalent modification in biological systems.
PNAS 1981, 78:68406844. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Melen GJ, Levy S, Barkai N, Shilo BZ: Threshold responses to morphogen gradients by zeroorder ultrasensitivity.
Mol Syst Biol 2005, 1:msb4100036E1msb4100036E11. Publisher Full Text

Legewie S, Blthgen N, Schfer R, Herzel H: Ultrasensitization: switchlike regulation of cellular signaling by transcriptional induction.
PLoS Comput Biol 2005, 1(5):e54. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wiley HS, Cunningham DD: A steady state model for analyzing the cellular binding, internalization and degradation of polypeptide ligands.
Cell 1981, 25(2):433440. PubMed Abstract  Publisher Full Text

Starbuck C: Mathematical model for the effects of epidermal growth factor receptor trafficking dynamics on fibroblast responses.
Biotechn Prog 1992, 8:132143. Publisher Full Text

Schoeberl B, EichlerJonsson C, Gilles ED, Mueller G: Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors.
Nat Biotechnol 2002, 20(4):3705. PubMed Abstract  Publisher Full Text

Lauffenburger D, Linderman J: Receptors. Models for Binding, Trafficking and Signaling. Douglas Lauffenburger; 1993.

Wiley HS, Cunningham D: The endocytotic rate constant. A cellular parameter for quantitating receptormediated endocytosis.
J Biol Chem 1982, 257(8):42224229. PubMed Abstract  Publisher Full Text

Kholodenko BN, Demin OV, Moehren G, Hoek JB: Quantification of short term signaling by the epidermal growth factor receptor.
J Biol Chem 1999, 274(42):3016981. PubMed Abstract  Publisher Full Text

Wiley HS, Shvartsman SY, Lauffenburger DA: Computational modeling of the EGFreceptor system: a paradigm for systems biology.
Trends Cell Biol 2003, 13:4350. PubMed Abstract  Publisher Full Text

Resat H, Ewald JA, Dixon DA, Wiley HS: An integrated model of epidermal growth factor receptor trafficking and signal transduction. [http://www.biophysj.org/cgi/content/full/85/2/730] webcite
Biophys J 2003, 85(2):73043. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hendriks BS, Wiley HS, Lauffenburger D: HER2mediated effects on EGFR endosomal sorting: analysis of biophysical mechanisms. [http://www.biophysj.org/cgi/content/full/85/4/2732] webcite
Biophys J 2003, 85(4):273245. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lund K, Opresko L, Starbuck C, Walsh B, Wiley H: Quantitative analysis of the endocytic system involved in hormone induced receptor internalization.
J Biol Chem 1990, 265(26):1571315723. PubMed Abstract  Publisher Full Text

Grant BD, Audhya A: The ins and outs of endocytic transport.
Nat Cell Biol 2005, 7(12):11511154. PubMed Abstract  Publisher Full Text

Coda L, Salcini AE, Confalonieri S, Pelicci G, Sorkina T, Sorkin A, Pelicci PG, Fiore PPD: Eps15R is a tyrosine kinase substrate with characteristics of a docking protein possibly involved in coated pitsmediated internalization.
J Biol Chem 1998, 273(5):30033012. PubMed Abstract  Publisher Full Text

Sorkina T, Bild A, Tebar F, Sorkin A: Clathrin, adaptors and eps15 in endosomes containing activated epidermal growth factor receptors.
J Cell Sci 1999, 112(Pt 3):317327. PubMed Abstract  Publisher Full Text

Tebar F, Sorkina T, Sorkin A, Ericsson M, Kirchhausen T: Eps15 is a component of clathrincoated pits and vesicles and is located at the rim of coated pits.
J Biol Chem 1996, 271(46):2872728730. PubMed Abstract  Publisher Full Text

Stang E, Blystad FD, Kazazic M, Bertelsen V, Brodahl T, Raiborg C, Stenmark H, Madshus IH: Cbldependent ubiquitination is required for progression of EGF receptors into clathrincoated pits. [http://dx.doi.org/10.1091/mbc.E04010041] webcite
Mol Biol Cell 2004, 15(8):35913604. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

de Melker AA, van der Horst G, Borst J: Ubiquitin Ligase Activity of cCbl Guides the Epidermal Growth Factor Receptor into Clathrincoated Pits by Two Distinct Modes of Eps15 Recruitment.
J Biol Chem 2004, 279(53):5546555473. PubMed Abstract  Publisher Full Text

Tyson JJ, Chen KC, Novak B: Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell.
Curr OpinCell Biol 2003, 15:221231. Publisher Full Text

Bhalla US, Iyengar R: Emergent properties of networks of biological signaling pathways.
Science 1999, 283:3817. PubMed Abstract  Publisher Full Text

Angeli D, Ferrell J, James E, Sontag ED: Detection of multistability, bifurcations, and hysteresis in a large class of biological positivefeedback systems.
PNAS 2004, 101:18221827. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Onn Brandman RL, Ferrell JamesE Jr, Meyer T: Interlinked Fast and Slow Positive Feedback Loops Drive Reliable Cell Decisions.
Science 2005, 310:496498. PubMed Abstract  Publisher Full Text

Barkai N, Leibler S: Robustness in simple biochemical networks.
Nature 1997, 387:913917. PubMed Abstract  Publisher Full Text

Hartwell LH: Theoretical biology: A robust view of biochemical pathways.
Nature 1997, 387:855. PubMed Abstract  Publisher Full Text

Kitano H: Biological Robustness.
Nature Reviews Genetics 2004, 5:826837. PubMed Abstract  Publisher Full Text

Hinrichsen L, Harborth J, Weber K, Ungewickell E: Effect of clathrin heavy chain and alphaadaptinspecific small inhibitory RNAs on endocytic accessory proteins and receptor trafficking in HeLa cells.
J Biol Chem 2003, 278(46):4516070. PubMed Abstract  Publisher Full Text

Grovdal LM, Stang E, Sorkin A, Madshus IH: Direct interaction of Cbl with pTyr 1045 of the EGF receptor (EGFR) is required to sort the EGFR to lysosomes for degradation.
Exp Cell Res 2004, 300(2):388395. PubMed Abstract  Publisher Full Text

Dekel E, Alon U: Optimality and evolutionary tuning of the expression level of a protein.
Nature 2005, 436:58892. PubMed Abstract  Publisher Full Text

Bray D: Protein molecules as computational elements in living cells.
Nature 1995, 376:30712. PubMed Abstract  Publisher Full Text

Hartwell LH, Hopfield JJ, Leibler S, Murray AW: From molecular to modular cell biology.
Nature 1999, 402:C47C52. PubMed Abstract  Publisher Full Text

Hill A: The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves.