/*
 * Mathematical model of human growth hormone (hGH)-stimulated
 * cell proliferation explains the efficacy of hGH variants as
 * receptor agonists or antagonists
 * 
 * Model Status
 * 
 * This CellML model runs in both COR and OpenCell. The units have
 * been checked and they are consistent. The CellML model may recreate
 * the results of the original published model but there is no
 * simple validation method as there are no "concentration against
 * time" figures in the paper. The CellML model is based on equations
 * A1a, A1b, A1c, A1d and A1e from the Appendix. Parameter values
 * were taken from table 1 in the paper and were also supplied
 * through correspondence with the original model author.
 * 
 * Model Structure
 * 
 * ABSTRACT: Human growth hormone (hGH) is a therapeutically important
 * endocrine factor that signals various cell types. Structurally
 * and functionally, the interactions of hGH with its receptor
 * have been resolved in fine detail, such that hGH and hGH receptor
 * variants can be practically engineered by either random or rational
 * approaches to achieve significant changes in the free energies
 * of binding. A somewhat unique feature of hGH action is its homodimerization
 * of two hGH receptors, which is required for intracellular signaling
 * and stimulation of cell proliferation, yet the potencies of
 * hGH mutants in cell-based assays rarely correlate with their
 * overall receptor-binding avidities. Here, a mathematical model
 * of hGH-stimulated cell signaling is posed, accounting not only
 * for binding interactions at the cell surface but induction of
 * receptor endocytosis and downregulation as well. Receptor internalization
 * affects ligand potency by imposing a limit on the lifetime of
 * an active receptor complex, irrespective of ligand-receptor
 * binding properties. The model thus explains, in quantitative
 * terms, the numerous published observations regarding hGH receptor
 * agonism and antagonism and challenges the interpretations of
 * previous studies that have not considered receptor trafficking
 * as a central regulatory mechanism in hGH signaling.
 * 
 * model diagram
 * 
 * [[Image file: haugh_2004.png]]
 * 
 * A schematic diagram of the kinetic model of human growth hormone
 * (hGH) receptor binding and trafficking. The extracellular ligand
 * (L), hGH, has two sites for binding the hGH receptor (R), and
 * these are shown in cerise and cyan in the diagram and are numbered
 * 1 and 2 respectively. Site 1 always binds to the receptor first
 * to form a 1:1 ligand-receptor complex (C). A 1:2 dimer (D) may
 * then form through the binding of a second receptor to site 2
 * of the ligand. Dimer dissociation can occur the uncoupling of
 * either hGH site 1 or site 2, but the 1:1 ligand-receptor complexes
 * bound through site 2 dissociate much faster than those bound
 * through site 1. Dimerised complexes (D) are internalised at
 * a higher rate than are free receptors (R) or 1:1 complexes (C),
 * and while internalised dimers are always degraded, internalised
 * 1:1 complexes or free receptors can either be degraded or recycled
 * to the cell surface. A steady state is maintained through de
 * novo receptor synthesis.
 * 
 * The original paper reference is cited below:
 * 
 * Mathematical model of human growth hormone (hGH)-stimulated
 * cell proliferation explains the efficacy of hGH variants as
 * receptor agonists or antagonists, Jason M. Haugh, 2004, Biotechnology
 * Progress, volume 20, issue 5, 1337-1344. PubMed ID: 15458315
 */

import nsrunit;
unit conversion on;
// unit nanomolar predefined
unit minute=60 second^1;
unit flux=1.6666667E-8 meter^(-3)*second^(-1)*mole^1;
unit first_order_rate_constant=.01666667 second^(-1);
unit second_order_rate_constant=1.6666667E4 meter^3*second^(-1)*mole^(-1);

math main {
	realDomain time minute;
	time.min=0;
	extern time.max;
	extern time.delta;
	real C(time) nanomolar;
	when(time=time.min) C=0.0;
	real kf1 second_order_rate_constant;
	kf1=0.1;
	real kr1 first_order_rate_constant;
	real kx2 second_order_rate_constant;
	kx2=4.83;
	real k_x2 first_order_rate_constant;
	k_x2=0.016;
	real kt first_order_rate_constant;
	kt=0.005;
	real D(time) nanomolar;
	when(time=time.min) D=0.0;
	real L nanomolar;
	L=0.01;
	real R(time) nanomolar;
	when(time=time.min) R=2000.0;
	real k_x1 first_order_rate_constant;
	real ke first_order_rate_constant;
	ke=0.10;
	real R_initial nanomolar;
	R_initial=2000.0;
	real krec first_order_rate_constant;
	krec=0.0;
	real Ri(time) nanomolar;
	when(time=time.min) Ri=200.0;
	real Vs flux;
	Vs=10.0;
	real kdeg first_order_rate_constant;
	kdeg=0.05;
	real signal(time) dimensionless;
	real kappaE dimensionless;
	kappaE=0.20;
	real KD nanomolar;
	KD=1.0;

	// <component name="environment">

	// <component name="C">
	C:time=(kf1*L*R+k_x2*D-(kr1+kx2*R+kt)*C);

	// <component name="D">
	D:time=(kx2*R*C-(k_x2+k_x1+ke)*D);

	// <component name="R">
	R:time=(Vs+kr1*C+(k_x2+2*k_x1)*D+krec*Ri-(kf1*L+kx2*C+kt)*R);

	// <component name="Ri">
	Ri:time=(kt*(R+C)-(krec+kdeg)*Ri);

	// <component name="signal">
	signal=(2*D/R_initial/(kappaE+2*D/R_initial));

	// <component name="model_parameters">
	kr1=(KD*kf1);
	k_x1=(.01*kr1);
}