/*
 * A mathematical model of luteinizing hormone release from ovine
 * pituitary cells in perifusion
 * 
 * Model Status
 * 
 * The CellML model presented here represents the third model from
 * the published paper which describes the basic dynamics of GnRH-receptor
 * binding in the pituitary and the subsequent release of LH (as
 * in the first core model) however in this third model the constant
 * addition of LH to the releasable pool is removed such that it
 * becomes exhaustable (as in the second model) and also the desensitised
 * receptor is removed. The CellML model runs in COR and OpenCell
 * to replicate the published results. The units have been checked
 * and they are consistent.
 * 
 * Model Structure
 * 
 * ABSTRACT: We model the effect of gonadotropin-releasing hormone
 * (GnRH) on the production of luteinizing hormone (LH) by the
 * ovine pituitary. GnRH, released by the hypothalamus, stimulates
 * the secretion of LH from the pituitary. If stimulus pulses are
 * regular, LH response will follow a similar pattern. However,
 * during application of GnRH at high frequencies or concentrations
 * or with continuous application, the pituitary delivers a decreased
 * release of LH (termed desensitization). The proposed mathematical
 * model consists of a system of nonlinear differential equations
 * and incorporates two possible mechanisms to account for this
 * observed behavior: desensitized receptor and limited, available
 * LH. Desensitization was provoked experimentally in vitro by
 * using ovine pituitary cells in a perifusion system. The model
 * was fit to resulting experimental data by using maximum-likelihood
 * estimation. Consideration of smaller models revealed that the
 * desensitized receptor is significant. Limited, available LH
 * was significant in three of four chambers. Throughout, the proposed
 * model was in excellent agreement with experimental data.
 * 
 * model diagram
 * 
 * [[Image file: heinze_1998c.png]]
 * 
 * Schematic diagram of the components and reactions involved in
 * the first model of luteinizing hormone (LH) release. kfb and
 * kbf are kinetic constants, F and B represent the free and bound
 * states of the gonadotropin-releasing hormone (GnRH) receptor,
 * while R represents releasable LH and B is bound LH. s is the
 * rate of the reaction, a1 determines the rate of basal LH secretion
 * and a2 is the rate of LH secretion in the presence of bound
 * receptor. There is no desensitised GnRH receptor in this model.
 * 
 * The original paper reference is cited below:
 * 
 * A mathematical model of luteinizing hormone release from ovine
 * pituitary cells in perifusion, K. Heinze, R. W. Keener, and
 * A. R. Midgley, Jr., 1998, American Journal of Physiology, 275,
 * E1061-E1071. PubMed ID: 9843750
 */

import nsrunit;
unit conversion on;
unit hour=3600 second^1;
unit ng=1E-12 kilogram^1;
unit ng_ml=1E-6 kilogram^1*meter^(-3);
// unit nanomolar predefined
unit first_order_rate_constant=2.7777778E-4 second^(-1);
unit second_order_rate_constant=2.7777778E2 meter^3*second^(-1)*mole^(-1);

math main {
	realDomain time hour;
	time.min=0;
	extern time.max;
	extern time.delta;
	real GnRH(time) nanomolar;
	real F(time) dimensionless;
	when(time=time.min) F=1.0;
	real kfb second_order_rate_constant;
	kfb=19.35;
	real kbf first_order_rate_constant;
	kbf=9.91;
	real B(time) dimensionless;
	when(time=time.min) B=0.0;
	real R(time) dimensionless;
	when(time=time.min) R=2115.0;
	real s first_order_rate_constant;
	s=6.80;
	real a1 first_order_rate_constant;
	a1=0.0328;
	real a2 first_order_rate_constant;
	a2=0.0830;
	real C(time) dimensionless;
	when(time=time.min) C=0.0;

	// <component name="environment">

	// <component name="GnRH">
	GnRH=(if ((time>=(0 hour)) and (time<(.0666667 hour))) (.5 nanomolar) else if ((time>=(.0666667 hour)) and (time<(.4 hour))) (0 nanomolar) else if ((time>=(.4 hour)) and (time<(.4666667 hour))) (.5 nanomolar) else if ((time>=(.4666667 hour)) and (time<(2.4666667 hour))) (0 nanomolar) else if ((time>=(2.4666667 hour)) and (time<(2.533333 hour))) (.5 nanomolar) else if ((time>=(2.533333 hour)) and (time<(2.6166667 hour))) (0 nanomolar) else if ((time>=(2.6166667 hour)) and (time<(2.6833333 hour))) (.5 nanomolar) else if ((time>=(2.6833333 hour)) and (time<(4.6833333 hour))) (0 nanomolar) else if ((time>=(4.6833333 hour)) and (time<(4.75 hour))) (.5 nanomolar) else if ((time>=(4.75 hour)) and (time<(4.916667 hour))) (0 nanomolar) else if ((time>=(4.916667 hour)) and (time<(4.983333 hour))) (.5 nanomolar) else if ((time>=(4.983333 hour)) and (time<(6.983333 hour))) (0 nanomolar) else if ((time>=(6.983333 hour)) and (time<(7.066667 hour))) (.5 nanomolar) else if ((time>=(7.066667 hour)) and (time<(7.733333 hour))) (0 nanomolar) else if ((time>=(7.733333 hour)) and (time<(7.8 hour))) (.5 nanomolar) else if ((time>=(7.8 hour)) and (time<(9.8 hour))) (0 nanomolar) else 0);

	// <component name="F">
	F:time=(kbf*B-kfb*F*GnRH);

	// <component name="B">
	B:time=(kfb*F*GnRH-kbf*B);

	// <component name="R">
	R:time=(s-(a1+a2*B)*R);

	// <component name="C">
	C:time=((a1+a2*B)*R);

	// <component name="model_parameters">
}