Exemplary Final Term Paper 5

INTRODUCTION

Background and Hypothesis

The paper, Cell Surface Receptors for Signal Transduction and Ligand Transport: A Design Principles Study by Harish Shankaran, Haluk Resat, and H. Steven Wiley, takes a mathematical modeling approach to understanding the dynamics of cell receptor systems. For their analysis, the authors selected four unique cellular receptor systems for which kinetic parameters were well known and had been measured in previous studies. The selected systems are EGFR, epidermal growth factor receptor, LDLR, low-density lipoprotein receptor, TfR, transferrin receptor, and VtgR, vitellogenin receptor. EGFR is a signaling receptor system that is often seen to be overactive in cancers, and consequently inhibiting this receptor system can be beneficial in slowing or even stopping cancer growth. LDLR is a cargo receptor system primarily linked to the transport of cholesterol. TfR is also a cargo receptor system except its primary function is the cellular uptake of iron. Finally, the VtgR system acts as an importer of Vitellogenin in oocytes, female egg cells in oviparous animals. After being internalized, the Vitellogenin molecule is broken up to create various nutrient molecules for the prospective embryo. These receptor systems perform their designated operations through the process of receptor mediated endocytosis which is depicted in the following diagram.

Figure 1. Receptor Mediated Endocytosis
1. Receptors (blue) bond with ligands (yellow).
2. A deep pocket forms in the cell wall, surrounding the receptor-ligand complexes.
3. The pocket pinches together, successfully internalizing the complexes.
4. Receptors and ligands dissociate, and the internalized cell wall begins to pinch together separating the two.
5. Receptors and ligands are separated from eachother. Receptors are sent back to the cell surface while the ligands are sent to the lysosomes, Golgi complex, or nucleus.
6. Pocket containting now free receptors merges with the cell wall, returning the receptors to the surface.

Due to the general commonality of the process of receptor mediated endocytosis, the same mathematical model can be applied to all four selected receptor systems, provided that the correct kinetic parameters are used for each. If nondimensionalized, the differential modeling equations can be expressed in terms of two new parameters, the specific avidity and the partition coefficient. These two new parameters are functions of the original kinetic parameters defined for the modeled process. Specific avidity describes the system's ability to form receptor ligand complexes while the partition coefficient is a measure of consumption, the system's ability to internalize receptor-ligand complexes. The paper then hypothesizes that, by parameter manipulation, the selected receptor systems may be classified as avidity sensitive, consumption sensitive, or a combination of the two which will, in turn, allow them to partially, if not completely, distinguish the selected receptor systems from one another. In extension to the original hypothesis, we hypothesize that other methods of varying the dimensionless parameters will produce qualitatively similar results, and that this approach to receptor analysis can viably be applied to other receptor systems. Through equilibrium analysis, we also intend to confirm our belief that the system can only have a single, stable equilibrium based on the irreversible directionality of the process of receptor mediated endocytosis as described in this model.

Significance

The authors undertook the analysis of these selected systems, first to better understand how, at the system level, a cell can adjust its kinetic parameters to control function, and secondly, they wished to highlight this concept of using a systems approach rather than a molecular approach to understanding cell function, a method they predict can have more diverse applications within the field of systems biology. Additionally, the results presented in the paper concerning identification of receptor systems through responsiveness to parameter variation could potentially be a useful tool when distinguishing cellular receptor systems. Furthermore, receptor systems play key roles in cell to cell communications such as those used for targeting immune responses, and also are commonly exploited as a means of entry into the cell by harmful viruses. Therefore, any gains in understanding of receptor systems are important steps towards understanding and better controlling or protecting our own biological function.

Brief Results and Model Extension

Based on our mathematical modeling analysis, we identified the EGFR system to express dual sensitivity to variations in both consumption and avidity, the VtgR system to be sensitive to consumption only, and the TfR and LDLR systems to be responsive to avidity only. This confirms the original paper's belief that parameter variation would allow the receptor systems to be classified and at least partially distinguished from one another. Expanding this analysis, we then applied the same parameter variations to two additional signaling receptor systems since the original selection had included three transport and only one signaling. The first was the Erythropoietin receptor system (EpoR) which controls red cell maturation. The second was the Granulocyte-colony stimulating factor receptor system (GCSFR) which promotes white blood cell production. Interestingly, our results show these two additional signaling receptors to exhibit the same dual sensitivity as EGFR, the original signaling receptor. Next, we examined an additional parameter pathway for varying the avidity, and we found results that qualitatively matched those from the initial pathway examined in the original model. Finally, our equilibrium analysis identified the system's single equilibrium as a stable sink equilibrium, all initial conditions eventually approach this equilibrium point in time.

MODEL DESCRIPTION

State Variables and Parameters

State Variables

The three state variables of the receptor mediated endocytosis model examined by Shankaran et. al. are the ligand, L[t], which bonds to the receptors, R[t], to form receptor-ligand complexes, C[t], which can subsequently be internalized by the cell.

Table 1: State Variables

Variable	Description	Initial Value
R	Number of Free Surface Receptors (molecules)	RT
C	Number of Receptor-Ligand Complexes (molecules)	0
L	Extracellular Concentration of the Ligand (nM)	0.01 KD

Parameters

The two steps of bonding and internalization can be seen to be analogous to chemical reactions and are therefore governed by corresponding kinetic rate constants, k_on, k_off, k_e, k_t. To create a functional model, the authors introduced the extracellular volume, V, a rate of introduction of new receptors on the cell surface, Q_R, and a time dependent function f[t] that could be used to describe the rate of introduction of new ligand to the system.

Table 2: Independent Parameters

Variable	Description	EGFR Value	TfR Value	LDLR Value	VtgR Value	EpoR Value	GCSFR Value
kon	Rate of formation of receptor-ligand complexes (/(min * nM))	0.097	0.0030	0.0028	5.4 E -5	0.48	0.111
koff	Rate of dissociation of or receptor-ligand complexes (/min)	0.24	0.09	0.04	0.07	0.01	0.03
ke	Rate of internalization of receptor-ligand complexes (/min)	0.15	0.6	0.195	0.108	0.06	0.104
kt	Rate of internalization of unbonded receptors (/min)	0.02	0.6	0.195	0.108	0.005	0.005
QR	Rate of creation of new receptors on the cell surface (molecules/min)	4000	15600	390	2.2 E10	9.9	25
V	Extracellular volume (liters)	4 E-10	4 E-10	4 E-10	4 E-10	4 E-10	4 E-10
f[t]	Introduction of new ligand molecules to the control volume (nM/min)	varies

