/*
 * Modelling fluctuation phenomena in the plasma cortisol secretion
 * system in normal man
 * 
 * Model Status
 * 
 * This CellML model runs in OpenCell and COR, and the units are
 * consistent. The model is an accurate representation of the published
 * model (equations 1a-1c), and runs to recreate the published
 * results (figure 4). Parameter values were taken from the legend
 * of figure 4.
 * 
 * Model Structure
 * 
 * ABSTRACT: A system of three non-linear differential equations
 * with exponential feedback terms is proposed to model the self-regulating
 * cortisol secretion system and explain the fluctuation patterns
 * observed in clinical data. It is shown that the model exhibits
 * bifurcation and chaos patterns for a certain range of parametric
 * values. This helps us to explain clinical observations and characterize
 * different dynamic behaviors of the self-regulative system.
 * 
 * The original paper reference is cited below:
 * 
 * Modelling fluctuation phenomena in the plasma cortisol secretion
 * system in normal man, Yongwimon Lenbury and Pariwatana Pacheenburawana,
 * 1991, BioSystems , 26, 117-125. PubMed ID: 1668715
 * 
 * model diagram
 * 
 * [[Image file: lenbury_1991.png]]
 * 
 * Schematic diagram of the self-regulatory system for corticol
 * secretion. CRH represents corticotropin-releasing hormone, ACTH
 * is adrenocorticotropic hormone, and F represents cortisol. The
 * red lines represent positive feedback pathways, while the blue
 * lines represent negative feedback loops.
 */

import nsrunit;
unit conversion on;
unit per_second=1 second^(-1);

math main {
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real x(time) dimensionless;
	when(time=time.min) x=1.5;
	real alpha per_second;
	alpha=0.5;
	real omega per_second;
	omega=2;
	real D per_second;
	D=0.8228;
	real a dimensionless;
	a=8.1252;
	real b dimensionless;
	b=1.091;
	real z(time) dimensionless;
	when(time=time.min) z=1.0;
	real y(time) dimensionless;
	when(time=time.min) y=0.95;
	real beta per_second;
	beta=0.38;
	real gamma per_second;
	gamma=0.6;

	// <component name="environment">

	// <component name="x">
	x:time=(alpha*exp(a*(1-z^2)+b*(1-y^2))+D*cos(omega*time)-alpha*x);

	// <component name="y">
	y:time=(beta*x*exp(b*(1-z^2))-beta*y);

	// <component name="z">
	z:time=(gamma*y-gamma*z);

	// <component name="model_constants">
}