/*
 * Mathematical modeling of the hypothalamic-pituitary-adrenal
 * system activity
 * 
 * Model Status
 * 
 * This CellML model runs in PCEnv and OpenCell, but not in COR
 * because the free variable (time) is dimensionless. The units
 * are consistent throughout, and it runs to reproduce figure 5
 * from the paper. In order to reproduce figure 5, equation 7 was
 * discarded and the variables were defined explicitly as described
 * in the legend for figure 5. An equation was also added describing
 * the SCN periodic input to PVN. The model must be run for roughly
 * 1 day before it stabilizes.
 * 
 * Model Structure
 * 
 * ABSTRACT: Mathematical modeling has proven to be valuable in
 * understanding of the complex biological systems dynamics. In
 * the present report we have developed an initial model of the
 * hypothalamic-pituitary-adrenal system self-regulatory activity.
 * A four-dimensional non-linear differential equation model of
 * the hormone secretion was formulated and used to analyze plasma
 * cortisol levels in humans. The aim of this work was to explore
 * in greater detail the role of this system in normal, homeostatic,
 * conditions, since it is the first and unavoidable step in further
 * understanding of the role of this complex neuroendocrine system
 * in pathophysiological conditions. Neither the underlying mechanisms
 * nor the physiological significance of this system are fully
 * understood yet.
 * 
 * The original paper reference is cited below:
 * 
 * Mathematical modeling of the hypothalamic-pituitary-adrenal
 * system activity, Smiljana Jelic, Zeljko Cupic, and Ljiljana
 * Kolar-Anic, 2005, Mathematical Biosciences, 197, pp173-187.
 * PubMed ID: 16112688
 * 
 * model diagram
 * 
 * [[Image file: jelic_2005.png]]
 * 
 * A schematic diagram of the hypothalamic-pituitary-adrenal system.
 * Corticotrophin-releasing hormone (CRH) and arginin-vasopressin
 * (AVP) stimulate adrenocorticotropin (ACTH) secretion from the
 * pituitary, followed by cortisol secretion from the adrenal cortex.
 * GR represents a glucocorticoid receptor and MR denotes a mineralocorticoid
 * receptor. + represents a postive feedback loop while - represents
 * a negative feedback loop.
 */

import nsrunit;
unit conversion on;
unit per_second=1 second^(-1);
// unit molar predefined
unit molar_per_second=1E3 meter^(-3)*second^(-1)*mole^1;
unit second_per_molar=.001 meter^3*second^1*mole^(-1);
unit per_molar2_per_second=1E-6 meter^6*second^(-1)*mole^(-2);

math main {
	realDomain tau dimensionless;
	tau.min=0;
	extern tau.max;
	extern tau.delta;
	real time(tau) second;
	real C_0 per_second;
	real A(tau) molar;
	real a(tau) dimensionless;
	when(tau=tau.min) a=4.39927;
	real alpha dimensionless;
	alpha=0.008;
	real beta dimensionless;
	beta=1.485;
	real KP(tau) dimensionless;
	real C_1 molar;
	real g(tau) dimensionless;
	when(tau=tau.min) g=1.96477;
	real G(tau) molar;
	real gamma dimensionless;
	gamma=11.385;
	real L dimensionless;
	L=0.1;
	real C_2 molar;
	real k2 per_second;
	k2=6e-4;
	real k3 per_second;
	k3=0.0000048;
	real k6 per_second;
	k6=0.000891;
	real k7 per_second;
	k7=0.006831;
	real k0 molar_per_second;
	k0=8.7831e-11;
	real k4 per_molar2_per_second;
	k4=2.1e12;
	real km molar_per_second;
	km=6.9001e-14;
	real K dimensionless;
	K=30;
	real P(tau) dimensionless;

	// <component name="environment">
	time=(tau/C_0);

	// <component name="a">
	A=(C_1*a);
	a:tau=(KP-((1+alpha+beta)*a+a*g^2));

	// <component name="g">
	G=(C_2*g);
	g:tau=((1-alpha)*a+a*g^2-(L+gamma*g));

	// <component name="reaction_constants">
	C_0=k2;
	C_1=sqrt(k2/k4);
	C_2=sqrt(k2/k4);
	KP=(K*P);
	P=(.735+.1563*sin(2*3.141592653*time/(24*(3600 second)))+.3127*abs(sin(3.141592653*time/(24*(3600 second)))));
}