/*
 * Controlling Cell Cycle Dynamics Using a Reversibly Binding Inhibitor
 * 
 * Model Status
 * 
 * This CellML model runs in both OpenCell and COR to reproduce
 * the published results. The units have been checked and they
 * are consistent.
 * 
 * Model Structure
 * 
 * ABSTRACT: We demonstrate, by using mathematical modeling of
 * cell division cycle (CDC) dynamics, a potential mechanism for
 * precisely controlling the frequency of cell division and regulating
 * the size of a dividing cell. Control of the cell cycle is achieved
 * by artificially expressing a protein that reversibly binds and
 * inactivates any one of the CDC proteins. In the simplest case,
 * such as the checkpoint-free situation encountered in early amphibian
 * embryos, the frequency of CDC oscillations can be increased
 * or decreased by regulating the rate of synthesis, the binding
 * rate, or the equilibrium constant of the binding protein. In
 * a more complex model of cell division, where size-control checkpoints
 * are included, we show that the same reversible binding reaction
 * can alter the mean cell mass in a continuously dividing cell.
 * Because this control scheme is general and requires only the
 * expression of a single protein, it provides a practical means
 * for tuning the characteristics of the cell cycle in vivo.
 * 
 * The original paper reference is cited below:
 * 
 * A theory for controlling cell cycle dynamics using a reversibly
 * binding inhibitor, Timothy S. Gardner, Milos Dolnik and James
 * J. Collins, 1998, Proceedings of the National Academy of Sciences,
 * 95, 14190-14195. PubMed ID: 9826676
 * 
 * cell diagram
 * 
 * [[Image file: gardner_1998.png]]
 * 
 * Control of the Goldbeter 1991 model with a cyclin inhibitor.
 * M represents cdc2 kinase, X represents the fraction of active
 * (phosphorylated) cyclin protease, and * represents the fraction
 * of inactive enzymes. The Goldbeter model is outlined by the
 * dashed box. Solid arrows indicate protein synthesis, degradation
 * or enzymatic conversion. Dashed arrows represent activation.
 */

import nsrunit;
unit conversion on;
unit minute=60 second^1;
unit per_minute=.01666667 second^(-1);
// unit micromolar predefined

math main {
	realDomain time minute;
	time.min=0;
	extern time.max;
	extern time.delta;
	real C(time) dimensionless;
	when(time=time.min) C=0.01;
	real Ctot(time) dimensionless;
	real vi per_minute;
	vi=0.1;
	real k1 per_minute;
	k1=0.5;
	real X(time) dimensionless;
	when(time=time.min) X=0.01;
	real K5 dimensionless;
	K5=0.02;
	real kd per_minute;
	kd=0.02;
	real Z(time) dimensionless;
	when(time=time.min) Z=0.28;
	real M(time) dimensionless;
	when(time=time.min) M=0.01;
	real M_ dimensionless;
	real V1(time) per_minute;
	real K1 dimensionless;
	K1=0.02;
	real V2 per_minute;
	V2=0.25;
	real K2 dimensionless;
	K2=0.02;
	real X_(time) dimensionless;
	real V3(time) per_minute;
	real K3 dimensionless;
	K3=0.1;
	real V4 per_minute;
	V4=0.1;
	real K4 dimensionless;
	K4=0.1;
	real Y(time) dimensionless;
	when(time=time.min) Y=0.01;
	real Ytot(time) dimensionless;
	real vs per_minute;
	vs=0.1;
	real d1 per_minute;
	d1=0.05;
	real a1 per_minute;
	a1=1.5;
	real a2 per_minute;
	a2=1.5;
	real alpha dimensionless;
	alpha=0.1;
	real BP(time) dimensionless;
	real Kd dimensionless;
	Kd=1;
	real K6 dimensionless;
	K6=0.3;
	real V1_dash per_minute;
	V1_dash=0.75;
	real V3_dash per_minute;
	V3_dash=0.3;

	// <component name="environment">

	// <component name="C">
	C:time=(vi-k1*X*C/(C+K5)-kd*C);
	Ctot=(C+Z);

	// <component name="M">
	M_=((-1)*1);
	M:time=(V1*M_/(M_+K1)-V2*M/(M+K2));

	// <component name="X">
	X_=(1-X);
	X:time=(V3*X_/(X_+K3)-V4*X/(X+K4));

	// <component name="Y">
	Y:time=(vs-d1*Y-a1*C*Y+(a2+alpha*kd)*Z);
	Ytot=(Y+Z);

	// <component name="Z">
	Z:time=(a1*C*Y-(a2+alpha*kd+alpha*d1)*Z);

	// <component name="BP">
	BP=(1+Kd*Ytot/(C+Kd)^2);

	// <component name="model_parameters">
	V1=(C*V1_dash/(C+K6));
	V3=(M*V3_dash);
}