Additional parameters can be calculated as functions of the original six independent parameters defined for each of the systems. These calculated parameters can help with physical understanding of the receptor systems.

Table 3: Calculated Parameters

Variable	Formula	Description	EGFR Value	TfR Value	LDLR Value	VtgR Value	EpoR Value	GCSFR Value
KD	koff/ kon	Dissociation constant (nM)	2.47	29.8	14.3	1300	0.021	0.27
RT	QR/ kt	Steady state receptor abundance (molecules)	2 E 5	2.6 E 4	2 E 4	2 E 11	2000	5000
D	ke/ kt	Degree of Downregulation	7.5	1.0	1.0	1.0	12	21
β	ke/ koff	Partition coefficient	0.63	6.67	5.51	1.44	6	3.5
$\gamma$	(QR kon)(109) /(V kt koff Nav)	Specific avididty	0.34	0.004	0.006	638.6	0.39	0.077

• • N_av = Avogadro's number

Model Assumptions

Biological systems are inherently complex, and therefore it was necessary for the authors to make a series of simplifying assumptions regarding receptor mediated endocytosis that would allow them to create a system that was manageable enough for analysis. These assumptions are as follows:

• Q_R includes both the rate of creation of new receptors along with accounting for old receptors being recycled.
• The cell is sufficiently large that cell wall surface area losses do not factor into limiting the the rate of internalization. The rate of consumption of cell wall through internalization roughly equals the rate of return at the end of the endocytosis process.
• The extracellular volume is always well-mixed, i.e. ligand concentration gradients caused by internalization are negligible.
• Processes after internalization are ignored.
• Once internalized, the ligand is consumed and never returns to the extracellular volume.
• Formation of receptor ligand-complexes occures in an even one to one ratio.
• Receptors are assumed to be uniformly distributed such that all have roughly the same probability of bonding with a ligand.
• The system behaves linearly for small input concentrations of ligand.

Model Equations

The equations used to describe the systems are derived from a series of chemical reaction formulas that have been created from a picture provided within the paper. This picture depicts the situation that the authors set out to model, and includes the independent parameters used for the model.

Figure 1. Visual Representation of the System to be Modeled

Table 4: Reaction Formulas

	Rate Constant	Formula
1	kon	$R + L \quad \longrightarrow \quad C$
2	koff	$C \quad \longrightarrow \quad R + L$
3	kt	$R \quad \longrightarrow \quad R_{internalized}$
4	ke	$C \quad \longrightarrow \quad C_{internalized}$

Equation 1:

$\frac{dR}{dt} = -k_{on} R(t) L(t) + k_{off} C(t) - k_t R + Q_R$

Chemical rate laws determine the rate of a given chemical reaction. The rate equals a constant multiplied by the concentrations of the reactants raised respectively to the power of their stoichiometric coefficient. Looking at Table 4, we see that the first term comes from the chemical rate law for the first formula and is negative corresponding to the reduction of R. We note that while the state variables are defined in number not concentration, the corresponding rate constant has been modified accordingly and the units work out appropriately. The positive second term then corresponds to the second formula, and the third term to the third formula. Finally, the fourth term, Q_R, describes the rate of creation of new receptors on the cell surface and therefore stands alone.

Equation 2:

$\frac{dC}{dt} = k_{on} R(t) L(t) - k_{off} C(t)- k_e C(t)$

The terms of Equation 2 are derived in the same manner as the first three terms of equation one, coming from formulas one, two, and four respectively.

Equation 3:

$\frac{dL}{dt} = \frac{-k_{on} R(t) L(t) + k_{off} C(t)}{N_{av} V / 10^{-9}} + f(t)$

First looking at the numerator of the first term, we see that once again the formulas were used with their corresponding rate constant to create both of the numerator terms. However, unlike R and C which have units of molecules, L has units of nM, and thus the denominator, Avogadro's number times volume, exists to convert molecules to M. Then, in what appears to be a slight discrepancy from the original article, we have introduced the 10^-9 factor in the denominator to convert M to nM, which are the units the article gives for L. We incorporated this conversion into the definition of $\gamma$ as well.

Nondimensionalization

Table 5: Nondimensionalization Definitions

New Term $\quad$	Definition $\quad$
t*	$\quad = k_{off} t$
R*	$\quad= \frac {k_t R}{Q_R}$
C*	$\quad = \frac {k_t C}{Q_R}$
L*	$\quad = \frac {k_{on} L}{k_{off}}$

Now, using these definitions that have been defined within the original paper, the initial differential equations can be nondimensionalized by observing what parameters the given derivative, $\frac{dY}{dt}$ must be multiplied by to transform it to a dimensionless form, $\frac{dY*}{dt*}$. Subsequent simplification will allow us to attain the dimensionless equations that will be used when examining the model.

Equation 1*:

$[\frac{dR}{dt} = -k_{on} R(t) L(t) + k_{off} C(t) - k_t R + Q_R]\frac {k_t}{Q_R k_{off}}$

$\frac{dR*}{dt*} = [-k_{on} R(t) L(t) + k_{off} C(t) - k_t R + Q_R]\frac {k_t}{Q_R k_{off}}$

$\frac{dR*}{dt*} = -(\frac{k_{on}}{k_{off}})(\frac{k_{t}}{Q_{R}}) R(t) L(t) + (\frac{k_{off}}{k_{off}})(\frac{k_{t}}{Q_{R}})C(t) - (\frac{k_{t}}{k_{off}})(\frac{k_{t}}{Q_{R}}) R + (\frac{k_{t}}{k_{off}})(\frac{Q_{R}}{Q_{R}})$

$\frac{dR*}{dt*} = - R^*(t) L^*(t) + C^*(t) - (\frac{k_{t}}{k_{off}})(R^*(t) - 1)$

Equation 2*:

$\frac{dC}{dt} = k_{on} R(t) L(t) - k_{off} C(t)- k_e C(t)$

$[\frac{dC}{dt} = k_{on} R(t) L(t) - k_{off} C(t)- k_e C(t)]\frac {k_t}{Q_R k_{off}}$

$\frac{dC*}{dt*} = [k_{on} R(t) L(t) - k_{off} C(t)- k_e C(t)]\frac {k_t}{Q_R k_{off}}$

