/*
 * Sharp Developmental Thresholds Defined Through Bistability By
 * Antagonistic Gradients of Retinoic Acid and FGF Signaling
 * 
 * Model Status
 * 
 * This CellML model runs in OpenCell and COR. It was created from
 * equations [1] to [5], and [12] and [13]. The model parameters
 * were taken from the caption of Figure 3. The CellML model cannot
 * recreate the figures from the paper due to constant values for
 * the parameters K_I and K_A. The units have been checked and
 * they are consistent.
 * 
 * Model Structure
 * 
 * ABSTRACT: The establishment of thresholds along morphogen gradients
 * in the embryo is poorly understood. Using mathematical modeling,
 * we show that mutually inhibitory gradients can generate and
 * position sharp morphogen thresholds in the embryonic space.
 * Taking vertebrate segmentation as a paradigm, we demonstrate
 * that the antagonistic gradients of retinoic acid (RA) and Fibroblast
 * Growth Factor (FGF) along the presomitic mesoderm (PSM) may
 * lead to the coexistence of two stable steady states. Here, we
 * propose that this bistability is associated with abrupt switches
 * in the levels of FGF and RA signaling, which permit the synchronized
 * activation of segmentation genes, such as mesp2, in successive
 * cohorts of PSM cells in response to the segmentation clock,
 * thereby defining the future segments. Bistability resulting
 * from mutual inhibition of RA and FGF provides a molecular mechanism
 * for the all-or-none transitions assumed in the "clock and wavefront"
 * somitogenesis model. Given that mutually antagonistic signaling
 * gradients are common in development, such bistable switches
 * could represent an important principle underlying embryonic
 * patterning.
 * 
 * The original paper reference is cited below:
 * 
 * Sharp developmental thresholds defined through bistability by
 * antagonistic gradients of retinoic acid and FGF signaling. Goldbeter
 * A, Gonze D, Pourquie O. 2007, Developmental Dynamics, 236, 1495-1508.
 * PubMed ID: 17497689
 * 
 * reaction diagram
 * 
 * [[Image file: goldbeter_2007.png]]
 * 
 * Scheme of the regulatory interactions between RA and FGF signaling.
 * Synthesis of RA is catalyzed by the enzyme RALDH2, while RA
 * is hydrolyzed by the enzyme CYP26. The inhibitory effect exerted
 * on RA by FGF arises from the induction of cyp26 expression by
 * FGF. The inhibition exerted by RA on FGF occurs through impeding
 * the rate of fgf8 mRNA translation. As shown by the model built
 * according to this regulatory scheme, bistability readily arises
 * from the mutual inhibition between RA and FGF.
 */

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 RA(time) nanomolar;
	when(time=time.min) RA=0.1;
	real v_s1 flux;
	real k_d1 second_order_rate_constant;
	k_d1=1;
	real RA.C nanomolar;
	RA.C=0.1;
	real k_d5 first_order_rate_constant;
	k_d5=0;
	real M_C(time) nanomolar;
	when(time=time.min) M_C=0.1;
	real V_0 flux;
	V_0=0.365;
	real V_sC flux;
	V_sC=7.1;
	real F(time) nanomolar;
	when(time=time.min) F=0.0001;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real n dimensionless;
//	n=2;
	real K_A nanomolar;
	K_A=0.2;
	real k_d3 first_order_rate_constant;
	k_d3=1;
	real C.C(time) nanomolar;
	when(time=time.min) C.C=0.1;
	real k_s2 first_order_rate_constant;
	k_s2=1;
	real k_d2 first_order_rate_constant;
	k_d2=0.28;
	real k_s3 first_order_rate_constant;
	k_s3=1;
	real M_F nanomolar;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real m dimensionless;
//	m=2;
	real K_I nanomolar;
	K_I=0.2;
	real k_d4 first_order_rate_constant;
	k_d4=1;
	real k_s1 first_order_rate_constant;
	k_s1=1;
	real RALDH2_0 nanomolar;
	RALDH2_0=7.1;
	real x dimensionless;
	x=15;
	real L dimensionless;
	L=50;
	real M_0 nanomolar;
	M_0=5;
	real alpha_1(time) dimensionless;
	real K_r1 nanomolar;
	K_r1=1;
	real alpha_2(time) dimensionless;
	real K_r2 nanomolar;
	K_r2=1;
	real rho(time) dimensionless;

	// <component name="environment">

	// <component name="RA">
	RA:time=(v_s1-k_d1*RA.C*RA-k_d5*RA);

	// <component name="M_C">
	M_C:time=(V_0+V_sC*F^2/(K_A^2+F^2)-k_d3*M_C);

	// <component name="C">
	C.C:time=(k_s2*M_C-k_d2*C.C);

	// <component name="F">
	F:time=(k_s3*M_F*K_I^2/(K_I^2+RA^2)-k_d4*F);

	// <component name="v_s1">
	v_s1=(k_s1*RALDH2_0*(1-x/L));

	// <component name="M_F">
	M_F=(M_0*x/L);

	// <component name="alpha_1">
	alpha_1=(RA/(RA+K_r1));

	// <component name="alpha_2">
	alpha_2=(F/(F+K_r2));

	// <component name="rho">
	rho=(alpha_2/alpha_1);

	// <component name="model_parameters">
}