/*
 * Modelling the Interellations Between Calcium Oscillations and
 * ER Membrane Potential Oscillations
 * 
 * Model Status
 * 
 * This is the original unchecked version of the model imported
 * from the previous CellML model repository, 24-Jan-2006.
 * 
 * Model Structure
 * 
 * Oscillations in the cytosolic calcium concentration ([Ca2+]i),
 * in the form of repetitive spikes, may be stimulated by hormones
 * or neurotransmitters in a variety of cell types, including hepatocytes,
 * oocytes, cardiomyocytes and fibroblasts. It is generally accepted
 * that the mechanism underlying these calcium oscillations is
 * related to the calcium-induced calcium release (CICR) mechanism.
 * Often, secondary messengers such as inositol 1,4,5-trisphosphate
 * (IP3) are also involved. Calcium oscillations play a role in
 * the regulation of cell excretion, muscle contraction, activation
 * of mammalian oocytes, call signalling and other cellular activities.
 * 
 * Calcium oscillations and their underlying molecular mechanisms
 * have been the subject of interest for several mathematical models,
 * including:
 * 
 * Li and Rinzel, IP3-Mediated Ca2+ Oscillations - A Simplified
 * Model, 1994;
 * 
 * D. Friel, [Ca2+]i Oscillations in in Sympathetic Neurons, 1995;
 * 
 * Chay et al., Intracellular Calcium Spikes in Non-Excitable Cells,
 * 1995; and
 * 
 * Keizer and Levine, RyR Adaptation and Ca2+ Oscillations, 1996.
 * 
 * In the Marhl et al. 1997 publication described here, the authors
 * present an electrochemical model which accounts for intracellular
 * calcium oscillations and their interrelations with oscillations
 * of the potential difference across the membrane of the endoplasmic
 * reticulum (ER) or other inracellular calium stores (see below).
 * The model shows tht when a calcium buffering system such as
 * calmodulin is included, calcium oscillations can arise without
 * a permanent influx of calium into the cell.
 * 
 * The model has been described here in CellML (the raw CellML
 * description of the Marhl et al. 1997 model can be downloaded
 * in various formats as described in ).
 * 
 * The complete original paper reference is cited below:
 * 
 * Modelling the interellations between calcium oscillations and
 * ER membrane potential oscillations, Marko Marhl, Stefan Schuster,
 * Milan Brumen, and Reinhart Heinrich, 1997, Biophysical Chemistry,
 * 63, 221-239. PubMed ID: 12362939
 * 
 * cell diagram
 * 
 * [[Image file: marhl_1997.png]]
 * 
 * Schematic diagram of the model system.
 */

import nsrunit;
// Warning: unit conversion turned off due to unit errors in 4 equation(s)
unit conversion off;
// unit micromolar predefined
// unit millivolt 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);
unit cm2=1E-4 meter^2;
unit micromolar_per_volt_second=1E-3 kilogram^(-1)*meter^(-5)*second^2*ampere^1*mole^1;
unit microS_per_cm2=.01 kilogram^(-1)*meter^(-4)*second^3*ampere^2;
unit joule_per_kelvin_mole=1 kilogram^1*meter^2*second^(-2)*kelvin^(-1)*mole^(-1);
unit coulomb_per_millimole=.001 second^(-1)*ampere^(-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:  Pr, Ca_cyt
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real CaPr(time) micromolar;
	real Pr(time) micromolar;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) Pr=0;
	real Pr_tot micromolar;
	Pr_tot=600.0;
	real k_plus second_order_rate_constant;
	k_plus=0.1;
	real k_minus first_order_rate_constant;
	k_minus=0.5;
	real Ca_cyt(time) micromolar;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) Ca_cyt=0;
	real J_ch(time) flux;
	real J_leak(time) flux;
	real J_pump(time) flux;
	real Ca_ER(time) micromolar;
	real Ca_tot micromolar;
	Ca_tot=45.0;
	real rho dimensionless;
	rho=0.01;
	real delta_psi(time) millivolt;
	real a(time) micromolar;
	real b(time) micromolar;
	real c(time) micromolar;
	real C_tot micromolar;
	C_tot=5.0E3;
	real A_tot micromolar;
	A_tot=3.89E3;
	real R joule_per_kelvin_mole;
	R=8.3143;
	real T kelvin;
	T=310.0;
	real F coulomb_per_millimole;
	F=96.4867;
	real E_Ca(time) millivolt;
	real S cm2;
	S=6.16E-3;
	real g_Ca(time) microS_per_cm2;
	real g_Ca_max microS_per_cm2;
	g_Ca_max=100.0;
	real k_ch(time) micromolar_per_volt_second;
	real K_Ca micromolar;
	K_Ca=5.0;
	real Vcyt litre;
	Vcyt=5.84E-11;
	real k_pump first_order_rate_constant;
	k_pump=76.0;
	real k_leak micromolar_per_volt_second;
	k_leak=10.0;

	// <component name="environment">

	// <component name="free_protein_calcium_binding_sites">
	Pr:time=(k_minus*CaPr-k_plus*Ca_cyt*Pr);
	CaPr=(Pr_tot-Pr);

	// <component name="cytosolic_calcium">
	Ca_cyt:time=(J_ch+J_leak+k_minus*CaPr-(J_pump+k_plus*Ca_cyt*Pr));

	// <component name="ER_calcium">
	Ca_ER=((Ca_tot-(Ca_cyt+(Pr_tot-Pr)))/rho);

	// <component name="potential_difference">
	delta_psi=(R*T/F*ln(((-1)*b-(b^2-4*a*c)^.5)/(2*a)));
	a=(rho*(2*Ca_cyt-(2*Pr+A_tot)));
	b=(C_tot-A_tot+2*(Ca_cyt-Pr)*(1+rho^2));
	c=(rho*(2*Ca_cyt-2*Pr+C_tot));

	// <component name="calcium_reversal_potential">
	E_Ca=(R*T/(2*F)*ln(Ca_ER/Ca_cyt));

	// <component name="flux_through_the_channels">
	J_ch=(k_ch*(E_Ca-delta_psi));
	k_ch=(g_Ca/(2*F*Vcyt));
	g_Ca=(g_Ca_max*S*(Ca_cyt/K_Ca)^2/(1+(Ca_cyt/K_Ca)^2));

	// <component name="calcium_pump">
	J_pump=(k_pump*Ca_cyt);

	// <component name="leak_flux">
	J_leak=(k_leak*(E_Ca-delta_psi));

	// <component name="model_parameters">
}