/*
 * IP3-Mediated Ca2+ Release
 * 
 * Model Status
 * 
 * This model is able to be integrated in OpenCell but is not currently
 * valid CellML.
 * 
 * Model Structure
 * 
 * Ca2+ is a ubiquitous intracellular secondary messenger, and
 * evidence from several different cell types suggests that an
 * important mode of signalling is through oscillations rather
 * than the maintenance of a steady state level. The oscillatory
 * behaviour of inositol 1,4,5-triphosphate (IP3)-mediated Ca2+
 * release has been modelled by Gary W. De Young and Joel Keizer.
 * Their 1992 paper is referenced fully below.
 * 
 * >A single-pool inositol 1,4,5-triphosphate-receptor-based model
 * for agonist-stimulated oscillations in Ca2+ concentration, Gary
 * W. De Young and Joel Keizer, 1992, Proc. Natl. Acad. Sci. USA
 * , 89, 9895-9899.PubMed ID: 1329108
 * 
 * Several mechanisms have been proposed to explain oscillations
 * of intracellular Ca2+ concentration in cells. In this study,
 * De Young and Keizer investigate the idea that a biphasic response
 * of the IP3 receptor/channel to cytosolic Ca2+ may alone be sufficient
 * to induce Ca2+ oscillations.
 * 
 * They constructed a simplified model of the IP3 receptor/channel
 * by assuming that three equivalent and independent subunits are
 * involved in Ca2+ conduction. Each subunit has three binding
 * sites: one for IP3, one for Ca2+ activation, and one for Ca2+
 * inactivation. Thus each subunit may exist in eight states with
 * transitions governed by second-order (association) and first-order
 * (dissociation) rate constants (see below). All three subunits
 * must be in the state S110 (one IP3 and one activating Ca2+ bound)
 * for the channel to be open and conducting.
 * 
 * A schematic diagram of the kinetics of an IP3 receptor/channel
 * subunit
 * 
 * [[Image file: deyoung_1992.png]]
 * 
 * A schematic diagram of the kinetics of an IP3 receptor/channel
 * subunit.
 */

import nsrunit;
// Warning: unit conversion turned off due to unit errors in 4 equation(s)
unit conversion off;
// unit micromolar predefined
// unit nanomolar 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:  Ca_i, IP3, x_000, x_001, x_010
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real Ca_i(time) micromolar;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) Ca_i=0;
	real J1(time) flux;
	real J2(time) flux;
	real v1 first_order_rate_constant;
	v1=6.0;
	real v2 first_order_rate_constant;
	v2=0.11;
	real v3 second_order_rate_constant;
	v3=0.9;
	real k3 micromolar;
	k3=0.1;
	real c1 dimensionless;
	c1=0.185;
	real Ca_ER(time) micromolar;
	real P_open(time) dimensionless;
	real c0 micromolar;
	c0=2.0;
	real IP3(time) micromolar;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) IP3=0;
	real k4 micromolar;
	k4=1.1;
	real alpha dimensionless;
	alpha=0.5;
	real Ir flux;
	Ir=1.0;
	real v4 first_order_rate_constant;
	v4=1.2;
	real d1 micromolar;
	d1=0.13;
	real d2 micromolar;
	d2=1.049;
	real d3 nanomolar;
	d3=943.4;
	real d4 nanomolar;
	real d5 nanomolar;
	d5=82.34;
	real a5 second_order_rate_constant;
	a5=20.0;
	real b5 first_order_rate_constant;
	real a2 second_order_rate_constant;
	a2=0.2;
	real b2 first_order_rate_constant;
	real a3 second_order_rate_constant;
	a3=400.0;
	real b3 first_order_rate_constant;
	real a4 second_order_rate_constant;
	a4=0.2;
	real b4 first_order_rate_constant;
	real IP3_cold nanomolar;
	IP3_cold=15.0;
	real K_d2 nanomolar;
	real K_d1 nanomolar;
	real x_000(time) micromolar;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) x_000=0;
	real V1(time) flux;
	real V3(time) flux;
	real x_001(time) micromolar;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) x_001=0;
	real V4(time) flux;
	real x_010(time) micromolar;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) x_010=0;
	real V2(time) flux;
	real x_011(time) micromolar;
	real a1 second_order_rate_constant;
	a1=400.0;
	real b1 first_order_rate_constant;

	// <component name="environment">

	// <component name="Ca_i">
	Ca_i:time=(J1-J2);
	J1=(c1*(v1*P_open+v2)*(Ca_ER-Ca_i));
	J2=(v3*Ca_i^2/(Ca_i^2+k3^2));

	// <component name="Ca_ER">
	Ca_ER=((c0-Ca_i)/c1);

	// <component name="IP3">
	IP3:time=(v4*((Ca_i+(1-alpha)*k4)/(Ca_i+k4))-Ir*IP3);

	// <component name="receptor_dissociation_constants">
	d1=(K_d1-IP3_cold);
	d2=(b2/a2);
	d3=((K_d2-IP3_cold)*(1+d2)-d1*d2);
	d4=(d1*d2/d3);
	d5=(b5/a5);

	// <component name="probability_of_an_open_IP3_receptor_channel">
	P_open=((Ca_i*IP3*d2/((Ca_i*IP3+IP3*d2+d1*d2+Ca_i*d3)*(Ca_i+d5)))^3);

	// <component name="x_000">
	x_000:time=((-1)*V1-V3);

	// <component name="x_001">
	x_001:time=(V1-V4);

	// <component name="x_010">
	x_010:time=(V3-V2);

	// <component name="x_011">
	x_011=(1-(x_000+x_001+x_010));

	// <component name="fluxes">
	V1=(a4*(Ca_i*x_000-d4*x_001));
	V2=(a4*(Ca_i*x_010-d4*x_011));
	V3=(a5*(Ca_i*x_000-d5*x_010));
	V4=(a5*(Ca_i*x_001-d5*x_011));

	// <component name="constants">
	b1=(d1*a1);
	b3=(d3*a3);
	b4=(d4*a4);
}