/*
 * Inclusion of the glucocorticoid receptor in a hypothalamic pituitary
 * adrenal axis model reveals bistability
 * 
 * Model Status
 * 
 * This particular version of the model has a pulsatile stress
 * and recreates the shapes of the graphs in figure 5 of the paper.
 * This CellML model runs in both COR and OpenCell. The units have
 * been checked and they are consistent. The equations are based
 * on the scaled set in the paper (9-12) and the concentrations
 * of the hormones are therefore dimensionless. For this reason
 * the CellML model doesn't quite recreate the graphs in figure
 * 5 (where the hormones are assigned units of picomolar and nanomolar).
 * 
 * Model Structure
 * 
 * ABSTRACT: BACKGROUND: The body's primary stress management system
 * is the hypothalamic pituitary adrenal (HPA) axis. The HPA axis
 * responds to physical and mental challenge to maintain homeostasis
 * in part by controlling the body's cortisol level. Dysregulation
 * of the HPA axis is implicated in numerous stress-related diseases.
 * RESULTS: We developed a structured model of the HPA axis that
 * includes the glucocorticoid receptor (GR). This model incorporates
 * nonlinear kinetics of pituitary GR synthesis. The nonlinear
 * effect arises from the fact that GR homodimerizes after cortisol
 * activation and induces its own synthesis in the pituitary. This
 * homodimerization makes possible two stable steady states (low
 * and high) and one unstable state of cortisol production resulting
 * in bistability of the HPA axis. In this model, low GR concentration
 * represents the normal steady state, and high GR concentration
 * represents a dysregulated steady state. A short stress in the
 * normal steady state produces a small perturbation in the GR
 * concentration that quickly returns to normal levels. Long, repeated
 * stress produces persistent and high GR concentration that does
 * not return to baseline forcing the HPA axis to an alternate
 * steady state. One consequence of increased steady state GR is
 * reduced steady state cortisol, which has been observed in some
 * stress related disorders such as Chronic Fatigue Syndrome (CFS).
 * CONCLUSION: Inclusion of pituitary GR expression resulted in
 * a biologically plausible model of HPA axis bistability and hypocortisolism.
 * High GR concentration enhanced cortisol negative feedback on
 * the hypothalamus and forced the HPA axis into an alternative,
 * low cortisol state. This model can be used to explore mechanisms
 * underlying disorders of the HPA axis.
 * 
 * The original paper reference is cited below:
 * 
 * Inclusion of the glucocorticoid receptor in a hypothalamic pituitary
 * adrenal axis model reveals bistability, Shakti Gupta, Eric Aslakson,
 * Brian M. Gurbaxani and Suzanne D. Vernon, 2007,Theoretical Biology
 * and Medical Modelling, volume 4, issue 8. PubMed ID: 17300722
 * 
 * model diagram
 * 
 * [[Image file: gupta_2007.png]]
 * 
 * The HPA axis model has three compartments, namely the hypothalamus,
 * pituitary and adrenals. F is an external stress that triggers
 * CRH (C) secretion by the hypothalamus. In turn, this signals
 * to the pituitary to release ACTH (A) which stimulates the adrenal
 * gland to release cortisol (O). Cortisol binds to the glucocorticoid
 * receptor (R) in the pituitary and has a negative feedback effect
 * on the release of CRH and ACTH. In addition, the cortisol-receptor
 * complex has a positive feedback effect, promoting further receptor
 * synthesis.
 */

import nsrunit;
unit conversion on;
unit hour=3600 second^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 c(time) dimensionless;
	when(time=time.min) c=0.6;
	real f(time) dimensionless;
	real kcd first_order_rate_constant;
	kcd=1.0;
	real ki1 dimensionless;
	ki1=0.1;
	real o(time) dimensionless;
	when(time=time.min) o=0.055;
	real a(time) dimensionless;
	when(time=time.min) a=0.055;
	real kad first_order_rate_constant;
	kad=10.0;
	real ki2 first_order_rate_constant;
	ki2=0.1;
	real r(time) dimensionless;
	when(time=time.min) r=0.01;
	real kcr first_order_rate_constant;
	kcr=0.05;
	real krd first_order_rate_constant;
	krd=0.9;
	real k dimensionless;
	k=0.001;

	// <component name="environment">

	// <component name="c">
	c:time=((1 first_order_rate_constant)*((1+f)/(1+o/ki1))-kcd*c);
	f=(if ((((((time>(0 hour)) and (time<(1 hour))) or ((time>(4 hour)) and (time<(5 hour)))) or ((time>(8 hour)) and (time<(9 hour)))) or ((time>(12 hour)) and (time<(13 hour)))) or ((time>(16 hour)) and (time<(17 hour)))) 1 else 0);

	// <component name="a">
	a:time=(c/((1 hour)+o*r/ki2)-kad*a);

	// <component name="r">
	r:time=((o*r)^2/((1 hour)*(k+(o*r)^2))+kcr-krd*r);

	// <component name="o">
	o:time=((1 first_order_rate_constant)*(a-o));

	// <component name="reaction_constants">
}