$\frac{dC*}{dt*} = (\frac{k_{on}}{k_{off}})(\frac{k_{t}}{Q_{R}}) R(t) L(t) - (\frac{k_{off}}{k_{off}})(\frac{k_{t}}{Q_{R}})C(t) - (\frac{k_{e}}{k_{off}})(\frac{k_{t}}{Q_{R}})C(t)$

$\frac{dC*}{dt*} = R^*(t) L^*(t) - C^*(t) - (\frac{k_{e}}{k_{off}})C^*(t)$

Note: $\beta = \frac{k_{e}}{k_{off}}$

$\frac{dC*}{dt*} = R^*(t) L^*(t) - C^*(t) - \beta C^*(t)$

Equation 3*:

$[\frac{dL}{dt} = \frac{-k_{on} R(t) L(t) + k_{off} C(t)}{N_{av} V \ 10^{-9} } + f(t)](\frac{k_{on}}{k^2_{off}})$

$\frac{dL*}{dt*} = [\frac{-k_{on} R(t) L(t) + k_{off} C(t)}{N_{av} V \ 10^{-9}} + f(t)](\frac{k_{on}}{k^2_{off}})$

$\frac{dL*}{dt*} = (\frac{1}{N_{av} V \ 10^{-9}})[-(\frac{k^2_{on}}{k^2_{off}})R(t) L(t) + (\frac{k_{on}}{k_{off}}) C(t)] + (\frac{k_{on}}{k^2_{off}})f(t)$

Note: Here $\quad R = \frac {R^* Q_R}{k_t} , \quad C = \frac {C^* Q_R}{k_t} , \quad L = \frac {L^* k_{off}}{k_{on}} \quad$ are substituted in before the equation is further simplified.

$\frac{dL*}{dt*} = (\frac{1}{N_{av}\ V \ 10^{-9}})[-(\frac{k_{on}}{k_{off}})(\frac{Q_{R}}{k_{t}})R^*(t) L^*(t) + (\frac{k_{on}\ Q_{R}}{k_t\ k_{off}}) C^*(t)] + (\frac{k_{on}}{k^2_{off}})f(t)$

Note: $\gamma = \frac{k_{on}\ Q_R}{N_{av}\ k_{off}\ k_t\ V\ 10^{-9}}$

$\frac{dL*}{dt*} = \gamma [-R^*(t) L^*(t) + C^*(t)] + (\frac{k_{on}}{k^2_{off}})f(t)$

Linearization

Using partial derivatives and first term Taylor expansion linearization allows for the following formula for the linearization of multivariable equations at a given point, p.

$f({\mathbf{x}}) \approx f({\mathbf{p}}) + \left. {\nabla f} \right|_{\mathbf{p}} \cdot ({\mathbf{x}} - {\mathbf{p}})$

From Table 1, the initial condition for $\frac{dL*}{dt*}$ is 0.01 when the proper Table 5 transformation is performed. This value is relatively close to 0 which, in turn, indicates that a reasonable approximation of the system can be obtained through the multivariate linearization about the point (R*_o = 1, C*_o = 0, L*_o = 0 ).

This linearization results in:

Equation 1*L:

$(\frac{dR*}{dt*})_{lin} = (-1*0 + 0 - \frac{k_t}{k_{off}} * (1-1)) + (-L^*-\frac{k_t}{k_{off}})|_{(1,0,0)} * (R^* - 1) + (1)|_{(1,0,0)} * (C^* - 0) - R^*|_{(1,0,0)} * (L^* - 0)$

$(\frac{dR*}{dt*})_{lin} = -\frac{k_t}{k_{off}}(R^*(t) - 1) + C^*(t) - L^*(t)$

Equation 2*L:

$(\frac{dC*}{dt*})_{lin} = (1*0 + 0 + \beta * 0) + L^*|_{(1,0,0)} * (R^* - 1) + (-1 - \beta)|_{(1,0,0)} * (C^* - 0) + R^*|_{(1,0,0)} * (L^* - 0)$

$(\frac{dC*}{dt*})_{lin} = -C^*(t) - \beta C^*(t) + L^*(t)$

Equation 3*L:

$(\frac{dL*}{dt*})_{lin} = (\gamma(1*0 + 0) + (\frac{k_{on}}{k^2_{off}})f(t)) + (- \gamma L^*|_{(1,0,0)} * (R^* - 1) + \gamma|_{(1,0,0)} * (C^* - 0) - \gamma R^*|_{(1,0,0)} * (L^* - 0)$

$(\frac{dL*}{dt*})_{lin} = (\frac{k_{on}}{k^2_{off}})f(t) + \gamma C^*(t) - -\gamma L^*(t)$

Significance: The interpolating capabilities of Mathematica are relatively powerful and efficient. However, it would be a particular advantage to be able to obtain symbolic solutions that have been completely decoupled from one another. This would allow for shorter processing time and more straightforward application of mathematics when studying the system and its response to parameter manipulation. Interestingly, the linearization of the nondimensionalized equations has accomplished the first element of decoupling. $(\frac{dC*}{dt*})_{lin}$ and $(\frac{dL*}{dt*})_{lin}$ are now expressed solely in terms of $\beta$, $\gamma$, C*(t), and L*(t). Now it will be possible for Mathematica to solve these two differential equations in symbolic form to obtain solutions in terms of $\beta$, $\gamma$, and t.

Internalized Ligand

The three equations that have been examined and manipulated so far describe the rates of change of the number of free receptors on the cell surface, the number of receptor-ligand complexes on the surface, and the extracellular ligand concentration. However, the authors are primarily interested in the change in efficiency of ligand uptake when certain key parameters are varied. Therefore, the concentration of the internalized ligand with respect to time would be very useful for such analysis. Recognizing this, the authors created a fourth differential equation to describe the rate of change of the concentration of the ligand that had been successfully internalized.

Equation 4.

$\frac{dL_i}{dt} = \frac{k_e * C(t)}{N_{av} * V}$

Note: Initial condition, (L_i)_o = 0

Looking back at Figure 1 and Table 4, the numerator of this rate equation comes from reaction 4 which describes the internalization of receptor-ligand complexes and thus the internalization of the ligand itself. The denominator acts in the same way as in equation 3 which is to convert a number of internalized ligands to a concentration. This is advantageous when assessing how much of the original extracellular ligand has been internalized at any given time. We remember that such comparative assessments only only hold significance when there is no introduction of new ligand into the extracellular volume with time.

Simulation

Figure 2 Recreation

A: Using the parameter values for the EGFR system, Mathematica's NDSolve function will be used to perform numerical integration and derive an interpolating function that describes C(t) for the desired range of t. The intial conditions are defined in Table 1, and f(t), the function defining the rate of introduction of new ligand molecules to the system, is set to 0. The interpolating function will then be multiplied by k_t and divided by Q_R to nondimensionlize it to C*(t). This is the transformation that has been defined in Table 5. For scaling reasons the function is then multiplied by 100 and finally plotted with respect to time. These steps are repeated using values of k_e that correspond to downregulation values, D, of 1, 2, 5, 7.5, 10, and 20.

B: The same general methodology is used as in A, except now f(t) has been defined as an impulse function of magnitude $\frac{0.5 * Q_R}{N_{av} * V}$ with lengths of both the pulse and the resting periods set to 360 min. This function is taken from Professor Chiel's text in Chapter 7. Other differences are that the functions are no longer multiplied by 100 for scaling. Instead, each curve is normalized by dividing it by its individual maximum value on the plotted time interval.

C: The same process as in B is followed except only for two downregulation curves. Also, the characteristics, such as amplitude and lengths of both the pulse and the resting periods, were randomized using a normal distribution in which the mean values were set to 400 min . For our recreation, estimation was used to create a pulse function that closely resembled the function presented in the text.

Figure 3 Recreation

Figure 3 from the text is primarily a comparative examination of the four selected systems for the impact of downregulation on their efficiency of ligand uptake. In order to judge efficiency of uptake, the amount of ligand in the system must be considered constant. Therefore f(t) must be set to zero, signifying that no new ligand is being introduced with time. Now, Equation 4 can be numerically integrated using NDSolve and the parameter values unique to each system. Varying k_e provides interpolating functions for the each of the desired downregulation values. The functions are then divided by the initial condition of L(t) and multiplied by a 100 so that they represent the percentage of ligand that has been internalized at given time. Plotting the curves completes the creation of Figure 3.

Figure 4 Recreation

Figure 4 is created in exactly the same manner as Figure 3 except that now the extracellular volume, V, is varied instead of k_e.

Figure 5 Recreation

A: First, using the linearized equations, Mathematica's DSolve function will be used to obtain a symbolic expression for C*(t*) in terms of $\beta$, $\gamma$, and t. In order to achieve this symbolic expression, the initial conditions L*(t*)_o = 0.01 and C*(t*)_o = 0 are used, and then the solution is multiplied by 100 for additional simplification which will have no impact on further calculations. $\tau$ is defined as the relaxation time of C*(t*), and equals the time at which C*(t*) falls to a value of $\frac{C*_{max}}{e}$. The value of $\tau$ can be obtained by differentiating the symbolic equation, solving for its root, plugging the root in the the original equation to define C*_max, and then using Mathematica's FindRoot function to find $\tau$. This method of finding $\tau$ can now be used to create the contour plot that is Figure 5a.

B: r_s is defined to be a system's relative sensitivity of relaxation time with respect to changes in $\beta$ and $\gamma$. Therefore, a value of r_s can be obtained by calculating $\tau$ at a given point and then increasing $\beta$ and $\gamma$ by 1% separately to calculate two new relaxation time values. These values can then be used to calculate numerical partial derivatives for the mathematical formulas used to define r_s. The original paper defines:

$r_s = \frac{\partial_{\gamma}\tau}{\partial_{\beta}\tau}$

$\partial_{\gamma}\tau = \frac{\gamma}{\tau} * \frac{\partial \tau}{\partial g} * \frac{1}{Ln[10]}$

$\partial_{\beta}\tau = \frac{\beta}{\tau} * \frac{\partial \tau}{\partial b} * \frac{1}{Ln[10]}$

Note: b and g are simply indices that are used to create a log log scaling for the plot. b = Log[$\beta$] and g = Log[$\gamma$].

With r_s defined, Mathematica, can now easily create a contour plot.

Figure 6 Recreation

Figure six in the original paper is a schematic that has been created using figure 5a along with a series of boundary conditions to divide the plotted range into regions based on the sensitivity expected to be displayed by a receptor with specific avidity and partition coefficient values within that region. We have decided not to recreate the schematic itself, but to instead create another contour plot that will be based on the magnitude of the sensitivity of the model to changes in $\beta$ and $\gamma$. The authors have defined this magnitude to be::

$m_s = \sqrt{\partial_{\beta}\tau^2 + \partial_{\gamma}\tau^2}$

Using the code from figure 5b, Mathmatica can now be used to create a contour plot that will quantify the total robustness of selected receptor systems.

Extension Simulations

The additional systems, GCSFR and EpoR, will be analyzed for responsiveness to both avidity and consumption by applying the same methodology used in the recreation of figures three and four. Then, the responsiveness to an alternate method of varying avidity will be analyzed by varying the internalization rate of unbound receptors, k_t. k_t was selected because it is reasonable to believe that, if the cell can control the internalization rate of complexes, k_e, then it can likely control the rate of free receptors as well. Furthermore, changing k_t, will vary downregulation as was done in Figure 3, but now we expect to see avidity sensitivities instead of consumption sensitivities. This will decouple sensitivity to downregulation and sensitivity to consumption. Finally, for equilibrium analysis, we will use Mathematica to solve for the equilibrium of the dimensionless equations symbolically, and then use partial derivatives to obtain eigenvalues that describe the stability of the equilibrium point. As added level of equilibrium analysis we will set each of the original equations equal to zero, and then solve for C[t] as functions of R[t] and L[t]. This creates a set of three planes where the rate of change of each respective state variable is zero. These are called nullplanes, and points at which they all intersect are equilibrium points where the system is unchanging in time.

RESULTS

Figure 2

Original Results	Our Results
Figure 2. Shankaran et. al.. EGFR system response to varying downregulation and ligand introduction.	Figure 2. Group Simulation. EGFR system response to varying downregulation and ligand introduction.

In Figure 2, the effect of downregulation on the receptor function was studied for the EGFR system. Downregulation is the reduction of a cellular component caused by some stimulus. In this case, the stimulus is an increased rate of internalization of receptor-ligand complexes, quantified as an increase in k_e. Of the four receptor systems studied in this paper, EGFR was selected for this analysis because it is a signaling receptor which makes its ability to quickly respond to ligand concentrations integral to its function. TfR, LDLR, and VtgR are transport receptors, so internalization rates are important but not quite as key to function as with EGFR. Figure 2a serves primarily as a first test of the model to confirm that it is exhibiting logical responses and to get a first look at the impact of downregulation on receptor function. The curves plotted are dimensionless functions for the number of receptor-ligand complexes, C*[t]. In response to a one time impulse of ligand that has been quantified as an initial condition of L[t], we see the number of complexes rise to a peak then fall once again to zero. Logically this makes sense since a sudden introduction of ligand would cause a bunch of complexes to form immediately after introduction, but as more and more of the ligand has been successfully internalized there soon becomes no ligands for the formation of complexes. Thus, the number of complexes degrades to zero in time. Additionally, as downregulation is increased by increasing k_e, we see the height of the maximum peak is reduced, which also makes sense because a higher rate of internalization of complexes would reduce the ligand concentration more quickly as well as preventing large numbers of complexes from building up on the cell surface. Another key impact of increasing downregulation is that the peak occures at an earlier time, indicating a quicker response to the initial ligand input. With these two trends in mind, the rate of introduction of new ligand was then switched to a pulse function with active and resting periods of 360 minutes and an amplitude small enough to keep the vast majority of the receptors unbound on the cell surface, $\frac {0.05 Q_R}{N_{av} * V * 10^{-9}}$. Also, the functions were normalized by dividing each curve by its respective maximum value. This allowed us to look primarily at the efficiency of the responses without being distracted by response magnitudes. Once again, higher downregulation, led to a closer fit to the ligand introduction curve, and thus indicated a quicker response to the ligand stimulus. Then, in Figure 2c, the pulse curve was randomized using a Gaussian distribution. These results serve only to confirm the positive impact of downregulation with respect to response efficiency to perturbations of the extracellular ligand concentration.

Figure 3

Original Results	Our Results
Figure 3. Shankaran et. al. Responsiveness of the selected receptor systems to changes in consumption.	Figure 3. Group Simulation. Responsiveness of the selected receptor systems to changes in consumption.

Remembering the dimensionless system of equations that was attained earlier, we note that $\beta$, the partition coefficient used to quantify a system's dependence on consumption, was defined as the ratio of the rate of complex internalization, k_e, to the rate of dissociation of receptor-ligand complexes, k_off. Pairing this with the study of downregulation from Figure 2, we see that by increasing k_e to increase downregulation we were subsequently increasing the partition coefficient as well, and therefore identifying whether or not the function of the given receptor system, EGFR is this case, was responsive to changes in consumption. Varying downregulation through increasing k_e has now been applied to the remaining systems to determine whether or not they are consumption controlled. One important change has been made in the analysis process however. Because the other receptor systems are transport receptors whose primary function is to internalize their designated ligand, the original authors decided that the percentage of the extracellular ligand that had been internalized at a given time would be a more appropriate response measure than the dimensionless or normalized number of receptor-ligand complexes used in Figure 2. Results of this analysis indicate that EGFR and VtgR systems are consumption controlled while TfR and LDLR systems are robust to changes in consumption. Looking at the parameter values in Table 3, it appears that lower values of $\beta$ indicated sensitivity to its variation. Also, it is interesting to note that, contrary to intuition, EGFR, the sole signaling receptor, did not exhibit the quickest and most efficient responses. In order of increasing efficiency/speed of response the receptors are VtgR, EGFR, TfR, and finally LDLR.

Figure 4

Original Results	Our Results
Figure 4. Shankaran et. al. Responsiveness of the selected receptor systems to changes in avidity.	Figure 4. Group Simulation. Responsiveness of the selected receptor systems to changes in avidity.

Now, the authors needed to address the second element of the hypothesis which discussed the receptor systems' responsiveness to changes in avidity. Specific avidity is a measure of the system's efficiency at forming receptor-ligand complexes, and it is quantified in the dimensionless equations as $\gamma = \frac{Q_R * k_{on} * 10^9}{k_t * k_{off} * V * N_{av}}$. After considering this mathematical definition of specific avidity, the authors decided that varying the extracellular volume would be the most biologically significant method of altering avidity. The upper volume limit, 4*10^-10 L, was calculated by observing that "a typical cell culture experiment would have 10 mL of media and a cell count of 2.5x 10^7 cells." Then, the lower limit, 4*10^-13 L, "corresponded to approximately one cell volume per cell." Results from this analysis reveal LDLR and TfR to be highly sensitive, EGFR to be moderately sensitive, and VtgR showed no discernible variation. Once again looking at the parameter values in Table 3, it appears that, similar to the partition coefficient, higher specific avidity values correspond to systems that are more robust to variation in avidity.

Figure 5

Original Results	Our Results
Figure 5. Shankaran et. al. Efficiency analysis with respect to avidity and consumption values.	Figure 5. Group Simulation. Efficiency analysis with respect to avidity and consumption values.

As receptor systems whose primary function is to internalize their complementary ligand, the time it takes these systems to do so becomes the primary measurement of the efficiency of the system as a whole. To quantify this response time, Shankaran et. al. used a relaxation time, $\tau$, which they defined as the time it takes the number of complexes, C[t], to decay to a value of 1/e times its maximum value. Therefore, an increase in $\tau$ indicates a decrease in efficiency of the system and vise versa. In order to attain an expression for $\tau$ in terms of $\beta$ and $\gamma$, the linearized functions presented in the Model Description were used as a good approximation due to the small initial concentration of ligand that was being studied. Thinking purely from a physical standpoint, since specific avidity describes the system's ability to form receptor-ligand complexes and the partition coeffficient describes the system's ability to internalize those complexes once formed, we would expect an increase in either component to lead to greater efficiency of the entire system, and in fact we do see this trend in Figure 5a. The highest values of $\tau$ are located in the lower left hand corner of the plot and the lowest values are located in the upper right hand corner. As a quick check of Figure 5a, we can estimate the log[$\tau$] value for the EGFR system to be about 1.1 which yields a dimensionless $\tau$ of 10^1.1 which is then divided by k_off as the reverse of the nondimensionalizing of time seen in Table 5 to yield a relaxation time of 52 min for the default EGFR system. This value can now be compared with the Figure 2a curve for a downregulation of 7.5 which is the default EGFR system. The curve has a maximum of roughly 0.45 which means $\tau$ corresponds to the dimensionless C*[t] value of 0.17. The time at which the plot reaches this value is approximately 50, thus demonstrating an agreement between the two figures.

Having seen that increasing either $\beta$ or $\gamma$ leads to a more efficient system, the logical next step is to find out which increase will cause the greatest gain in efficiency for a given receptor system. To do this, Shankaran et. al. defined r_s, the relative sensitivity of $\tau$, as the change in $\tau$ with respect to a 1% increase in $\gamma$ divided by the change in $\tau$ with respect to a 1% increase in $\beta$. Therefore, high values indicate that increasing the specific avidity, $\gamma$, is the best option while low values indicate that increasing consumption, $\beta$, is the best way to increase efficiency. In Figure 5b, we see that the contour plot of r_s has divided the plotted area into three distinct regions. The red region is most responsive to avidity, the blue to consumption, and the green/yellow region indicates relatively equal sensitivities to both.

Figure 5 Linear Approximation Test

Test: Confirmation of the linear approximation used to create Figure 5

To create Figure 5 we used a linearization of the dimensionless equations at the point ( R* = 1 , C* = 0, L* = 0 ) that would allow us to solve symbolically. However, we needed to make sure that this approximation was appropriate, so we set up a Manipulate function in Mathematica that would allow us to compare the symbolic equation for C*[t*] attained from the linearization to the interpolating function Mathematica's NDSolve created from the non-linearized equations. Our results show a near perfect approximation for all but the lowest values of specific avidity, and even then the difference is relatively slight. We conclude that the approximation used for the creation of Figure 5 was appropriate.

Figure 6

Original Results	Our Results
Figure 6. Shankaran et. al. Schematic representing the qualitative results of Figure 5 along with additional robustness analysis otherwise not included in the paper.	Figure 6. Group Simulation. Recreation of the additional robustness analysis used to create Figure 6 from the original paper.

In Figure 5a, we saw that increasing $\beta$ and $\gamma$ lead to more efficient receptor systems, and in Figure 5b we discovered the regions for which increasing $\beta$ had the greatest impact, increasing $\gamma$ had the greatest impact, and where both increases has relatively the same impact. Now, the next step is to quantify the magnitude of these responses to determine the robustness of a system with given $\beta$ and $\gamma$ values. To quantify the magnitude of the robustness of the system, Shankaran et. al. simply added in quadrature the change in $\beta$ with respect to a 1% increase in $\gamma$ and the change in $\gamma$ with respect to 1% increase in $\beta$. This allows for the creation of the contour plot displayed as our result. Low values indicate robustness while high values indicate sensitivity, or relatively large increases in efficiency for increases in $\beta$ or $\gamma$. This plot shows us that while increasing $\beta$ and $\gamma$ does produce gains in efficiency, it does so at a declining rate such that eventually the increases in efficiency will become negligible. Looking at the schematic provided as Figure 6 in the original paper, we note that it is simply a combination of the results of Figures 5a and 5b along with the robustness plot we have provided. Figure 5a is used as the background, Figure 5b is used to mark the avidity controlled, consumption controlled, and dual sensitivity regions, and then our robustness plot is used to mark the robust and sensitive regions of their schematic.

Hypothesis

Remembering from the introduction...

"Specific avidity describes the system's ability to form receptor ligand complexes while the partition coefficient is a measure of consumption, the system's ability to internalized receptor-ligand complexes. The paper then hypothesizes that, by parameter manipulation, the selected receptor systems may be classified as avidity controlled, consumption controlled, or a combination of the two which will in turn all them to be partially, if not completely, distinguished from one another."

With the figures originally presented within the paper that we have now recreated and confirmed, this hypothesis has been addressed. First, the results in Figure 2 confirmed a working and reasonable model. Next, in Figure 3, we varied the internalization rate of receptor-ligand complexes, k_e, which in turn varied the partition coefficient $\beta = k_e / k_{off}$. The corresponding plots have shown the EGFR and VtgR systems to be particularly sensitive to changes in the partition coefficient, while TfR and LDLR are generally unresponsive. Then, in Figure 4, we varied the extracellular volume which caused variation in the specific avidity, $\gamma = \frac{Q_R * k_{on} * 10^9}{V * k_t * k_{off}}$. This time, TfR and LDLR were both very responsive, EGFR was reasonably responsive, and VtgR showed no change. Now we see that EGFR has been isolated as responsive to both avidity and consumption, Vtgr has been isolated as responsive only to consumption, and TfR and LDLR were unable to be separated from one another, both displaying responsiveness to only changes in avidity. Additionally, TfR and LDLR even expressed nearly identical magnitudes of responsiveness to varying $\beta$ and $\gamma$. In any case, we have still managed to identify each system as avidity controlled, consumption controlled, or both as well as using these identifications to at least partially separate the receptor systems from one another based on qualitative differences.

Discrepancies

Correction 1

From the original text:

"For computing these sensitivity indices, we generated a linear grid for the variables (b, g). Logarithmic b and c values at a particular (b_i, g_j) grid point are given by $\beta_i = log_{10}(b_i)$ and $\gamma_j = log_{10}(g_j)$. The i and j indices vary from 1 to 100; that is, the size of the grid was 100x100. The linear grid spanned the values from -2 to 2 and -3 to 3 for the b and g variables, respectively."

This description is given in the Materials and Methods section of the original text. It is describing a portion of the methodology used to create Figure 5a. However, a mistake has been made in how b_i and g_j are defined with respect to $\beta$ and $\gamma$. Looking at Figure 5a, we see that the contour plot is a log log plot of $\beta$ and $\gamma$. Therefore, the arguments of the logarithms in the definitions above should be $\beta$ and $\gamma$ respectively instead of bi and g_j. This mistake is made more obvious by the ranges given for b_i and g_j, "-2 to 2 and -3 to 3," because the log of a negative number is undefined. Our recreations above were made using the correct definitions of $bi = log_{10}(\beta_i)$ and $gj = log_{10}(\gamma_j)$.

Correction 2

Our first attempt at recreating Figure 5b resulted in a contour plot that was scaled appropriately, but was in fact a 180 degree rotation of the correct contour plot shown in the paper.

Incorrect Figure 5b: Rotated version of Figure 5b caused by a discrepancy in original text.

After carefully checking our work and still being unable to find the discrepancy between our code and the mathematical definitions given in the paper, Professor Chiel suggested that we contact the authors of the original paper to see if they had any insight as to what we were doing wrong. A special thanks is due to the original authors of the paper for their quick response, within one business day. In his response to our email, Harish Shankaran pointed out that there had been a typo in the original paper.

Email response from Harish Shankaran describing discrepancy

With this typo corrected, we were then able to successfully recreate Figure 5b as well as the robustness plot we created for Figure 6.

Extension

Additional Receptor Systems

Figure 7: Receptor Analysis for two Additional Signaling Receptors

The additional analysis of new signaling receptor systems yielded some very interesting results. The original paper only sought to establish that the selected receptor systems could in fact be categorized based on sensitivities to avidity and consumption. They did not however look for any trends in the results of their analysis. With the addition of two more signaling receptor systems, EpoR and GCSFR, we have now uncovered the beginnings of a meaningful trend. In the figure above we see that both GCSFR and EpoR exhibit the same dual sensitivity as we saw with the EGFR system. This leads us to believe that perhaps dual sensitivity is a trait of most if not all signaling receptor systems.

Alternate Method of Varying Avidity

Figure 8: Results from varying downregulation by changing k_t instead of k_e

The rate of internalization of free receptors, k_t, appears in the denominator of the expressions for both specific avidity and downregulation (See Table 3). Therefore, increasing downregulation would require a decrease in k_t which would cause a subsequent increase in specific avidity. In this segment of our extension, we expected to observe the same systems displaying sensitivity to avidity as in Figure 4. We were happy to see that our results clearly match those of Figure 4 with LDLR and TfR showing the greatest degrees of sensitivity, EGFR a more moderate sensitivity, and VtgR no perceptible sensitivity.

Equilibrium Analysis

Using the algebra computing capabilities of Mathematica along with the dimensionless equations obtained earlier, we calculated the eigenvalues of the system. Eigenvalues are related to straight line solutions to the differential equations, and a negative value indicates a progression towards the equilibrium with time while a positive value indicates a progression away from the equilibrium. For our analysis, the introduction of new ligand, f[t], was set to zero first because it is set to zero for all figures except 2b and 2c, but also because Mathematica was unable to calculate the eigenvalues symbolically when f[t] had value. Mathematica output the following eigenvalues;

$\lambda_1 = - k_t/ k_{off}$

$\lambda_2 = -(1+\beta+\gamma) - \sqrt{(1+\beta)^2- 2\beta \gamma + \gamma^2 + 2}$

$\lambda_3 = -(1+\beta+\gamma) + \sqrt{(1+\beta)^2- 2\beta \gamma + \gamma^2 + 2}$

Remembering that all parameters examined in this model are positive, we immediately see that $\lambda_1$ will be negative. Looking just slightly more deeply shows that $\lambda_2$ will at least have a negative real part since a square root cannot output a negative value. Then we get to $\lambda_3$ whose sign is a bit less obvious. In order to help us figure out the sign we expanded the argument of the square root, and we squared and expanded the portion of the number that is outside of the square root.

Squared outside: $1 + 2 \beta + \beta^2 + 2 \gamma + 2\beta \gamma + \gamma^2$

Expanded square root argument: $1 + 2 \beta + \beta^2 + 2 \gamma - 2\beta \gamma + \gamma^2$

Now, we can see that because of the difference in sign of the $2\beta \gamma$ term between the two expressions, $\lambda_3$ must be negative as well. We know this because the argument of the square root is less than the square of what is outside of the square root which means that it will still be less after the square root is taken and hence not of large enough magnitude to make $-(1+\beta+\gamma)$ positive. Since we now know that all of our eigenvalues are negative, we can classify this equilibrium as sink which means that all sets of initial conditions for R[t], L[t] and C[t] will eventually progress to a steady state value at this equilibrium.

Next, we sought to locate the equilibrium point for the four receptor systems originally selected for study. As earlier described, we accomplished this by setting each of the original differential equations to zero and then solving for C'[t] to obtain the equations for the null planes of the system. These planes were then plotted, and their intersection was located as the sink equilibrium that we had just identified. We would also like to note that Mathematica was used to symbolically solve the the system of original differential equations to confirm that there was in fact only one equilibrium point for the system.

Figure 9: Equilibrium Identification using NullPlanes

Thinking logically, the results of this equilibrium identification make perfect sense. As modeled, there is only one direction for the system to go. Any extracellular ligand is essentially trapped and must eventually bind with a receptor and be internalized. Once internalized, the ligand has been effectively removed from the system and never returns to the extracellular volume. Therefore the L[t] must go to zero in time and consequently take C[t] to zero along with it. Now looking at the equilibrium value of R[t] we remember the calculated parameter R_T. R_T is defined as the steady state receptor abundance and hence we see that its value matches up with the R[t] value we see at the nullplane identified equilibrium points.

DISCUSSION

Hypothesis Analysis

There are two dimensionless parameters that play a major role in the function of both signaling and transport cell surface receptors. The first is specific avidity which is used to quantify the receptor system's ability to form receptor-ligand complexes. The next is a measure of consumption, the partition coefficient, which quantifies a receptor system's ability to internalize receptor-ligand complexes before they dissociate. The analysis undertaken in the original paper that has been recreated here sought understanding of how these two functional parameters could be used to predict receptor function and sensitivity to alterations in these parameters. Four receptor systems, EGFR, TfR, LDLR, and VtgR, were studied in hopes of gaining insight as to the classification of and distinctions between one another based solely on responsiveness to changes in avidity, changes in consumption, or a dual responsiveness to variation in both parameters. Analysis confirms the possibility of distinguishing receptor systems from one another based on responsiveness to changes in avidity and consumption. Additional, extension of the analysis to incorporate two additional signaling receptors, GCSFR and EpoR, has identified the beginnings of trend that shows signaling receptors to display dual sensitivity while transport receptors show a sensitivity to either avidity or consumption. Of the systems studied, TfR and LDLR were found to be responsive to avidity, VtgR was responsive to consumption, and EGFR, GCSFR, and EpoR displayed sensitivity to both parameters. As another extension, we expanded analysis by incorporating another method of varying avidity. Results from this analysis matched up nearly perfectly with those earlier in the paper, confirming the sensitivities to avidity and not some other parameter that simply happened to change with avidity in the original analysis. Furthermore, it was discovered that increases in either parameter led to higher levels of receptor system function and efficiency. Then, by inspecting the mathematical definitions for these parameters, further insight was gained about optimal strategy to improve receptor system efficiency. Specific avidity is defined as $\gamma = \frac{Q_R*K_{on}*10^9}{k_t*k_{off}*V*N_{av}}$. In the paper, avidity was increased by simply decreasing the extracellular volume, V, but we feel it is unlikely that in practice cells will be able to exert control over extracellular volume. On the other hand, a more reasonable approach for the cell to increasing its avidity would be to increase its balance between the rate of introduction of new surface receptors, Q_R, and the rate of internalization of receptors that haven't yet bonded with a ligand, k_t. This ratio is defined as R_T, the steady state receptor abundance, and other studies have found the LDLR and TfR systems to have receptor levels that can readily be controlled by growth factor signals. The VtgR system was found to be consumption sensitive which was quantified through variation of the partition coefficient, $\beta = k_e/k_{off}$, by increasing k_e. Findings from Opresko et. al. support this method of variation, finding that hormones were able to significantly impact the internalization rate of vitellogenin in frog oocytes. Now, we examine the EGFR system which displayed a dual sensitivity to both avidity and consumption, and will consequently experience the greatest gain in efficiency if both avidity and consumption are simultaneously increased. The original authors present that the optimal method of causing this increase would be a to both increase the rate of internalization of bonded receptors, k_e, while simultaneously decreasing the the rate of internalization of free receptors, k_t. Research supports this decoupled control over the two internalization rates for the EGFR system. However, we believe that another likely option would be for the receptor system to develop a stronger bonding mechanism with the relatively small signaling ligands which would reduce k_off, one of the factors in the denominator of both $\beta$ and $\gamma$. Perhaps the biological pressure of increased signaling efficiency will eventually lead to this adaption in the EGFR system. At this point, we see that this mathematical modeling study has allowed us to classify and distinguish between the selected receptor systems as well as providing biologically relevant results that agree with other work in the field of receptor mediated endocytosis.

Recreation Discrepancies

My partner and I selected to recreate the differential model for receptor mediated endocytosis as used by Shankaran et. al. to study the sensitivity and categorize four different receptor systems. We have been able to entirely recreate results presented in the original paper and have found no final differences between our figures and those presented in the paper. Our only difficulties in recreating the original work stemmed from two typos we encountered in the Materials and Methods section of the paper. The first was an error in defining a pair of intermediate terms used for the log log scaling in the contour plots. This mistake was mathematically obvious and we were able to notice and correct it easily. Unfortunately, the second typo proved more difficult to recognize. This error led to a 180 degree rotation of the contour plot presented as Figure 5b. After carefully checking our work multiple times, we contacted to the original authors for insight as to what our mistake might be. Thankfully, the authors responded quickly with an answer to our problem. A second typo had been made, again with a pair of intermediate terms used in the creation of the contour plots. After making this simple correction in our Mathematica code, we attained a correct plot for Figure 5b and were subsequently able to create the robustness contour plot. Once again, we are happy with our completely accurate recreation and the timely help provided to us by the original authors.

Limitations

While the receptor mediated endocytosis model studied by Shankaran et. al. has clearly produced meaningful results that are supported by other biological evidence, it is subject to simplifications that may somewhat limit its application. One such limitation is the model's lack of account for receptor recycling after internalization, a process that is know to occur in both LDLR and TfR. Instead, the model merely pushes all recycling effects into a single constant creation rate of new free receptors on the cell surface, an approximation that has not be investigated for appropriateness. Stemming from the disregard for receptor recycling, the model also ignores all processes post internalization such as receptor ligand dissociation, receptor and ligand degradation, and in some cases ligand recycling. As with chemical reactions, its possible that one of these steps could be rate-limiting for the entire endocytosis process and consequently lead to qualitatively different results. Another possible rate limiting process that was not addressed was the formation of clatherin coated pits in the cell wall that eventually close off to complete the internalization of the receptor-ligand complexes. Possible limitation could stem from either a shortage of clatherin or an over-reduction of the cell wall available for internalization. Lund et. al. displays support for such possible rate limitations with their Satin plots which show a saturation of internalization rate with increased ligand concentrations. However, while these possible limitations are clearly relevant to the modeling or receptor mediated endocytosis, they do not necessarily retract from the significance of the results obtained by Shankaran et. al. We remember that the original paper emphasizes low ligand concentrations and likely avoids running into such limitation by doing so. Nonetheless, it is likely to be an important development to this field of study when a model is developed that incorporates the total kinetics of the process of receptor mediated endocytosis.

Comparison to Other Models

After recognizing some of the limitations experienced by the model that we had selected to reproduce, we searched for similar results in models that had incorporated some of the disregarded kinetics in our model. In one paper we found, Radhakrishnan et. al. examined the impact of varying the ligand affinity of receptor systems, k_off. Their model includes elements of recycling and both receptor and ligand degradation absent from our model. We remember the mathematical definitions of specific avidity and the partition coefficient from our model possessed k_off terms as factors of their respective denominators. Therefore, as we suggested for the EGFR system, a reduction in k_off, symbolizing an increase in bonding affinity, would increase both parameters and hence would be expected to increase efficiency of the system. Consequently, results from Radhakrishnan et. al. show reduced k_off correlating to higher initial receptor system activity and greater quantities of internalized ligands which is consistent with our expectations. A similar model used by Krippendorff et. al. that also incorporates post internalization processes looked at the impact of increasing the complex internalization rate, k_e. K_e is the numerator of the partition coefficient in our model, so we would expect an increase to result in at least some degree of increase in efficiency depending on how robust the given system is from the start. Once again, true to our expectation, their results showed increased efficiency with increased k_e. Then Krippendorff et. al. looked at the same changes in k_e but reduced both k_on and k_off by factors of 100. Again remembering k_off is a term in the denominators of both avidity and partition coefficient, we would expect this change to make the system much more robust to changes in k_e, which is in fact what their results show. Having identified these similar results in other papers that examined more complete models for receptor mediated endocytosis, we have come to believe that the assumptions made in our model have not had significant impact on its qualitative results.

MATHEMATICA NOTEBOOKS

Figure 2a

Figure 2b

Figure 2c

Figure 3

Figure 4

Figure 5a

Figure 5b

Figure 5 Approx. Test

Figure 6

Equilibrium Extension

Varying kt Extension

New Receptors Extension

REFERENCES