/*
 * Agonist-Dependent Phosphorylation Of The Inositol 1,4,5-Triphosphate
 * Receptor
 * 
 * Model Structure
 * 
 * The production of the intracellular signalling factor inositol
 * 1,4,5-triphosphate (IP3) and the subsequent release of Ca2+
 * stored in intracellular organelles is a fundamental cellular
 * signalling function. The inositol triphosphate receptor (IPR)
 * is an IP3-activated Ca2+ channel in the ER. The properties of
 * IP3-dependent intracellular calcium oscillations in pancreatic
 * acinar cells depend on the agonist used to stimulate them.
 * 
 * This agonist-dependency is captured in Andrew P. LeBeau et al's
 * 1999 mathematical model of agonist-specific intracellular calcium
 * oscillations in pancreatic acinar cells. They assume that the
 * complete IPR is composed of four functionally identical, independent
 * subunits, (see below). IP3 must be bound to all four subunits
 * for the receptor to be in the conducting state.
 * 
 * The complete original paper reference is cited below:
 * 
 * Agonist-dependent Phosphorylation of the Inositol 1,4,5-Triphosphate
 * Receptor. A Possible Mechanism for Agonist-specific Calcium
 * Oscillations in Pancreatic Acinar Cells, Andrew P. LeBeau, David
 * I. Yule, Guy E. Groblewski and James Sneyd, 1999, The Journal
 * Of General Physiology, 113, 851-871. PubMed ID: 10352035
 * 
 * The raw CellML description of the dynamic model of the type-2
 * inositol triphosphate receptor can be downloaded in various
 * formats as described in .
 * 
 * A simplified diagram of the IPR model
 * 
 * [[Image file: lebeau_1999.png]]
 * 
 * A diagram of the receptor states of the model of the IP3 receptor.
 * S denotes the fraction of the subunits in the shut state. Binding
 * of IP3 causes the receptor to be converted to the open state
 * O. O is a relatively unstable state and the subunits will progress
 * through to the more stable I1 (inactivated) state in which IP3
 * is still bound but the channels do not conduct. I2 represents
 * a second inactivated state of the receptor in which IP3 is no
 * longer bound.
 */

import nsrunit;
// Warning: unit conversion turned off due to unit errors in 7 equation(s)
unit conversion off;
// unit micromolar predefined
unit flux=1E-3 meter^(-3)*second^(-1)*mole^1;
unit first_order_rate_constant=1 second^(-1);
unit second_order_rate_constant=1E3 meter^3*second^(-1)*mole^(-1);

math main {
	//Warning:  the following variables were set 'extern' or given
	//  an initial value of '0' because the model would otherwise be
	//  underdetermined:  O, I1, I2, c
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real O(time) micromolar;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) O=0;
	real S(time) micromolar;
	real I1(time) micromolar;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) I1=0;
	real I2(time) dimensionless;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) I2=0;
	real p micromolar;
	p=3.0;
	real k1(time) first_order_rate_constant;
	real k1_ first_order_rate_constant;
	k1_=0.88;
	real k2 first_order_rate_constant;
	k2=0.5;
	real k3 first_order_rate_constant;
	k3=0.5;
	real k4 first_order_rate_constant;
	real k5 first_order_rate_constant;
	k5=0.02;
	real c(time) micromolar;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) c=0;
	real J_rel(time) flux;
	real J_pump(time) flux;
	real J_influx flux;
	J_influx=0.4;
	real V_p flux;
	V_p=2.6;
	real K_p micromolar;
	K_p=0.54;
	real k_flux flux;
	k_flux=28.0;
	real alpha_1 flux;
	alpha_1=40.0;
	real beta_1 micromolar;
	beta_1=0.8;
	real alpha_4 flux;
	alpha_4=0.06;
	real beta_4 micromolar;
	beta_4=0.01;

	// <component name="environment">

	// <component name="O">
	O:time=(k1*p*S-(k1_*O+k2*O));
	S=(4-(O+I1+I2));

	// <component name="I1">
	I1:time=(k2*O-(k3+k4)*I1);

	// <component name="I2">
	I2:time=(k4*I1-k5*I2);

	// <component name="c">
	c:time=(J_rel-J_pump+J_influx);
	J_rel=(k_flux*O^4);
	J_pump=(V_p*c^2/(K_p^2+c^2));

	// <component name="constants">
	k1=(alpha_1*c^3/(beta_1^3+c^3));
	k4=(alpha_4*p/(beta_4+p));
}