/*
 * Regulation of the G1-S Transition of the Mammalian Cell Cycle
 * 
 * Model Status
 * 
 * This CellML model has been recoded to remove the reaction element.
 * The model opens in PCEnv and COR but unfortunately it does not
 * run due to there being 'circular arguments' (COR) or the model
 * being 'under-constrained' (PCEnv). At the moment this is a limitation
 * of the simulation software because it is unable to process nonlinear
 * algebraic equations that must be solved iteratively at each
 * time step.
 * 
 * Model Structure
 * 
 * In order to better understand developmental and tumour biology,
 * extensive research on the control of the mammalian cell cycle
 * is being carried out. At the protein level, several molecules
 * play a role in the control of the G1-S transition of the cell
 * cycle. Specifically, G1 cyclins such as cyclin E are known to
 * form phosphorylated complexes with cyclin dependent kinases
 * (cdk's, such as cdk2). These complexes, which can be inhibited
 * by other proteins, can phosphorylate pocket proteins such as
 * pRb. Hyperphosphorylation of pRb occurs at the time of E2F-1
 * release. This transcription factor initiates the transcription
 * of genes involved in S phase activities, including DNA polymerase.
 * 
 * Hatzimanikatis, Lee and Bailey capture the components and interactions
 * of cyclin E, cdk2, pRb, inhibitor and E2F-1 during the G1-S
 * transition of the cell cycle in their mathematical model (see
 * the figure below).
 * 
 * The complete original paper reference is cited below:
 * 
 * A Mathematical Description of Regulation of the G1-S Transition
 * of the Mammalian Cell Cycle, V. Hatzimanikatis, K. H. Lee and
 * J. E. Bailey, 1999, Biotechnology and Bioengineering, 65, 631-637.
 * PubMed ID: 10550769
 * 
 * diagram of the cell cycle
 * 
 * [[Image file: hatzimanikatis_1999.png]]
 * 
 * A schematic diagram of the reactions involved during the G1-S
 * transition of the mammalian cell cycle.
 */

import nsrunit;
unit conversion on;

math main {
	realDomain time dimensionless;
	time.min=0;
	extern time.max;
	extern time.delta;
	real C(time) dimensionless;
	when(time=time.min) C=0.38;
	real gamma dimensionless;
	gamma=1.0;
	real Vs(time) dimensionless;
	real V1(time) dimensionless;
	real V2(time) dimensionless;
	real Vd(time) dimensionless;
	real K(time) dimensionless;
	when(time=time.min) K=0.1;
	real RP(time) dimensionless;
	when(time=time.min) RP=1.0;
	real V3(time) dimensionless;
	real V4(time) dimensionless;
	real E(time) dimensionless;
	real sigma dimensionless;
	sigma=10;
	real RE(time) dimensionless;
	real KP(time) dimensionless;
	real KPI(time) dimensionless;
	real thetaI dimensionless;
	thetaI=1.0;
	real I(time) dimensionless;
	real lambda dimensionless;
	lambda=1.0;
	real thetaE dimensionless;
	thetaE=0.01;
	real R(time) dimensionless;
	real VCs dimensionless;
	VCs=1.0;
	real Vsm dimensionless;
	Vsm=1.0;
	real KsE dimensionless;
	KsE=0.1;
	real V1m dimensionless;
	V1m=50.0;
	real K1C dimensionless;
	K1C=0.1;
	real K1 dimensionless;
	K1=0.0001;
	real V2m dimensionless;
	V2m=40;
	real K2 dimensionless;
	K2=0.0001;
	real V3m dimensionless;
	V3m=3000;
	real K3 dimensionless;
	K3=0.0001;
	real V4m dimensionless;
	V4m=3.0;
	real K4 dimensionless;
	K4=0.0001;
	real VdEm dimensionless;
	VdEm=1000.0;
	real KdC dimensionless;
	KdC=0.005;

	// <component name="environment">

	// <component name="C">
	C:time=(Vs+gamma*V2-(gamma*V1+Vd));

	// <component name="K">
	K:time=(V2-V1);

	// <component name="RP">
	RP:time=(V3-V4);

	// <component name="E">
	E=(1-sigma*RE);

	// <component name="KP">
	KP=(1-(KPI+K));

	// <component name="KPI">
	KPI=(thetaI*KP*I);

	// <component name="I">
	I=(1-lambda*KPI);

	// <component name="RE">
	RE=(thetaE*R*E);

	// <component name="R">
	R=(1-(RP+RE));

	// <component name="Vs">
	Vs=(VCs+Vsm*(E/(KsE+E)));

	// <component name="V1">
	V1=(V1m*(C/(K1C+C))*(K/(K1+K)));

	// <component name="V2">
	V2=(V2m*(KP/(K2+KP)));

	// <component name="V3">
	V3=(V3m*KP*(RE/(K3+RE)));

	// <component name="V4">
	V4=(V4m*(RP/(K4+RP)));

	// <component name="Vd">
	Vd=(C+VdEm*E*(C/(KdC+C)));
}