/*
 * Na/Ca Exchange in Models for Pancreatic Beta-Cells
 * 
 * Model Status
 * 
 * This CellML model represents Model III from the published paper.
 * The CellML matches the published equations and the model runs
 * in OpenCell and COR but it does not replicate the published
 * results.
 * 
 * Model Structure
 * 
 * ABSTRACT: In the presence of an insulinotropic glucose concentration,
 * beta-cells, in intact pancreatic islets, exhibit periodic bursting
 * electrical activity consisting of an alternation of active and
 * silent phases. The fraction of time spent in the active phase
 * over a period is called the plateau fraction and is correlated
 * with the rate of insulin release. However, the mechanisms that
 * regulate the plateau fraction remain unclear. In this paper
 * we investigate the possible role of the plasma membrane Na+/Ca2+
 * exchange of the beta-cell in controlling the plateau fraction.
 * We have extended different single-cell models to incorporate
 * this Ca2+-activated electrogenic Ca2+ transporter. We find that
 * the Na+/Ca2+ exchange can provide a physiological mechanism
 * to increase the plateau fraction as the glucose concentration
 * is raised. In addition, we show theoretically that the Na+/Ca2+
 * exchanger is a key regulator of the cytoplasmic calcium concentration
 * in clusters of heterogeneous cells with gap-junctional electrical
 * coupling.
 * 
 * The original paper reference is cited below:
 * 
 * Effect of Na/Ca Exchange on Plateau Fraction and [Ca2+]i in
 * Models for Bursting in Pancreatic Beta-Cells, David Gall and
 * Isabella Susa, 1999, Biophysical Journal, 77, 45-53. PubMed
 * ID: 10388739
 * 
 * cell schematic for the model
 * 
 * [[Image file: gall_1999.png]]
 * 
 * Schematic diagram of the pancreatic beta-cell plasma membrane
 * showing the ionic currents captured by the three single cell
 * mathematical models.
 */

import nsrunit;
unit conversion on;
// unit millisecond predefined
unit per_millisecond=1E3 second^(-1);
unit picoA=1E-12 ampere^1;
unit femtoF=1E-15 kilogram^(-1)*meter^(-2)*second^4*ampere^2;
// unit millimolar predefined
// unit micromolar predefined
// unit millivolt predefined
unit picoS=1E-12 kilogram^(-1)*meter^(-2)*second^3*ampere^2;
unit mole_per_microlitre_coulomb=1E9 meter^(-3)*second^(-1)*ampere^(-1)*mole^1;

math main {
	realDomain time millisecond;
	time.min=0;
	extern time.max;
	extern time.delta;
	real V(time) millivolt;
	when(time=time.min) V=-76.0;
	real Cm femtoF;
	Cm=5310.0;
	real i_Ca(time) picoA;
	real i_K(time) picoA;
	real i_K_Ca(time) picoA;
	real i_Na_Ca(time) picoA;
	real V_K millivolt;
	V_K=-75.0;
	real g_K picoS;
	g_K=2700.0;
	real n(time) dimensionless;
	when(time=time.min) n=0.1;
	real n_infinity(time) dimensionless;
	real lamda dimensionless;
	lamda=0.85;
	real tau_n(time) millisecond;
	real V_n millivolt;
	V_n=-15.0;
	real S_n millivolt;
	S_n=5.6;
	real a millivolt;
	a=65.0;
	real b millivolt;
	b=20.0;
	real c millisecond;
	c=6.0;
	real V_ millivolt;
	V_=-75.0;
	real V_Ca millivolt;
	V_Ca=25.0;
	real g_Ca picoS;
	g_Ca=1000.0;
	real m_infinity(time) dimensionless;
	real V_m millivolt;
	V_m=-20.0;
	real S_m millivolt;
	S_m=12.0;
	real g_K_Ca picoS;
	g_K_Ca=30000.0;
	real K_d micromolar;
	K_d=70.0;
	real Ca_i(time) micromolar;
	when(time=time.min) Ca_i=0.52;
	real g_Na_Ca picoS;
	g_Na_Ca=1000.0;
	real K_1_2 micromolar;
	K_1_2=1.5;
	real V_Na_Ca(time) millivolt;
	real RT_F millivolt;
	RT_F=26.54;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real nH dimensionless;
//	nH=5.0;
	real Ca_o micromolar;
	Ca_o=2600.0;
	real Na_i millimolar;
	Na_i=10.0;
	real Na_o millimolar;
	Na_o=140.0;
	real Ca_ret(time) micromolar;
	when(time=time.min) Ca_ret=0.7;
	real f dimensionless;
	f=0.001;
	real k_Ca per_millisecond;
	k_Ca=0.64;
	real k_rel per_millisecond;
	k_rel=0.0006;
	real k_pump per_millisecond;
	k_pump=0.2;
	real alpha mole_per_microlitre_coulomb;
	alpha=0.00006;

	// <component name="environment">

	// <component name="membrane">
	V:time=((-1)*(i_K+i_Ca+i_K_Ca+i_Na_Ca)/Cm);

	// <component name="rapidly_activating_K_current">
	i_K=(g_K*n*(V-V_K));

	// <component name="rapidly_activating_K_current_n_gate">
	n:time=(lamda*((n_infinity-n)/tau_n));
	n_infinity=(1/(1+exp((V_n-V)/S_n)));
	tau_n=(c/(exp((V-V_)/a)+exp((V_-V)/b)));

	// <component name="calcium_current">
	i_Ca=(g_Ca*m_infinity*(V-V_Ca));

	// <component name="calcium_current_m_gate">
	m_infinity=(1/(1+exp((V_m-V)/S_m)));

	// <component name="calcium_activated_K_current">
	i_K_Ca=(g_K_Ca*(Ca_i/(K_d+Ca_i))*(V-V_K));

	// <component name="Na_Ca_exchanger_current">
	i_Na_Ca=(g_Na_Ca*(Ca_i^5/(K_1_2^5+Ca_i^5))*(V-V_Na_Ca));
	V_Na_Ca=(RT_F*(3*ln(Na_o/Na_i-ln(Ca_o/Ca_i))));

	// <component name="ionic_concentrations">
	Ca_i:time=(f*((-1)*alpha*(i_Ca-2*i_Na_Ca)-k_Ca*Ca_i)+k_rel*(Ca_ret-Ca_i)-k_pump*Ca_i);
	Ca_ret:time=((-1)*k_rel*(Ca_ret-Ca_i)+k_pump*Ca_i);
}