/*
 * Modelling the glucose-insulin regulatory system and ultradian
 * insulin secretory oscillations
 * 
 * Model Status
 * 
 * This CellML model runs in both COR and OpenCell however it does
 * not recreate the published results. Glucose and insulin concentrations
 * do not oscillate over time. This is most likely to be due to
 * the lack of time delays in the CellML model (which are currently
 * impossible to represent in CellML, but are required for equation
 * 1 in the paper). The units are consistent. Parameter values
 * have been taken from table 1 and di=0.06 and Gin=0.54 - as in
 * figure 4.
 * 
 * Model Structure
 * 
 * ABSTRACT: In the glucose-insulin regulatory system, ultradian
 * insulin secretory oscillations are observed to have a period
 * of 50-150 min. After pioneering work traced back to the 1960s,
 * several mathematical models have been proposed during the last
 * decade to model these ultradian oscillations as well as the
 * metabolic system producing them. These currently existing models
 * still lack some of the key physiological aspects of the glucose-insulin
 * system. Applying the mass conservation law, we introduce two
 * explicit time delays and propose a more robust alternative model
 * for better understanding the glucose-insulin endocrine metabolic
 * regulatory system and the ultradian insulin secretory oscillations
 * for the cases of continuous enteral nutrition and constant glucose
 * infusion. We compare the simulation profiles obtained from this
 * two time delay model with those from the other existing models.
 * As a result, we notice many unique features of this two delay
 * model. Based on our intensive simulations, we suspect that one
 * of the possibly many causes of ultradian insulin secretion oscillations
 * is the time delay of the insulin secretion stimulated by the
 * elevated glucose concentration.
 * 
 * model diagram
 * 
 * [[Image file: li_2006.png]]
 * 
 * Schematic diagram of the glucose-insulin regulatory system model.
 * The dotted lines indicate that elevated glucose concentration
 * stimulates insulin synthesis and secretion by the pancreatic
 * beta-cells, and also insulin promotes glucose utilisation in
 * muscle, fat and other tissues.
 * 
 * The original paper reference is cited below:
 * 
 * Modeling the glucose-insulin regulatory system and ultradian
 * insulin secretory oscillations with two explicit time delays,
 * Jiaxu Li, Yang Kuang, and Clinton C. Mason, 2006, Journal of
 * Theoretical Biology, 242, 722-735. PubMed ID: 16712872
 */

import nsrunit;
unit conversion on;
unit minute=60 second^1;
unit U = fundamental;
unit mU_min=1.6666667E-5 second^(-1)*U^1;
unit per_l=1E3 meter^(-3);
unit l_min=1.6666667E-5 meter^3*second^(-1);
unit l_mU=1 meter^3*U^(-1);
unit mg_l_min=1.6666667E-5 kilogram^1*meter^(-3)*second^(-1);
unit mg_l=1E-3 kilogram^1*meter^(-3);
unit mU_l=1 meter^(-3)*U^1;
unit mg_min=1.6666667E-8 kilogram^1*second^(-1);
unit first_order_rate_constant=.01666667 second^(-1);

math main {
	realDomain time minute;
	time.min=0;
	extern time.max;
	extern time.delta;
	real G(time) mg_l;
	when(time=time.min) G=1000.0;
	real Vg litre;
	Vg=10.0;
	real Gin mg_l_min;
	Gin=0.54;
	real f2_G(time) mg_min;
	real f3_G(time) per_l;
	real f4_I(time) mg_min;
	real f5_I(time) mg_min;
	real C2 mg_l;
	C2=144.0;
	real C3 mg_l;
	C3=1000.0;
	real C4 mU_l;
	C4=80.0;
	real C5 mU_l;
	C5=26.0;
	real U0 mg_min;
	U0=40.0;
	real Um mg_min;
	Um=940.0;
	real Ub mg_min;
	Ub=72.0;
	real beta dimensionless;
	beta=1.77;
	real Rg mg_min;
	Rg=180.0;
	real alpha l_mU;
	alpha=0.29;
	real ti minute;
	ti=100;
	real E l_min;
	E=0.2;
	real Vp litre;
	Vp=3.0;
	real Vi litre;
	Vi=11.0;
	real I(time) mU_l;
	when(time=time.min) I=9000.0;
	real f1_G(time) mU_min;
	real Rm mU_min;
	Rm=210.0;
	real C1 mg_l;
	C1=2000.0;
	real a1 mg_l;
	a1=300.0;
	real di first_order_rate_constant;
	di=0.06;

	// <component name="environment">

	// <component name="G">
	G:time=(Gin+(1 per_l)*f5_I+(-1)*((1 per_l)*f2_G+f3_G*f4_I));
	f2_G=(Ub*(1-exp((1 litre)*((-1)*G/(C2*Vg)))));
	f3_G=(G/(C3*Vg));
	f4_I=(U0+(Um-U0)/(1+exp((-1)*beta*ln((1 per_l)*(I/(C4*(1/Vi+1/(E*ti))))))));
	f5_I=(Rg/(1+exp(alpha*(I/((1 per_l)*Vp-(1 l_mU)*C5)))));

	// <component name="I">
	I:time=((1 per_l)*f1_G-di*I);
	f1_G=(Rm/(1+exp((C1-G)/((1 per_l)*Vg)/a1)));
}