Exemplary Model Description Draft 5

State Variables and Parameters

State Variables

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

Table 2: Independent Parameters

Variable	Description	EGFR Value	TfR Value	LDLR Value	VtgR Value
kon	Rate of formation of receptor-ligand complexes (/(min * nM))	0.097	0.0030	0.0028	5.4 E -5
koff	Rate of dissociation of or receptor-ligand complexes (/min)	0.24	0.09	0.04	0.07
ke	Rate of internalization of receptor-ligand complexes (/min)	0.15	0.6	0.195	0.108
kt	Rate of internalization of unbonded receptors (/min)	0.02	0.6	0.195	0.108
QR	Rate of creation of new receptors on the cell surface (molecules/min)	4000	15600	390	2.2 E10
V	Extracellular volume (liters)	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

Table 3: Calculated Parameters

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

• • N_av = Avogadro's number

Model Assumptions

• 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 [ref class text?]. 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.

$f({\mathbf{x}}) \approx f({\mathbf{p}}) + \left. {\nabla f} \right|_{\mathbf{p}} \cdot ({\mathbf{x}} - {\mathbf{p}})$[not sure where to find a reference for this]

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 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 in 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 the characteristics of the standard pulse used in B. 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 an 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 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, 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. The original paper defines:

$r_s = \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.

Reference