/*
 * IP3-Mediated Ca2+ Oscillations - A Simplified Model
 * 
 * Model Status
 * 
 * This is the original unchecked version of the model imported
 * from the previous CellML model repository, 24-Jan-2006.
 * 
 * Model Structure
 * 
 * Cytosolic calcium ([Ca2+]i) shows complex spatio-temporal oscillatory
 * profiles in a large number of cell types. A number of mathematical
 * models have been published that attempt to model the temporal
 * [Ca2+]i oscillations. Some, such as the Keizer and De Young
 * model (see The De Young-Keizer model, 1992) assume that IP3
 * is necessarily oscillatory. They constructed a 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.
 * 
 * In 1994, Yue-Xian Li and John Rinzel analysed the nine-variable
 * De Young-Keizer model (1992) and reduced it to a two-variable
 * system. Their reduced system is analogous in form to the classic
 * Hodgkin-Huxley equations (see The Hodgkin-Huxley Squid Axon
 * Model, 1952) for plasma membrane electrical excitability. [Ca2+]i
 * dynamics in this model thus involve endoplasmic reticulum (ER)
 * membrane-associated excitability. Assuming further that the
 * binging of IP3 does not depend on Ca2+ occupancy at the inactivation
 * site, Li and Rinzel obtain a "minimal" model which still retains
 * the ability to reproduce experimental observations.
 * 
 * The complete original paper reference is cited below:
 * 
 * Equations for InsP3 Receptor-mediated [Ca2+]i Oscillations Derived
 * from a Detailed Kinetic Model: A Hodgkin-Huxley Like Formalism,
 * Yue-Xian Li and John Rinzel, 1994, Journal of Theoretical Biology,
 * 166, 461-473. PubMed ID: 8176949
 * 
 * In addition to the original version of the Li and Rinzel mathematical
 * model, the system described in equation nine of the paper has
 * also been implemented in CellML as a separate model (see for
 * the downloads of both models in various formats). This is a
 * set of three coupled differential equations, namely:
 * 
 * One equation for describing the calcium dynamics;
 * 
 * One equation which represents the gate-shutting variable h;
 * and
 * 
 * One equation which describes the influx of extracellular calcium
 * into the cell cytosol through the plasma membrane
 * 
 * The advantage of reducing the model to this set of three equations
 * is that it more simple, yet it still retains the ability to
 * explain oscillatory behaviours caused by the calcium dynamics.
 * All the model parameters have been set according to those which
 * are provided in Figure 5(a) of the paper.
 * 
 * The raw CellML descriptions of the Li and Rinzel 1994 model
 * can be downloaded in various formats as described in .
 * 
 * A schematic diagram of the kinetics of an IP3 receptor/channel
 * subnuit
 * 
 * [[Image file: li_1994.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 10 equation(s)
unit conversion off;
// unit micromolar predefined
// unit nanomolar predefined
unit flux=1 meter^(-3)*second^(-1)*mole^1;
unit first_order_rate_constant=1E3 second^(-1);
unit second_order_rate_constant=1E6 meter^3*second^(-1)*mole^(-1);
unit ms=.001 second^1;

math main {
	//Warning:  the following variables had initial values which were
	//  suppressed because the model would otherwise be overdetermined:
	//   factor1, factor2, v3_summand, v4_summand
	realDomain time ms;
	time.min=0;
	extern time.max;
	extern time.delta;
	real C(time) micromolar;
	when(time=time.min) C=0.2;
	real a first_order_rate_constant;
	a=1.;
	real v1 first_order_rate_constant;
	v1=40.;
	real v2(time) first_order_rate_constant;
	when(time=time.min) v2=0.02;
	real v3(time) second_order_rate_constant;
	when(time=time.min) v3=0.6;
	real v4 second_order_rate_constant;
	v4=1.8;
	real d_ip3 micromolar;
	d_ip3=0.2;
	real d_act micromolar;
	d_act=0.4;
	real d_inh micromolar;
	d_inh=0.4;
	real k_er micromolar;
	k_er=0.18;
	real k_pl micromolar;
	k_pl=0.1;
	real c0(time) micromolar;
	when(time=time.min) c0=2.0;
	real IP3 micromolar;
	IP3=0.6;
	real j_in(time) flux;
	when(time=time.min) j_in=1.2;
	real factor1(time) flux;
	//Warning: CellML initial value suppressed to prevent overdetermining model.  Original initial value: factor1=0.0;
	real factor2(time) flux;
	//Warning: CellML initial value suppressed to prevent overdetermining model.  Original initial value: factor2=0.0;
	real v3_summand(time) flux;
	//Warning: CellML initial value suppressed to prevent overdetermining model.  Original initial value: v3_summand=0.0;
	real v4_summand(time) flux;
	//Warning: CellML initial value suppressed to prevent overdetermining model.  Original initial value: v4_summand=0.0;
	real epsilon dimensionless;
	epsilon=0.01;
	real c1 dimensionless;
	c1=0.185;
	real h(time) dimensionless;
	when(time=time.min) h=0.8;

	// <component name="environment">

	// <component name="C">
	v2:time=(0 flux);
	v3:time=(0 flux);
	j_in:time=(0 flux);
	c0:time=(epsilon*.001*(j_in-v4*C^2/(k_pl^2+C^2)));
	h:time=(.001*(a*d_inh-h*a*(C+d_inh)));
	factor1=(.001*(v1*(IP3/(IP3+d_ip3))^3*(C/(C+d_act))^3*h^3+v2));
	factor2=(c0-C*(1+c1));
	v3_summand=(v3*.001*C^2/(k_er^2+C^2));
	v4_summand=(.001*epsilon*(j_in-v4*C^2/(k_pl^2+C^2)));
	C:time=(factor1*factor2-v3_summand+v4_summand);
}