/*
* 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;
//
//
H1:time=((-1)*(beta_I*H1)+A);
//
H2:time=(beta_I*H1-beta_C*H2);
//
A=(if ((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));
//
beta_I=(ln(2)/tau_I);
beta_C=(ln(2)/tau_C);
alpha=(ln(2)/tau_alpha);
lamda=(ln(2)/tau_lamda);
}