/*
 * Model Status
 * 
 * This CellML model runs in OpenCell and COR and reproduce the
 * published results. The units have been checked and they are
 * consistent.
 * 
 * Model Structure
 * 
 * ABSTRACT: The hypothesis that calcium release from the sarcoplasmic
 * reticulum in cardiac muscle is induced by rises in free cytosolic
 * calcium (Fabiato 1983, Am. J. Physiol 245) allows the possibility
 * that the release could be at least partly regenerative. There
 * would then be a non-linear relation between calcium current
 * and calcium release. We have investigated this possibility in
 * a single-cell version of the rabbit-atrial model developed by
 * Hilgemann & Noble (1987, Proc. R. Soc. Lond. B 230). The model
 * predicts different voltage ranges of activation for calcium-dependent
 * processes (like the sodium-calcium exchange current, contraction
 * or Fura-2 signals) and the calcium current, in agreement with
 * the experimental results obtained by Earm et al. (1990, Proc.
 * R. Soc. Lond. B 240) on exchange current tails, Cannell et al.
 * (1987, Science, Wash. 238) by using Fura-2 signals, and Fedida
 * et al. (1987, J. Physiol., Lond. 385) and Talo et al. (1988,
 * Biology of isolated adult cardiac myocytes) by using contraction.
 * However, when the Fura-2 concentration is sufficiently high
 * (greater than 200 microM) the activation ranges become very
 * similar as the buffering properties of Fura-2 are sufficient
 * to remove the regenerative effect. It is therefore important
 * to allow for the buffering properties of calcium indicators
 * when investigating the correlation between calcium current and
 * calcium release.
 * 
 * schematic diagram of a human atrial myocyte
 * 
 * [[Image file: earm_1990.png]]
 * 
 * A schematic representation of the Earm and Noble 1990 model.
 * 
 * The original paper reference is cited below:
 * 
 * A model of the single atrial cell: relation between calcium
 * current and calcium release, Y. E. Earm and D. Noble, 1990,
 * Proceedings of the Royal Society of London. Series B, Biological
 * Sciences , 240, 83-96. PubMed ID: 1972993
 */

import nsrunit;
unit conversion on;
unit per_second=1 second^(-1);
// unit millivolt predefined
unit per_millivolt=1E3 kilogram^(-1)*meter^(-2)*second^3*ampere^1;
unit per_millivolt_second=1E3 kilogram^(-1)*meter^(-2)*second^2*ampere^1;
unit microS=1E-6 kilogram^(-1)*meter^(-2)*second^3*ampere^2;
unit microF=1E-6 kilogram^(-1)*meter^(-2)*second^4*ampere^2;
unit nanoA=1E-9 ampere^1;
unit mA_nA=1E6  dimensionless;
// unit millimolar predefined
unit millimolar4=1 meter^(-12)*mole^4;
unit nanoA_per_millimolar=1E-9 meter^3*ampere^1*mole^(-1);
unit millimolar_per_second=1 meter^(-3)*second^(-1)*mole^1;
unit per_millimolar_second=1 meter^3*second^(-1)*mole^(-1);
// unit micrometre predefined
unit joule_per_kilomole_kelvin=.001 kilogram^1*meter^2*second^(-2)*kelvin^(-1)*mole^(-1);
unit coulomb_per_mole=1 second^1*ampere^1*mole^(-1);
unit micrometre3=1E-18 meter^3;
unit litre_micrometre3=1E15  dimensionless;

