/*
 * A computer simulation of the hypothalamic-pituitary-adrenal
 * axis
 * 
 * Model Status
 * 
 * This model runs in OpenCell and COR. The units have been checked
 * and they are consistent. Please note this CellML translation
 * was based on a an online version of the original paper - a version
 * which contains more defined parameter values than the paper
 * linked to from PubMed (the more detailed description). Also
 * there is a typo in this description in the ACTH Module. There
 * should be a minus sign rather than a plus sign infront of (k4*[CRH]+k5)
 * (We thank the original model author for highlighting this).
 * The correct equation is as follows: d[ACTH]/dt = k4*[CRH]+k5
 * - k6*[ACTH]- (k4*[CRH]+k5)*Imax*[cortisol]/(Kd+[cortisol]).
 * Unfortunately there are no published figures of the results
 * of the original model so we cannot confirm the CellML translation
 * recreates the original published results.
 * 
 * Model Structure
 * 
 * ABSTRACT: This paper describes the construction of a computer
 * model that simulates the hypothalamic-pituitary-adrenal axis
 * (HPA axis) regulation of cortisol production. It is presented
 * to illustrate the process of physiological modeling using standard
 * "off the shelf" technologies. The model simulates components
 * of the HPA axis involved in the continuous secretion and elimination
 * of cortisol, adrenocorticotropin (ACTH), and corticotropin releasing
 * hormone (CRH). The physiological relations of these component
 * pieces were modeled based on the current knowledge of their
 * functioning. Rate constants, half lives, and receptor affinities
 * were assigned values derived from the experimental literature.
 * At its current level of development the model is able to accurately
 * simulate the timing, magnitude and decay of the ACTH and cortisol
 * concentration peaks resulting from the ovine-CRH stimulation
 * test in normal and hypercortisolemic patients. The model will
 * be used to predict the effects of lesions in different components
 * of the HPA axis on the time course of cortisol and ACTH levels.
 * We plan to use the model to explore the experimental conditions
 * required to distinguish mechanisms underlying various disorders
 * of the HPA axis, particularly depression. Efforts are currently
 * underway to validate the model for a large variety of normal
 * and pathological perturbations of the HPA axis.
 * 
 * The complete paper reference is cited below:
 * 
 * A computer simulation of the hypothalamic-pituitary-adrenal
 * axis, J. Gonzalez-Heydrich, R.J. Steingard and I Kohane, 1994,
 * Eighteenth Annual Proceedings of the Symposium for Computer
 * Applications in Medical Care, pp1010. PubMed ID: 7949852
 * 
 * model diagram
 * 
 * [[Image file: gonzalezheidrich_1994.png]]
 * 
 * Schematic diagram of the self-regulatory system for corticol
 * secretion. CRH represents corticotropin-releasing hormone, and
 * ACTH is adrenocorticotropic hormone.
 */

import nsrunit;
unit conversion on;
unit mcg_ml=1E-3 kilogram^1*meter^(-3);
unit minute=60 second^1;
unit flux=1.6666667E-5 kilogram^1*meter^(-3)*second^(-1);
unit first_order_rate_constant=.01666667 second^(-1);

math main {
	realDomain time minute;
	time.min=0;
	extern time.max;
	extern time.delta;
	real cortisol(time) mcg_ml;
	when(time=time.min) cortisol=0.0;
	real k1 first_order_rate_constant;
	k1=5E5;
	real k2 flux;
	k2=0.01;
	real k3 first_order_rate_constant;
	k3=0.01;
	real ACTH(time) mcg_ml;
	when(time=time.min) ACTH=0.0;
	real k4 first_order_rate_constant;
	k4=10.0;
	real k5 flux;
	k5=4E-3;
	real k6 first_order_rate_constant;
	k6=0.035;
	real Kd mcg_ml;
	Kd=0.004;
	real Imax dimensionless;
	Imax=0.99;
	real CRH(time) mcg_ml;
	when(time=time.min) CRH=50.0;
	real k7 flux;
	k7=1E-6;
	real k8 first_order_rate_constant;
	k8=0.01;
	real pulse flux;
	pulse=50.0;

	// <component name="environment">

	// <component name="cortisol">
	cortisol:time=(k1*ACTH+k2-k3*cortisol);

	// <component name="ACTH">
	ACTH:time=(k4*CRH+k5-(k6*ACTH+(k4*CRH+k5)*(Imax*cortisol/(Kd+cortisol))));

	// <component name="CRH">
	CRH:time=(pulse+k7-k8*CRH);
}