/*
 * Complex Calcium Oscillations and the Role of Mitochondria and
 * Cytosolic Proteins
 * 
 * Model Status
 * 
 * This CellML model runs in both OpenCell and COR to recreate
 * the published results. The units have been checked and are consistent.
 * 
 * Model Structure
 * 
 * ABSTRACT: Intracellular calcium oscillations, which are oscillatory
 * changes of cytosolic calcium concentration in response to agonist
 * stimulation, are experimentally well observed in various living
 * cells. Simple calcium oscillations represent the most common
 * pattern and many mathematical models have been published to
 * describe this type of oscillation. On the other hand, relatively
 * few theoretical studies have been proposed to give an explanation
 * of complex intracellular calcium oscillations, such as bursting
 * and chaos. In this paper, we develop a new possible mechanism
 * for complex calcium oscillations based on the interplay between
 * three calcium stores in the cell: the endoplasmic reticulum
 * (ER), mitochondria and cytosolic proteins. The majority ( approximately
 * 80%) of calcium released from the ER is first very quickly sequestered
 * by mitochondria. Afterwards, a much slower release of calcium
 * from the mitochondria serves as the calcium supply for the intermediate
 * calcium exchanges between the ER and the cytosolic proteins
 * causing bursting calcium oscillations. Depending on the permeability
 * of the ER channels and on the kinetic properties of calcium
 * binding to the cytosolic proteins, different patterns of complex
 * calcium oscillations appear. With our model, we are able to
 * explain simple calcium oscillations, bursting and chaos. Chaos
 * is also observed for calcium oscillations in the bursting mode.
 * 
 * The original paper reference is cited below:
 * 
 * Complex calcium oscillations and the role of mitochondria and
 * cytosolic proteins, Marko Marhl, Thomas Haberichter, Milan Brumen
 * and Reinhart Heinrich, 2000, Biosystems, 57, 75-86. PubMed ID:
 * 11004387
 * 
 * cell schematic for the model
 * 
 * [[Image file: marhl_2000.png]]
 * 
 * Schematic diagram of the model sysytem.
 */

import nsrunit;
unit conversion on;
// 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 {
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real Ca_Pr(time) micromolar;
	real Ca_tot micromolar;
	Ca_tot=90.0;
	real rho_ER dimensionless;
	rho_ER=0.01;
	real beta_ER dimensionless;
	beta_ER=0.0025;
	real rho_m dimensionless;
	rho_m=0.01;
	real beta_m dimensionless;
	beta_m=0.0025;
	real Ca_cyt(time) micromolar;
	when(time=time.min) Ca_cyt=0.05;
	real Ca_ER(time) micromolar;
	when(time=time.min) Ca_ER=1.0;
	real Ca_m(time) micromolar;
	when(time=time.min) Ca_m=0.4;
	real Pr(time) micromolar;
	real Pr_tot micromolar;
	Pr_tot=120.0;
	real k_plus second_order_rate_constant;
	k_plus=0.1;
	real k_minus first_order_rate_constant;
	k_minus=0.01;
	real J_ch(time) flux;
	real J_leak(time) flux;
	real J_pump(time) flux;
	real J_out(time) flux;
	real J_in(time) flux;
	real k_pump first_order_rate_constant;
	k_pump=20.0;
	real k_ch first_order_rate_constant;
	k_ch=4100.0;
	real K1 micromolar;
	K1=5.0;
	real k_leak first_order_rate_constant;
	k_leak=0.05;
	real k_in flux;
	k_in=300.0;
	real K2 micromolar;
	K2=0.8;
	real k_out first_order_rate_constant;
	k_out=125.0;
	real k_m first_order_rate_constant;
	k_m=0.00625;
	real K3 micromolar;
	K3=5.0;

	// <component name="environment">

	// <component name="total_calcium">
	Ca_Pr=(Ca_tot-(Ca_cyt+rho_ER/beta_ER*Ca_ER+rho_m/beta_m*Ca_m));

	// <component name="total_protein">
	Pr=(Pr_tot-Ca_Pr);

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

	// <component name="ER_calcium">
	Ca_ER:time=(beta_ER/rho_ER*(J_pump-(J_ch+J_leak)));

	// <component name="mitochondrial_calcium">
	Ca_m:time=(beta_m/rho_m*(J_in-J_out));

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

	// <component name="Ca_efflux_from_the_ER">
	J_ch=(k_ch*(Ca_cyt^2/(K1^2+Ca_cyt^2))*(Ca_ER-Ca_cyt));

	// <component name="Ca_leak_flux_from_the_ER">
	J_leak=(k_leak*(Ca_ER-Ca_cyt));

	// <component name="mitochondrial_Ca_uptake">
	J_in=(k_in*(Ca_cyt^8/(K2^8+Ca_cyt^8)));

	// <component name="mitochondrial_Ca_release">
	J_out=((k_out*(Ca_cyt^2/(K3^2+Ca_cyt^2))+k_m)*Ca_m);
}