/*
 * A synthetic oscillatory network of transcriptional regulators.
 * 
 * Model Status
 * 
 * This CellML version of this model is not able to reproduce the
 * results shown in Fig1c of the original publication, as the initial
 * conditions for protein concentrations were not provided. The
 * original published equations were scaled and modified with reference
 * to the same model on the Biomodels database (BIOMD0000000012
 * - Elowitz2000_Repressilator). Once the model comes to equilibrium
 * (t > 400 minutes,) its output is correct.
 * 
 * The units have been checked in this model and are consistent.
 * 
 * Model Structure
 * 
 * ABSTRACT: Networks of interacting biomolecules carry out many
 * essential functions in living cells, but the 'design principles'
 * underlying the functioning of such intracellular networks remain
 * poorly understood, despite intensive efforts including quantitative
 * analysis of relatively simple systems. Here we present a complementary
 * approach to this problem: the design and construction of a synthetic
 * network to implement a particular function.
 * 
 * We used three transcriptional repressor systems that are not
 * part of any natural biological clock to build an oscillating
 * network, termed the repressilator, in Escherichia coli. The
 * network periodically induces the synthesis of green fluorescent
 * protein as a readout of its state in individual cells. The resulting
 * oscillations, with typical periods of hours, are slower than
 * the cell-division cycle, so the state of the oscillator has
 * to be transmitted from generation to generation. This artificial
 * clock displays noisy behaviour, possibly because of stochastic
 * fluctuations of its components. Such 'rational network design'
 * may lead both to the engineering of new cellular behaviours
 * and to an improved understanding of naturally occurring networks.
 * 
 * The complete original paper reference is cited below:
 * 
 * A synthetic oscillatory network of transcriptional regulators,
 * Michael B. Elowitz and Stanislas Leibler, 2000,Nature: International
 * Weekly Journal of Science, 403, 335-338. PubMed ID: 10659856
 * 
 * Figure 1a
 * 
 * [[Image file: elowitz_2000.png]]
 * 
 * The repressilator network.
 */

import nsrunit;
unit conversion on;
unit minute=60 second^1;
unit first_order_rate_constant=.01666667 second^(-1);
unit per_sec=1 second^(-1);
unit second_per_minute=.01666667  dimensionless;

math main {
	realDomain time minute;
	time.min=0;
	extern time.max;
	extern time.delta;
	real n dimensionless;
	n=2;
	real alpha_0 dimensionless;
	real alpha dimensionless;
	real beta dimensionless;
	real K_m dimensionless;
	K_m=40;
	real efficiency dimensionless;
	efficiency=20;
	real mRNA_halflife minute;
	mRNA_halflife=2;
	real t_ave minute;
	real kd_prot first_order_rate_constant;
	real kd_mRNA first_order_rate_constant;
	real prot_halflife minute;
	prot_halflife=10;
	real tps_repr per_sec;
	tps_repr=5e-4;
	real tps_active per_sec;
	tps_active=0.5;
	real k_tl first_order_rate_constant;
	real a_tr first_order_rate_constant;
	real a0_tr first_order_rate_constant;
	real M_lacl(time) dimensionless;
	when(time=time.min) M_lacl=0;
	real P_cl(time) dimensionless;
	when(time=time.min) P_cl=0;
	real M_tetR(time) dimensionless;
	when(time=time.min) M_tetR=20;
	real P_lacl(time) dimensionless;
	when(time=time.min) P_lacl=0;
	real M_cl(time) dimensionless;
	when(time=time.min) M_cl=0;
	real P_tetR(time) dimensionless;
	when(time=time.min) P_tetR=0;

	// <component name="environment">

	// <component name="parameters">
	t_ave=(mRNA_halflife/ln(2));
	kd_prot=(ln(2)/prot_halflife);
	kd_mRNA=(ln(2)/mRNA_halflife);
	k_tl=(efficiency*kd_mRNA);
	a_tr=((tps_active-tps_repr)*(60 second_per_minute));
	a0_tr=(tps_repr*(60 second_per_minute));
	alpha=(a_tr*efficiency/(kd_prot*K_m));
	alpha_0=(a0_tr*efficiency/(kd_prot*K_m));
	beta=(kd_prot/kd_mRNA);

	// <component name="M_lacl">
	M_lacl:time=(a0_tr+a_tr*K_m^n/(K_m^n+P_cl^n)-kd_mRNA*M_lacl);

	// <component name="M_tetR">
	M_tetR:time=(a_tr*K_m^n/(K_m^n+P_lacl^n)+a0_tr-kd_mRNA*M_tetR);

	// <component name="M_cl">
	M_cl:time=(a_tr*K_m^n/(K_m^n+P_tetR^n)+a0_tr-kd_mRNA*M_cl);

	// <component name="P_lacl">
	P_lacl:time=(k_tl*M_lacl-kd_prot*P_lacl);

	// <component name="P_tetR">
	P_tetR:time=(k_tl*M_tetR-kd_prot*P_tetR);

	// <component name="P_cl">
	P_cl:time=(k_tl*M_cl-kd_prot*P_cl);
}