math main {
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real V(time) millivolt;
	when(time=time.min) V=-91.6;
	real R joule_per_kilomole_kelvin;
	R=8314.472;
	real T kelvin;
	T=310;
	real F coulomb_per_mole;
	F=96485.3415;
	real RTONF millivolt;
	real C_m microF;
	C_m=4e-5;
	real i_b_K(time) nanoA;
	real i_K1(time) nanoA;
	real i_to(time) nanoA;
	real i_b_Na(time) nanoA;
	real i_b_Ca(time) nanoA;
	real i_NaK(time) nanoA;
	real i_NaCa(time) nanoA;
	real i_Na(time) nanoA;
	real i_Ca_L(time) nanoA;
	real i_Stim(time) nanoA;
	real stim_start second;
	stim_start=0.1;
	real stim_end second;
	stim_end=100000;
	real stim_period second;
	stim_period=1;
	real stim_duration second;
	stim_duration=0.002;
	real stim_amplitude nanoA;
	stim_amplitude=-1.3;
	real g_Na microS;
	g_Na=0.5;
	real E_mh(time) millivolt;
	real Na_o millimolar;
	Na_o=140;
	real Na_i(time) millimolar;
	when(time=time.min) Na_i=6.48;
	real K_c millimolar;
	K_c=4;
	real K_i(time) millimolar;
	when(time=time.min) K_i=140;
	real m(time) dimensionless;
	when(time=time.min) m=0.076;
	real h(time) dimensionless;
	when(time=time.min) h=0.015;
	real alpha_m(time) per_second;
	real beta_m(time) per_second;
	real delta_m millivolt;
	delta_m=1e-5;
	real E0_m(time) millivolt;
	real alpha_h(time) per_second;
	real beta_h(time) per_second;
	real g_to microS;
	g_to=0.01;
	real E_K(time) millivolt;
	real g_to_s dimensionless;
	g_to_s=0;
	real r(time) dimensionless;
	when(time=time.min) r=0;
	real s(time) dimensionless;
	when(time=time.min) s=1;
	real alpha_s(time) per_second;
	real beta_s(time) per_second;
	real i_NaK_max nanoA;
	i_NaK_max=0.14;
	real K_mK millimolar;
	K_mK=1;
	real K_mNa millimolar;
	K_mNa=40;
	real E_Na(time) millivolt;
	real g_b_Na microS;
	g_b_Na=0.00012;
	real E_Ca(time) millivolt;
	real g_b_Ca microS;
	g_b_Ca=5e-5;
	real Ca_o millimolar;
	Ca_o=2;
	real Ca_i(time) millimolar;
	when(time=time.min) Ca_i=1e-5;
	real k_NaCa nanoA;
	k_NaCa=0.0001;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real n_NaCa dimensionless;
//	n_NaCa=3;
	real d_NaCa dimensionless;
	d_NaCa=0.0001;
	real gamma dimensionless;
	gamma=0.5;
	real g_b_K microS;
	g_b_K=0.0017;
	real g_K1 microS;
	g_K1=0.017;
	real K_m_K1 millimolar;
	K_m_K1=10;
	real i_Ca_L_Ca(time) nanoA;
	real i_Ca_L_K(time) nanoA;
	real i_Ca_L_Na(time) nanoA;
	real P_Ca_L nanoA_per_millimolar;
	P_Ca_L=0.05;
	real d(time) dimensionless;
	when(time=time.min) d=0.0011;
	real f_Ca(time) dimensionless;
	when(time=time.min) f_Ca=0.785;
	real CaChon(time) dimensionless;
	real alpha_d(time) per_second;
	real beta_d(time) per_second;
	real E0_d(time) millivolt;
	real alpha_f_Ca(time) per_second;
	real beta_f_Ca(time) per_second;
	real CaChoff(time) dimensionless;
	real E0_f(time) millivolt;
	real i_up(time) millimolar_per_second;
	real K_1 dimensionless;
	real K_2(time) millimolar;
	real K_cyca millimolar;
	K_cyca=0.0003;
	real K_xcs dimensionless;
	K_xcs=0.4;
	real K_srca millimolar;
	K_srca=0.5;
	real alpha_up millimolar_per_second;
	alpha_up=3;
	real beta_up millimolar_per_second;
	beta_up=0.23;
	real Ca_up(time) millimolar;
	when(time=time.min) Ca_up=0.3;
	real i_rel(time) millimolar_per_second;
	real VoltDep(time) dimensionless;
	real RegBindSite(time) dimensionless;
	real ActRate(time) per_second;
	real InactRate(time) per_second;
	real K_leak_rate per_second;
	K_leak_rate=0;
	real K_m_rel per_second;
	K_m_rel=250;
	real PrecFrac(time) dimensionless;
	real ActFrac(time) dimensionless;
	when(time=time.min) ActFrac=0;
	real ProdFrac(time) dimensionless;
	when(time=time.min) ProdFrac=0;
	real ProdFracRate per_second;
	ProdFracRate=1;
	real Ca_rel(time) millimolar;
	when(time=time.min) Ca_rel=0.3;
	real i_trans(time) millimolar_per_second;
	real alpha_tr per_second;
	alpha_tr=50;
	real V_i micrometre3;
	real Ca_Calmod(time) millimolar;
	when(time=time.min) Ca_Calmod=0.0005;
	real Ca_Trop(time) millimolar;
	when(time=time.min) Ca_Trop=0.0015;
	real Calmod millimolar;
	Calmod=0.02;
	real Trop millimolar;
	Trop=0.15;
	real alpha_Calmod per_millimolar_second;
	alpha_Calmod=100000;
	real beta_Calmod per_second;
	beta_Calmod=50;
	real alpha_Trop per_millimolar_second;
	alpha_Trop=100000;
	real beta_Trop per_second;
	beta_Trop=200;
	real radius micrometre;
	radius=0.01;
	real length micrometre;
	length=0.08;
	real V_Cell micrometre3;
	real V_i_ratio dimensionless;
	real V_rel_ratio dimensionless;
	V_rel_ratio=0.1;
	real V_e_ratio dimensionless;
	V_e_ratio=0.4;
	real V_up_ratio dimensionless;
	V_up_ratio=0.01;

	// <component name="environment">

	// <component name="membrane">
	RTONF=(R*T/F);
	i_Stim=(if (((time>=stim_start) and (time<=stim_end)) and ((time-stim_start-floor((time-stim_start)/stim_period)*stim_period)<=stim_duration)) stim_amplitude else (0 nanoA));
	V:time=((-1)*1/C_m*(i_Stim+i_b_K+i_K1+i_to+i_b_Na+i_b_Ca+i_NaK+i_NaCa+i_Na+i_Ca_L));

	// <component name="fast_sodium_current">
	E_mh=(RTONF*ln((Na_o+.12*K_c)/(Na_i+.12*K_i)));
	i_Na=(g_Na*m^3*h*(V-E_mh));

	// <component name="fast_sodium_current_m_gate">
	E0_m=(V+(41 millivolt));
	alpha_m=(if (abs(E0_m)<delta_m) (2E3 per_second) else (200 per_millivolt_second)*E0_m/(1-exp((-1)*(.1 per_millivolt)*E0_m)));
	beta_m=((8E3 per_second)*exp((-1)*(.056 per_millivolt)*(V+(66 millivolt))));
	m:time=(alpha_m*(1-m)-beta_m*m);

	// <component name="fast_sodium_current_h_gate">
	alpha_h=((20 per_second)*exp((-1)*(.125 per_millivolt)*(V+(75 millivolt))));
	beta_h=((2E3 per_second)/(1+320*exp((-1)*(.1 per_millivolt)*(V+(75 millivolt)))));
	h:time=(alpha_h*(1-h)-beta_h*h);

	// <component name="transient_outward_current">
	E_K=(RTONF*ln(K_c/K_i));
	i_to=(g_to*(g_to_s+s*(1-g_to_s))*r*(V-E_K));

	// <component name="transient_outward_current_r_gate">
	r:time=((333 per_second)*(1/(1+exp((-1)*(V+(4 millivolt))/(5 millivolt)))-r));

	// <component name="transient_outward_current_s_gate">
	alpha_s=((.033 per_second)*exp((-1)*V/(17 millivolt)));
	beta_s=((33 per_second)/(1+exp((-1)*(.125 per_millivolt)*(V+(10 millivolt)))));
	s:time=(alpha_s*(1-s)-beta_s*s);

	// <component name="sodium_potassium_pump">
	i_NaK=(i_NaK_max*K_c/(K_mK+K_c)*Na_i/(K_mNa+Na_i));

	// <component name="sodium_background_current">
	E_Na=(RTONF*ln(Na_o/Na_i));
	i_b_Na=(g_b_Na*(V-E_Na));

	// <component name="calcium_background_current">
	E_Ca=(.5*RTONF*ln(Ca_o/Ca_i));
	i_b_Ca=(g_b_Ca*(V-E_Ca));

	// <component name="Na_Ca_exchanger">
	i_NaCa=(k_NaCa*(exp(gamma*(3-2)*V/RTONF)*Na_i^3*Ca_o-exp((gamma-1)*(3-2)*V/RTONF)*Na_o^3*Ca_i)/(((1 millimolar4)+d_NaCa*(Ca_i*Na_o^3+Ca_o*Na_i^3))*(1+Ca_i/(.0069 millimolar))));

	// <component name="potassium_background_current">
	i_b_K=(g_b_K*(V-E_K));

	// <component name="time_independent_potassium_current">
	i_K1=(g_K1*K_c/(K_c+K_m_K1)*(V-E_K)/(1+exp((V-E_K-(10 millivolt))*2/RTONF)));

	// <component name="L_type_calcium_current">
	i_Ca_L=(i_Ca_L_Ca+i_Ca_L_K+i_Ca_L_Na);
	i_Ca_L_Ca=(4*P_Ca_L*d*CaChon*(V-(50 millivolt))/RTONF/(1-exp((-1)*(V-(50 millivolt))*2/RTONF))*(Ca_i*exp((100 millivolt)/RTONF)-Ca_o*exp((-1)*(V-(50 millivolt))*2/RTONF)));
	i_Ca_L_K=(.002*P_Ca_L*d*CaChon*(V-(50 millivolt))/RTONF/(1-exp((-1)*(V-(50 millivolt))/RTONF))*(K_i*exp((50 millivolt)/RTONF)-K_c*exp((-1)*(V-(50 millivolt))/RTONF)));
	i_Ca_L_Na=(.01*P_Ca_L*d*CaChon*(V-(50 millivolt))/RTONF/(1-exp((-1)*(V-(50 millivolt))/RTONF))*(Na_i*exp((50 millivolt)/RTONF)-Na_o*exp((-1)*(V-(50 millivolt))/RTONF)));

	// <component name="L_type_calcium_current_d_gate">
	E0_d=(V+(19 millivolt));
	alpha_d=(if (abs(E0_d)<(1E-4 millivolt)) (120 per_second) else (30 per_millivolt_second)*E0_d/(1-exp((-1)*E0_d/(4 millivolt))));
	beta_d=(if (abs(E0_d)<(1E-4 millivolt)) (120 per_second) else (12 per_millivolt_second)*E0_d/(exp(E0_d/(10 millivolt))-1));
	d:time=(alpha_d*(1-d)-beta_d*d);

	// <component name="L_type_calcium_current_f_Ca_gate">
	E0_f=(V+(34 millivolt));
	alpha_f_Ca=(if (abs(E0_f)<(1E-4 millivolt)) (25 per_second) else (6.25 per_millivolt_second)*E0_f/(exp(E0_f/(4 millivolt))-1));
	beta_f_Ca=((12 per_second)/(1+exp((-1)*E0_f/(4 millivolt))));
	f_Ca:time=((120*(1-f_Ca)*CaChoff+(1-f_Ca)*(1-CaChoff))*beta_f_Ca-alpha_f_Ca*f_Ca);
	CaChoff=(Ca_i/((.001 millimolar)+Ca_i));
	CaChon=((1-f_Ca)*(1-CaChoff));

	// <component name="sarcoplasmic_reticulum_calcium_pump">
	K_1=(K_cyca*K_xcs/K_srca);
	K_2=(Ca_i+Ca_up*K_1+K_cyca*K_xcs+K_cyca);
	i_up=(Ca_i/K_2*alpha_up-Ca_up*K_1/K_2*beta_up);

	// <component name="calcium_release">
	PrecFrac=(1-ActFrac-ProdFrac);
	VoltDep=exp((.08 per_millivolt)*(V-(40 millivolt)));
	RegBindSite=((Ca_i/(Ca_i+(5E-4 millimolar)))^2);
	ActRate=((600 per_second)*VoltDep+(500 per_second)*RegBindSite);
	InactRate=((60 per_second)+(500 per_second)*RegBindSite);
	ActFrac:time=(PrecFrac*ActRate-ActFrac*InactRate);
	ProdFrac:time=(ActFrac*InactRate-ProdFracRate*ProdFrac);
	i_rel=(((ActFrac/(ActFrac+.25))^2*K_m_rel+K_leak_rate)*Ca_rel);

	// <component name="calcium_translocation">
	i_trans=((Ca_up-Ca_rel)*alpha_tr);

	// <component name="extracellular_sodium_concentration">

	// <component name="intracellular_sodium_concentration">
	Na_i:time=((-1)*(1 mA_nA)/((1 litre_micrometre3)*V_i*F)*(i_Na+i_b_Na+3*i_NaK+3*i_NaCa+i_Ca_L_Na));

	// <component name="extracellular_calcium_concentration">

	// <component name="extracellular_potassium_concentration">

	// <component name="intracellular_potassium_concentration">
	K_i:time=((-1)*(1 mA_nA)/((1 litre_micrometre3)*V_i*F)*(i_K1+i_Ca_L_K+i_to+i_b_K-2*i_NaK));

	// <component name="intracellular_calcium_concentration">
	V_Cell=(3.141592654*radius^2*length);
	V_i_ratio=(1-V_e_ratio-V_up_ratio-V_rel_ratio);
	V_i=(V_Cell*V_i_ratio);
	Ca_i:time=((-1)*(1 mA_nA)/(2*(1 litre_micrometre3)*V_i*F)*(i_Ca_L_Ca+i_b_Ca-2*i_NaCa)+i_rel*V_rel_ratio/V_i_ratio-Ca_Calmod:time-Ca_Trop:time-i_up);
	Ca_up:time=(V_i_ratio/V_up_ratio*i_up-i_trans);
	Ca_rel:time=(V_up_ratio/V_rel_ratio*i_trans-i_rel);
	Ca_Calmod:time=(alpha_Calmod*Ca_i*(Calmod-Ca_Calmod)-beta_Calmod*Ca_Calmod);
	Ca_Trop:time=(alpha_Trop*Ca_i*(Trop-Ca_Trop)-beta_Trop*Ca_Trop);
}