/*
 * Mathematical modeling and qualitative analysis of insulin therapies
 * 
 * Model Status
 * 
 * This CellML model runs in PCenv, COR and OpenCell to recreate
 * the published results. The initial conditions were not explicitly
 * defined in the paper, but were extrapolated from Table 1 and
 * Figure 2. The model was then run until it reached a steady state,
 * and new intial values were taken from the steady state.
 * 
 * Model Structure
 * 
 * ABSTRACT: We propose a dynamical model for the pulsatile secretion
 * of the hypothalamo-pituitary-adrenal axis. This model takes
 * into account both the binding of hormone with proteins in the
 * plasma and tissues and mutual interactions of the hormones in
 * the system. The deductions from this model are in good agreement
 * with experiment results.
 * 
 * model diagram
 * 
 * [[Image file: liu_1999.png]]
 * 
 * Schematic diagram of the hormone flow and interactions between
 * the different components of the hypothalamic-pituitary-adrenal
 * (HPA) axis which are described by the mathematical model. CRH
 * represents corticotropin releasing hormone, ACTH represents
 * corticotropin, and CBG is corticosteroid binding globulin. Red
 * arrows represent hormone disposal, while black arrows represent
 * flow between the different glands and plasma. + represents a
 * stimulatory effect while - represents an inhibitory effect.
 * 
 * The original paper reference is cited below:
 * 
 * A dynamical model for the pulsatile secretion of the hypothalamo-pituitary-adrenal
 * axis, Yi-Wei Liu, Zhi-Hong Hu, Jian-Hua Peng, and Bing-Zheng
 * Liu, 1999, Mathematical and Computer Modelling, 29 (Issue 4),
 * 103-110. (No Pubmed ID)
 */

import nsrunit;
unit conversion on;
unit minute=60 second^1;
unit microg_l=1E-6 kilogram^1*meter^(-3);
unit per_microg_l=1E6 kilogram^(-1)*meter^3;
unit per_microg_l2=1E12 kilogram^(-2)*meter^6;
unit flux=1.6666667E-8 kilogram^1*meter^(-3)*second^(-1);
unit first_order_rate_constant=.01666667 second^(-1);
unit second_order_rate_constant=1.6666667E4 kilogram^(-1)*meter^3*second^(-1);

math main {
	realDomain time minute;
	time.min=0;
	extern time.max;
	extern time.delta;
	real x1(time) microg_l;
	when(time=time.min) x1=0.01067;
	real lambda_1 first_order_rate_constant;
	lambda_1=0.059;
	real x3(time) microg_l;
	when(time=time.min) x3=6.51;
	real a1 flux;
	a1=0.000017;
	real a2 flux;
	a2=0.0023;
	real a3 first_order_rate_constant;
	a3=0.6;
	real a4 second_order_rate_constant;
	a4=45;
	real a5 per_microg_l;
	a5=36;
	real a6 per_microg_l2;
	a6=216;
	real a7 per_microg_l;
	a7=0.28;
	real a8 per_microg_l2;
	a8=0.36;
	real x2(time) microg_l;
	when(time=time.min) x2=0.04665;
	real lambda_2 first_order_rate_constant;
	lambda_2=0.028;
	real a9 flux;
	a9=0.0003;
	real a10 first_order_rate_constant;
	a10=0.18;
	real a11 second_order_rate_constant;
	a11=150;
	real a12 per_microg_l;
	a12=18;
	real a13 per_microg_l2;
	a13=460;
	real a14 per_microg_l;
	a14=0.46;
	real a15 per_microg_l2;
	a15=0.1;
	real lambda_3_ first_order_rate_constant;
	real x4(time) microg_l;
	when(time=time.min) x4=60.61;
	real x5(time) microg_l;
	when(time=time.min) x5=12.61;
	real a16 flux;
	a16=0.04;
	real a17 first_order_rate_constant;
	a17=150;
	real a18 second_order_rate_constant;
	a18=3800;
	real a19 first_order_rate_constant;
	a19=57;
	real a20 second_order_rate_constant;
	a20=2600;
	real a21 per_microg_l;
	a21=200;
	real a22 per_microg_l2;
	a22=9400;
	real a23 per_microg_l;
	a23=10;
	real a24 per_microg_l2;
	a24=320;
	real a25 first_order_rate_constant;
	a25=0.04;
	real a26 first_order_rate_constant;
	a26=0.00097;
	real lambda_4_ first_order_rate_constant;
	real a27 first_order_rate_constant;
	a27=0.57;
	real lambda_5_ first_order_rate_constant;
	real a28 first_order_rate_constant;
	a28=0.0017;
	real lambda_3 first_order_rate_constant;
	lambda_3=0.0986;
	real lambda_4 first_order_rate_constant;
	lambda_4=0.024;
	real lambda_5 first_order_rate_constant;
	lambda_5=3e-5;

	// <component name="environment">

	// <component name="x1">
	x1:time=(a1+(a2+a3*x1+a4*x1^2)/(1+a5*x1+a6*x1^2+a7*x3+a8*x3^2)-lambda_1*x1);

	// <component name="x2">
	x2:time=((a9+a10*x1+a11*x1^2)/(1+a12*x1+a13*x1^2+a14*x3+a15*x3^2)-lambda_2*x2);

	// <component name="x3">
	x3:time=(a16+(a17*x1+a18*x1^2+a19*x2+a20*x2^2)/(1+a21*x1+a22*x1^2+a23*x2+a24*x2^2)+a25*x4+a26*x5-lambda_3_*x3);

	// <component name="x4">
	x4:time=(a27*x3-lambda_4_*x4);

	// <component name="x5">
	x5:time=(a28*x3-lambda_5_*x5);

	// <component name="model_parameters">
	lambda_3_=(lambda_3+a27+a28);
	lambda_4_=(lambda_4+a25);
	lambda_5_=(lambda_5+a26);
}