/*
 * Computational model for effects of ligand/receptor binding properties
 * on interleukin-2 trafficking dynamics and T cell proliferation
 * response
 * 
 * Model Status
 * 
 * The model is known to run in OpenCell and COR and the results
 * replicate the published paper (Figure 5), as generated by the
 * original model. The units have been checked and they are consistent.
 * 
 * Model Structure
 * 
 * ABSTRACT: Multisubunit cytokine receptors such as the heterotrimeric
 * receptor for interleukin-2 (IL-2) are ubiquitous in hematopoeitic
 * cell types of importance in biotechnology and are crucial regulators
 * of cell proliferation and differentiation behavior. Dynamics
 * of cytokine/receptor endocytic trafficking can significantly
 * impact cell responses through effects of receptor down-regulation
 * and ligand depletion, and in turn are governed by ligand/receptor
 * binding properties. We describe here a computational model for
 * trafficking dynamics of the IL-2 receptor (IL-2R) system, which
 * is able to predict T cell proliferation responses to IL-2. This
 * model comprises kinetic equations describing binding, internalization,
 * and postendocytic sorting of IL-2 and IL-2R, including an experimentally
 * derived dependence of cell proliferation rate on these properties.
 * Computational results from this model predict that IL-2 depletion
 * can be reduced by decreasing its binding affinity for the IL-2R
 * betagamma subunit relative to the alpha subunit at endosomal
 * pH, as a result of enhanced ligand sorting to recycling vis-a-vis
 * degradation, and that an IL-2 analogue with such altered binding
 * properties should exhibit increased potency for stimulating
 * the T cell proliferation response. These results are in agreement
 * with our recent experimental findings for the IL-2 analogue
 * termed 2D1 [Fallon, E. M. et al. J. Biol. Chem. 2000, 275, 6790-6797].
 * Thus, this type of model may enable prediction of beneficial
 * cytokine/receptor binding properties to aid development of molecular
 * design criteria for improvements in applications such as in
 * vivo cytokine therapies and in vitro hematopoietic cell bioreactors.
 * 
 * The original paper reference is cited below:
 * 
 * Computational Model for Effects of Ligand/Receptor Binding Properties
 * on Interleukin-2 Trafficking Dynamics and T Cell Proliferation
 * Response, Eric M. Fallon and Douglas A. Lauffenburger, 2000,
 * Biotechnology Progress, 16, 905-916. PubMed ID: 11027188
 * 
 * diagram of the ligand-receptor binding model
 * 
 * [[Image file: fallon_2000.png]]
 * 
 * A schematic diagram of Fallon and Lauffenburger's 2000 computational
 * model of IL-2-receptor binding and trafficking. The model follows
 * the path of an extracellular ligand as it is bound to a cell
 * surface receptor, internalised, degraded or recycled.
 */

import nsrunit;
unit conversion on;
// unit picomolar predefined
unit minute=60 second^1;
unit flux=1.6666667E-11 meter^(-3)*second^(-1)*mole^1;
unit first_order_rate_constant=.01666667 second^(-1);
unit second_order_rate_constant=1.6666667E7 meter^3*second^(-1)*mole^(-1);
unit cell = fundamental;
unit number = fundamental;
unit number_per_cell=1 cell^(-1)*number^1;
unit number_per_picomole=1E12 mole^(-1)*number^1;
unit number_per_cell_minute=.01666667 second^(-1)*cell^(-1)*number^1;
unit L_per_cell=.001 meter^3*cell^(-1);
unit cell_per_L_minute=16.66666667 meter^(-3)*second^(-1)*cell^1;
unit cell_per_L=1E3 meter^(-3)*cell^1;

math main {
	realDomain time minute;
	time.min=0;
	extern time.max;
	extern time.delta;
	real Rs(time) number_per_cell;
	when(time=time.min) Rs=1500;
	real L(time) picomolar;
	when(time=time.min) L=10;
	real Cs(time) number_per_cell;
	when(time=time.min) Cs=1;
	real Vs number_per_cell_minute;
	Vs=11;
	real kf second_order_rate_constant;
	real kr first_order_rate_constant;
	kr=0.0138;
	real kt first_order_rate_constant;
	kt=0.007;
	real ksyn first_order_rate_constant;
	ksyn=0.0011;
	real ke first_order_rate_constant;
	ke=0.04;
	real Ri(time) number_per_cell;
	when(time=time.min) Ri=300;
	real Li(time) picomolar;
	when(time=time.min) Li=0.01;
	real Ci(time) number_per_cell;
	when(time=time.min) Ci=1;
	real kfe second_order_rate_constant;
	real kre first_order_rate_constant;
	real kh first_order_rate_constant;
	kh=0.035;
	real kx first_order_rate_constant;
	kx=0.15;
	real Ve L_per_cell;
	Ve=1e-14;
	real NA number_per_picomole;
	NA=6.022e11;
	real Ld(time) number_per_cell;
	when(time=time.min) Ld=1;
	real Y(time) cell_per_L;
	when(time=time.min) Y=2.5e8;
	real IL2 dimensionless;
	IL2=1;

	// <component name="environment">

	// <component name="Rs">
	Rs:time=((kr+ksyn)*Cs+Vs-(kf*L*Rs+kt*Rs));

	// <component name="Cs">
	Cs:time=(kf*L*Rs-(kr+ke)*Cs);

	// <component name="Ri">
	Ri:time=(kre*Ci+kt*Rs-(kfe*Li*Ri+kh*Ri));

	// <component name="Ci">
	Ci:time=(kfe*Li*Ri+ke*Cs-(kre+kh)*Ci);

	// <component name="Li">
	Li:time=((kre*Ci-kfe*Li*Ri)/(Ve*NA)-kx*Li);

	// <component name="Ld">
	Ld:time=(kh*Ci);

	// <component name="L">
	L:time=((kr*Cs+kx*Li*Ve*NA-kf*L*Rs)*Y/NA);

	// <component name="Y">
	Y:time=(if ((600*Cs/((250 number_per_cell)+Cs)-200)>0) (600*Cs/((250 number_per_cell)+Cs)-200)*(1E3 cell_per_L_minute) else (0 cell_per_L_minute));

	// <component name="model_parameters">
	kf=(if (IL2=1) kr/(11.1 picomolar) else kr/(8.2 picomolar));
	kre=(if (IL2=1) kr*8 else kr*5);
	kfe=(if (IL2=1) kre/(1E3 picomolar) else kre/(3E3 picomolar));
}