An ultrasensitive sorting mechanism for EGF Receptor Endocytosis

Background The Epidermal Growth Factor (EGF) receptor has been shown to internalize via clathrin-independent 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 pre-defined 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 EGF-concentration dependent response of CIE. The described mechanism involves a saturation effect of the dominant clathrin-dependent endocytosis pathway and implies distinct steady-states 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 steady-states 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 clear-cut, 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 [1][2][3][4]. 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 AP-2 [6,7]. Via direct or indirect binding, adaptors recruit the receptor to special membrane regions which are characterized by a particular composition of cage-proteins 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 [10][11][12]. 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].
A study addressing the sorting between Clathrin-vs lipid raft/Caveolae-mediated 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 EGF-concentration 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).
From a modeling point of view, the behaviour of the clathrin-independent 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 [22][23][24], an effect that can be enhanced by consecutive arrangement in the form of cascades [25][26][27][28][29] which has also been validated experimentally [30].
CDE and CIE pathways of EGF receptor 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 CIE-adaptors, respectively. See list of abbreviations.
Other models of ultrasensitivity have been derived for Michaelis-Menten 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 zero-order 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 [34][35][36][37][38][39][40][41][42]. 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 (mono-Ub) of the EGF Receptor could only be observed at high EGF concentrations, raising the question whether mono-Ub 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,[45][46][47][48][49].
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 switch-like response of CIE may result simply from a saturation effect of the CDE pathway. Together with the observation of EGF-concentration dependence of CIE, this analysis invites attention to an ultrasensitive regulatory mechanism for endocytic sorting. We give an analytical derivation of the switch-effect 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 post-ligand 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 clathrin-dependent 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 clathrin-dependent or -independent endocytosis.
In order to simulate the entry of ligand-bound receptor into an endocytosis pathway, we introduced variables CDE and CIE, representing adaptors for clathrin-dependent 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 cage-proteins. We assume that the affinity of activated receptors R_EGF is significantly higher for the CDE-pathway (k cde Ŭ k cie ).
To quantify the fraction of receptor going either pathway, we introduced variables R i_cde and R i_cie . The steady-state 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].

CIE-internalization depends on number of receptors and strength of EGF-stimulation
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 internalization-behavior of R_EGF, the activated receptor. They represent assumptions on the relative quantities of EGF-molecules, unbound receptors and endocytosis adaptors.
The first class of initial values (Figure 2A) represents the case that either the initial number of unbound receptors (R 0 ) or EGF-molecules (EGF 0 ) is lower than the capacity of the CDE-pathway (CDE 0 ) (Generally, X 0 denotes the initial value of variable X). This corresponds to an experimental setting where cells are stimulated with low EGF- concentrations, 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 switch-effect of CIEinternalization
The simulations shown in Figure 2 suggest conditions on the receptor number under which an EGF-dependent switch of CIE-internalization will occur. If a cell possesses less receptors than CDE-adaptors, CIE-internalization will be low independent of EGF-stimulation (cf. Figure 2A). If the cell exhibits more receptors than adaptors for CDE, but less than for both pathways, then, for EGF-stimulations 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, CIE-internalization will be switched on equally strong as CDE-internalization for EGF-stimulations 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 CIE-adaptors 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 .

Steady-State behaviour of internalized receptors
As predicted, for receptor levels lower than CDE 0 ( Figure  3A), CIE-internalization stays close to zero, independently of EGF-stimulation. If the initial amount of receptors is greater than the capacity of CDE, CIEinternalization sets in abruptly, albeit to a moderate degree compared to CDE-internalization, for EGF-stimulations 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 CDE-internalization ( Figure 3C).
We have thus derived an ultrasensitive response of CIEinternalization with respect to EGF-stimulation, 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 CDE-pathway (CDE 0 < R 0 , moderate switch) or both pathways (CDE 0 + CIE 0 < R 0 , maximal switch).

