/*
 * A Model For Circadian PER Oscillations In Drosophila
 * 
 * Model Status
 * 
 * This CellML model runs in both OpenCell and COR to reproduce
 * Figure 2 in the original published paper. The units have been
 * checked and they are consistent.
 * 
 * Model Structure
 * 
 * ABSTRACT: The mechanism of circadian oscillations in the period
 * protein (PER) in Drosophila is investigated by means of a theoretical
 * model. Taking into account recent experimental observations,
 * the model for the circadian clock is based on multiple phosphorylation
 * of PER and on the negative feedback exerted by PER on the transcription
 * of the period (per) gene. This minimal biochemical model provides
 * a molecular basis for circadian oscillations of the limit cycle
 * type. During oscillations, the peak in per mRNA precedes by
 * several hours the peak in total PER protein. The results support
 * the view that multiple PER phosphorylation introduces times
 * delays which strengthen the capability of negative feedback
 * to produce oscillations. The analysis shows that the rhythm
 * only occurs in a range bounded by two critical values of the
 * maximum rate of PER degradation. A similar result is obtained
 * with respect to the rate of PER transport into the nucleus.
 * The results suggest a tentative explanation for the altered
 * period of per mutants, in terms of variations in the rate of
 * PER degradation.
 * 
 * The original paper reference is cited below:
 * 
 * A Model for Circadian Oscillations in the Drosophila Period
 * Protein (PER), Albert Goldbeter, 1995, Proceedings of the Royal
 * Society of London, Series B, Biological Sciences, 261, 319-324.
 * PubMed ID: 8587874
 * 
 * reaction diagram
 * 
 * [[Image file: goldbeter_1995.png]]
 * 
 * Schematic diagram of the model for circadian oscillations in
 * PER protein and per mRNA.
 */

import nsrunit;
unit conversion on;
// unit micromolar predefined
unit hour=3600 second^1;
unit flux=2.7777778E-7 meter^(-3)*second^(-1)*mole^1;
unit first_order_rate_constant=2.7777778E-4 second^(-1);

math main {
	realDomain time hour;
	time.min=0;
	extern time.max;
	extern time.delta;
	real M(time) micromolar;
	when(time=time.min) M=0.6;
	real Km micromolar;
	Km=0.5;
	real vs flux;
	vs=0.76;
	real vm flux;
	vm=0.65;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real n dimensionless;
//	n=4.0;
	real KI micromolar;
	KI=1.0;
	real PN(time) micromolar;
	when(time=time.min) PN=1.1;
	real P0(time) micromolar;
	when(time=time.min) P0=0.5;
	real ks first_order_rate_constant;
	ks=0.38;
	real P1(time) micromolar;
	when(time=time.min) P1=0.6;
	real K1 micromolar;
	K1=2.0;
	real V1 flux;
	V1=3.2;
	real K2 micromolar;
	K2=2.0;
	real V2 flux;
	V2=1.58;
	real P2(time) micromolar;
	when(time=time.min) P2=0.6;
	real K3 micromolar;
	K3=2.0;
	real V3 flux;
	V3=5.0;
	real K4 micromolar;
	K4=2.0;
	real V4 flux;
	V4=2.5;
	real Kd micromolar;
	Kd=0.2;
	real vd flux;
	vd=0.95;
	real k1 first_order_rate_constant;
	k1=1.9;
	real k2 first_order_rate_constant;
	k2=1.3;
	real Pt(time) micromolar;

	// <component name="environment">

	// <component name="M">
	M:time=(vs*(KI^4/(KI^4+PN^4))-vm*(M/(Km+M)));

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

	// <component name="P1">
	P1:time=(V1*(P0/(K1+P0))-(V2*(P1/(K2+P1))+V3*(P1/(K3+P1)))+V4*(P2/(K4+P2)));

	// <component name="P2">
	P2:time=(V3*(P1/(K3+P1))-(V4*(P2/(K4+P2))+k1*P2)+k2*PN-vd*(P2/(Kd+P2)));

	// <component name="PN">
	PN:time=(k1*P2-k2*PN);

	// <component name="Pt">
	Pt=(P0+P1+P2+PN);

	// <component name="parameters">
}