/*
 * A mathematical model of diurnal variations in human plasma melatonin
 * levels
 * 
 * Model Status
 * 
 * This CellML model runs in OpenCell to recreate the published
 * results. The units have been checked and they are consistent.
 * The model also runs in COR but due to the time being in minutes
 * and the simulation time being 3000 minutes long the model is
 * not ideal for running in COR.
 * 
 * Model Structure
 * 
 * ABSTRACT: Studies in animals and humans suggest that the diurnal
 * pattern in plasma melatonin levels is due to the hormone's rates
 * of synthesis, circulatory infusion and clearance, circadian
 * control of synthesis onset and offset, environmental lighting
 * conditions, and error in the melatonin immunoassay. A two-dimensional
 * linear differential equation model of the hormone is formulated
 * and is used to analyze plasma melatonin levels in 18 normal
 * healthy male subjects during a constant routine. Recently developed
 * Bayesian statistical procedures are used to incorporate correctly
 * the magnitude of the immunoassay error into the analysis. The
 * estimated parameters [median (range)] were clearance half-life
 * of 23.67 (14.79-59.93) min, synthesis onset time of 2206 (1940-0029),
 * synthesis offset time of 0621 (0246-0817), and maximum N-acetyltransferase
 * activity of 7.17(2.34-17.93) pmol x l(-1) x min(-1). All were
 * in good agreement with values from previous reports. The difference
 * between synthesis offset time and the phase of the core temperature
 * minimum was 1 h 15 min (-4 h 38 min-2 h 43 min). The correlation
 * between synthesis onset and the dim light melatonin onset was
 * 0.93. Our model provides a more physiologically plausible estimate
 * of the melatonin synthesis onset time than that given by the
 * dim light melatonin onset and the first reliable means of estimating
 * the phase of synthesis offset. Our analysis shows that the circadian
 * and pharmacokinetics parameters of melatonin can be reliably
 * estimated from a single model.
 * 
 * model diagram
 * 
 * [[Image file: brown_1997.png]]
 * 
 * Schematic diagram of the components and variables in a mathematical
 * model of human melatonin. Diurnal rhythm of melatonin is observed
 * in the plasma compartment. A represents N-acetylserotonin activity
 * in the pineal gland, H1 is the concentration of melatonin in
 * the pineal gland, beta_I is the plasma melatonin infusion rate,
 * H2 is the plasma concentration of melatonin, and beta_C is the
 * clearance rate of plasma melatonin.
 * 
 * The original paper reference is cited below:
 * 
 * A mathematical model of diurnal variations in human plasma melatonin
 * levels, Emery N. Brown, Yong Choe, Theresa L. Shanahan, and
 * Charles A. Czeisler, 1997, American Journal of Physiology; Endocrinology
 * and Metabolism, 272, E506-E516. PubMed ID: 9124558
 */

import nsrunit;
unit conversion on;
unit minute=60 second^1;
// unit picomolar predefined
unit flux=1.6666667E-11 meter^(-3)*second^(-1)*mole^1;
unit first_order_rate_constant=.01666667 second^(-1);

math main {
	realDomain time minute;
	time.min=0;
	extern time.max;
	extern time.delta;
	real H1(time) picomolar;
	when(time=time.min) H1=0.05;
	real A(time) flux;
	real beta_I first_order_rate_constant;
	real H2(time) picomolar;
	when(time=time.min) H2=1.0;
	real beta_C first_order_rate_constant;
	real t_on minute;
	t_on=1316.0;
	real t_off minute;
	t_off=1792.0;
	real A_max flux;
	A_max=6.51;
	real alpha first_order_rate_constant;
	real lamda first_order_rate_constant;
	real tau_I minute;
	tau_I=2.82;
	real tau_C minute;
	tau_C=23.67;
	real tau_alpha minute;
	tau_alpha=25.92;
	real tau_lamda minute;
	tau_lamda=24.04;

	// <component name="environment">

	// <component name="H1">
	H1:time=((-1)*(beta_I*H1)+A);

	// <component name="H2">
	H2:time=(beta_I*H1-beta_C*H2);

	// <component name="A">
	A=(if ((time<t_off) and (time>=t_on)) A_max*((1-exp((-1)*lamda*(time-t_on)))/(1-exp((-1)*lamda*(t_off-t_on)))) else if (time>=t_off) A_max*exp((-1)*alpha*(time-t_off)) else (0 flux));

	// <component name="model_parameters">
	beta_I=(ln(2)/tau_I);
	beta_C=(ln(2)/tau_C);
	alpha=(ln(2)/tau_alpha);
	lamda=(ln(2)/tau_lamda);
}