/*
 * Minimal model for signal-induced Ca2+ oscillations and for their
 * frequency encoding through protein phosphorylation
 * 
 * Model Status
 * 
 * Runs in PCEnv and COR. Reproduces published output (figures
 * 2a, 2b when the value of 'k' is changed to 6, 4a, 4b when the
 * value of 'beta' is changed to 0.644)
 * 
 * Model Structure
 * 
 * ABSTRACT: In a variety of cells, hormonal or neurotransmitter
 * signals elicit a train of intracellular Ca2+ spikes. The analysis
 * of a minimal model based on Ca2+-induced Ca2+ release from intracellular
 * stores shows how sustained oscillations of cytosolic Ca2+ may
 * develop as a result of a rise in inositol 1,4,5-trisphosphate
 * (InsP3) triggered by external stimulation. This rise elicits
 * the release of a certain amount of Ca2+ from an InsP3-sensitive
 * intracellular store. The subsequent rise in cytosolic Ca2+ in
 * turn triggers the release of Ca2+ from a second store insensitive
 * to InsP3. In contrast to the model proposed by Meyer and Stryer
 * [Meyer, T. & Stryer, L. (1988) Proc. Matl. Acad. Sci. USA 85,
 * 5051-5055], the present model, which contains only two variable,
 * predicts the occurrence of periodic Ca2+ spikes in the absence
 * of InsP3 oscillations. Such results indicate that repetitive
 * Ca2+ spikes evoked by external stimuli do not necessarily require
 * the concomitant, periodic variation of InsP3. The model is closely
 * related to that proposed by Kuba and Takeshita [Kuba, K. & Takeshita,
 * S. (1981) J. Theor. Biol. 93, 1009-1031] for Ca2+ oscillations
 * in sympathetic neurones, based on Ca2+-induced Ca2+ release.
 * We extend their results by showing the minimal conditions in
 * which the latter process gives rise to periodic behaviour and
 * take into account the role of the rise in InsP3 caused by external
 * stimulation. The analysis further shows how signal-induced Ca2+
 * oscillations might be effectively encoded in terms of their
 * frequency through the phospohorylation of a cellular substrate
 * by a protein kinase activated by cytosolic Ca2+.
 * 
 * Minimla model for signal-induced Ca2+ oscillations and for their
 * frequency encoding through protein phosphorylation, Albert Goldbeter,
 * Genevieve Dupont, Michael J. Berridge, 1990, Biophysics, 87,
 * 1461-1465. PubMed ID: 2304911
 * 
 * cell diagram
 * 
 * [[Image file: goldbeter_1990.png]]
 * 
 * Schematic representation of the mechanism generating Ca2+ oscillations,
 * based on the self-amplified release of Ca2+ from intracellular
 * stores.
 */

import nsrunit;
unit conversion on;
unit s=1 second^1;
unit per_s=1 second^(-1);
unit uM=1E-3 meter^(-3)*mole^1;
unit uM_per_s=1E-3 meter^(-3)*second^(-1)*mole^1;

math main {
	realDomain time s;
	time.min=0;
	extern time.max;
	extern time.delta;
	real Z(time) uM;
	when(time=time.min) Z=0.5;
	real Y(time) uM;
	when(time=time.min) Y=1.75;
	real v_0 uM_per_s;
	v_0=1;
	real v_1 uM_per_s;
	v_1=7.3;
	real v_2(time) uM_per_s;
	real v_3(time) uM_per_s;
	real k per_s;
	k=10;
	real k_f per_s;
	k_f=1;
	real beta dimensionless;
	beta=0.301;
	real V_M2 uM_per_s;
	V_M2=65;
	real v_2.K_2 uM;
	v_2.K_2=1;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real n dimensionless;
//	n=2;
	real V_M3 uM_per_s;
	V_M3=500;
	real K_R uM;
	K_R=2;
	real K_A uM;
	K_A=0.9;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real m dimensionless;
//	m=2;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real p dimensionless;
//	p=4;
	real W_star(time) dimensionless;
	when(time=time.min) W_star=0;
	real W_T uM;
	W_T=1;
	real v_P uM_per_s;
	v_P=5;
	real v_K(time) uM_per_s;
	real K_1 dimensionless;
	K_1=0.1;
	real W_star.K_2 dimensionless;
	W_star.K_2=0.1;
	real V_MK uM_per_s;
	V_MK=40;
	real K_a uM;
	K_a=2.5;

	// <component name="environment">

	// <component name="Ca">
	Z:time=(v_0+v_1*beta+(-1)*v_2+v_3+k_f*Y+(-1)*(k*Z));
	Y:time=(v_2+(-1)*v_3+(-1)*(k_f*Y));

	// <component name="v_2">
	v_2=(V_M2*(Z^2/(v_2.K_2^2+Z^2)));

	// <component name="v_3">
	v_3=(V_M3*(Y^2/(K_R^2+Y^2))*(Z^4/(K_A^4+Z^4)));

	// <component name="W_star">
	W_star:time=(v_P/W_T*(v_K/v_P*((1-W_star)/(K_1+1+(-1)*W_star))-W_star/(W_star.K_2+W_star)));

	// <component name="v_K">
	v_K=(V_MK*(Z/(K_a+Z)));
}