Correspondence to distinct classes of steady-state
For dynamical systems with multiple steady-states, a certain steady-state will be reached depending on whether the system starts in the corresponding basin of attraction [23,51]. Thus, a switch between steady-states 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 switch-effect of CIE-internalization corresponds to such a transition between distinct steady-states of the system. Analytically, one derives two classes of steady-states (see Methods for derivation of conditions): where X* denotes the steady-state concentration of the respective component.
Note that classes of steady-states 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. Steady-state 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 steady-state 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 steady-state revealed initial conditions under which each steady-state class is reached. We found that if both EGF-stimulation 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 steady-state class II. Otherwise, steady-state class I will be reached. In this case, if EGF-stimulation 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 steady-state value of ligand-bound receptor (R_EGF*) is plotted as a function of EGF-stimulation 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 steady-state class I independently of EGF-stimulation, as seen from R_EGF* = 0.
Applying these derived conditions on the initial values, we can also show that the steady-states are stable. A steady-state is stable, if, for small perturbations, the system returns to this steady-state. Consider steady-state 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 steady-state, 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 steady-state class I. An analogous argument can be used to show stability of steady-state 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 steady-state 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 EGF-stimulations, corresponds to a transition of the system between steady-state classes I and II (see Table 1).

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 multisite-modification [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 Hill-type reaction (see Methods). Figure 5 shows the reaction velocity V of the Hill-formula, compared to -production in our model (stimulus-response curve) as a function of EGFstimulation. To generate the stimulus-response curve, we chose the same parameter set as for Fig. 3C as a reference. From this curve we extracted the Hill-coefficienth, V max and K m to compute the corresponding Hill-curve, which will be used as a reference curve later on. The Hill-coeffcient 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]. Stimulus-response curves with Hill-coefficients of 5 or higher are generally considered ultrasensitive [25,26,28]. The Hill-coefficient obtained for the stimulus-response curve shown in Figure 5 is 7.5.

R i_cie *
Approximability of switch-effect by Hill-curve   Table 1

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 clear-cut 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 CIE-adaptors 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 Hill-coefficient h as a measure of steepness (see Methods). The Hill-coefficients 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 stimulus-response curve remains approximable by a Hill-curve.
In Figure 8A we plotted the steady-state values of R i_cde and R i_cie (V max values of the stimulus-response curves) as a function of initial adaptor values CDE 0 and CIE 0 . Here, R 0 and EGF 0 were chosen 1.5.

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 line- The threshold of CIE-internalization (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 well-known that ligand-induced receptor modifications in the form of phosphorylation and/or ubiquitination play a functional role in signaling and contribute to the specificity of adaptor-binding. 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 Robustness of switch-effect for parameter variations Figure 6 Robustness of switch-effect for parameter variations. Comparison between stimulus-response curve (with Hill-coefficient h) and corresponding Hill-curve for selected r k values. Here, k cde = 1.0 and k cie was varied. Plotted are steady-state values of R i_cie (blue) or reaction velocity V (Hill-kinetics, magenta). Here, k onCbl and k offCbl are the rate constants for the association and dissociation of Cbl to ligand-bound 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 CIE-internalization to increasing EGF-stimulation. In Figure 9 we plotted the steady-state 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.

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, EGF-concentration 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 ligand- Deviation from Hill-curve Figure 7 Deviation from Hill-curve. Mean deviation from reference Hill-curve and Hill-coefficient 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.

Discussion
Our analysis addresses the experimentally observed dependence of endocytosis on EGF-concentration [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 steady-states 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 ubiquitin-ligase Cbl during the ubiquitination reaction has been proposed to be necessary for the observed switch-effect 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 switch-effect. Importantly, it is pointed out how by varying the abundance levels of active receptors or endocytic adaptors, cells may modulate their response to incoming EGF-stimulations: 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).

Inclusion of Receptor Modifications
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,[18][19][20][21]58]. Generally, the problem is considered at the 'single-molecule-level' : a single receptor is envisioned, which is thought to enter either the CDE-or the CIE-route (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 CDE-pathway qualitatively alters the system behaviour by enforcing an alternative steadystate (Fig. 3, Table 1). Our model states that an abrupt, switch-like 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 switch-effect '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 CDE-adaptors. 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 clathrin-independent endocytosis upon saturation of the clathrin pathway might have been proposed without mathematical modeling, the unexpected steepness of this switching-behaviour 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 systems-level regulation of endocytic sorting. The natural extensibility of the model to more than two cases may prove applicable in other signaling contexts.

Numerical Simulations and Steady-State Analysis
To numerically solve the systems of equations (1.1 -1.7) we used the MATLAB ODE15s function. The existence of distinct classes of steady-states 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 steady-state, we simultaneously set these equations equal to zero and solved for the variables.

Approximability by Hill-curve
In order to assess the steepness of the derived switcheffect, we compared the obtained stimulus-response curve to the switch-effect described by the Hill-formula. 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 half-maximal velocity is reached and h, the Hill-coefficient, represents the steepness of the response.
In order to investigate the sensitivity of the switch-effect of CIE-internalization 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 stimulus-response curve and extracted the Hill-coefficient 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 stimulus-response curve by the Hill-equation, we computed the mean deviation between the obtained stimulus-response curves and a reference Hill-curve. In this way we formally determined a broad range of parameter values for which an ultrasensitive response of CIE-internalization occurs ( Fig. 5 and Fig.  6).