/*
 * Hilemann-Noble Atrial Cell Model 1987
 * 
 * Model Status
 * 
 * This model is known to run in OpenCell and COR and reproduce
 * the published results.
 * 
 * Model Structure
 * 
 * ABSTRACT: Interactions of electrogenic sodium-calcium exchange,
 * calcium channel and sarcoplasmic reticulum in the mammalian
 * heart have been explored by simulation of extracellular calcium
 * transients measured with tetramethylmurexide in rabbit atrium.
 * The approach has been to use the simplest possible formulations
 * of these mechanisms, which together with a minimum number of
 * additional mechanisms allow reconstruction of action potentials,
 * intracellular calcium transients and extracellular calcium transients.
 * A 3:1 sodium-calcium exchange stoichiometry is assumed. Calcium-channel
 * inactivation is assumed to take place by a voltage-dependent
 * mechanism, which is accelerated by a rise in intracellular calcium;
 * intracellular calcium release becomes a major physiological
 * regulator of calcium influx via calcium channels. A calcium
 * release mechanism is assumed, which is both calcium- and voltage-sensitive,
 * and which undergoes prolonged inactivation. 200 microM cytosolic
 * calcium buffer is assumed. For most simulations only instantaneous
 * potassium conductances are simulated so as to study the other
 * mechanisms independently of time- and calcium-dependent outward
 * current. Thus, the model reconstructs extracellular calcium
 * transients and typical action-potential configuration changes
 * during steady-state and non-steady-state stimulation from the
 * mechanisms directly involved in trans-sarcolemmal calcium movements.
 * The model predicts relatively small trans-sarcolemmal calcium
 * movements during regular stimulation (ca. 2 mumol kg-1 fresh
 * mass per excitation); calcium current is fully activated within
 * 2 ms of excitation, inactivation is substantially complete within
 * 30 ms, and sodium-calcium exchange significantly resists repolarization
 * from approximately -30 mV. Net calcium movements many times
 * larger are possible during non-steady-state stimulation. Long
 * action potentials at premature excitations or after inhibition
 * of calcium release can be supported almost exclusively by calcium
 * current (net calcium influx 5-30 mumol kg-1 fresh mass); action
 * potentials during potentiated post-stimulatory contractions
 * can be supported almost exclusively by sodium-calcium exchange
 * (net calcium efflux 4-20 mumol kg-1 fresh mass). Large calcium
 * movements between the extracellular space and the sarcoplasmic
 * reticulum can take place through the cytosol with virtually
 * no contractile activation. The simulations provide integrated
 * explanations of electrical activity, contractile function and
 * trans-sarcolemmal calcium movements, which were outside the
 * explanatory range of previous models.
 * 
 * The original paper reference is cited below:
 * 
 * Excitation-contraction coupling and extracellular calcium transients
 * in rabbit atrium: reconstruction of the basic cellular mechanisms,
 * Hilemann, D.W. and Noble, D. 1987,Proc. R. Soc. Lond., B230,
 * 163-205. PubMed ID: 2884668
 * 
 * cell diagram of the Hilemann-Noble model showing ionic currents
 * across the sarcolemma
 * 
 * [[Image file: hilgemann_noble_1987.png]]
 * 
 * A schematic diagram describing the current flows across the
 * cell membrane that are captured in the Hilemann-Noble model.
 */

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=-88;
	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=0.006;
	real i_K1(time) nanoA;
	real i_b_Na(time) nanoA;
	real i_b_Ca(time) nanoA;
	real i_b_K(time) nanoA;
	real i_NaK(time) nanoA;
	real i_NaCa(time) nanoA;
	real i_Na(time) nanoA;
	real i_si(time) nanoA;
	real i_Stim(time) nanoA;
	real stim_start second;
	stim_start=0.1;
	real stim_end second;
	stim_end=10000;
	real stim_period second;
	stim_period=1;
	real stim_duration second;
	stim_duration=0.002;
	real stim_amplitude nanoA;
	stim_amplitude=-200;
	real g_Na microS;
	g_Na=50;
	real E_mh(time) millivolt;
	real Na_o millimolar;
	Na_o=140;
	real Na_i(time) millimolar;
	when(time=time.min) Na_i=6.5;
	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 i_NaK_max nanoA;
	i_NaK_max=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.012;
	real E_Ca(time) millivolt;
	real g_b_Ca microS;
	g_b_Ca=0.005;
	real Ca_o(time) millimolar;
	when(time=time.min) Ca_o=2;
	real Ca_i(time) millimolar;
	when(time=time.min) Ca_i=1e-5;
	real k_NaCa nanoA;
	k_NaCa=0.01;
//	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 E_K(time) millivolt;
	real g_b_K microS;
	g_b_K=0.17;
	real g_K1 microS;
	g_K1=1.7;
	real K_m_K1 millimolar;
	K_m_K1=10;
	real i_siCa(time) nanoA;
	real i_siK(time) nanoA;
	real i_siNa(time) nanoA;
	real P_si nanoA_per_millimolar;
	P_si=5;
	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 delta_d millivolt;
	delta_d=0.0001;
	real E0_d(time) millivolt;
	real alpha_f_Ca(time) per_second;
	real beta_f_Ca(time) per_second;
	real CaChoff(time) dimensionless;
	real delta_f millivolt;
	delta_f=0.0001;
	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 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 Cab millimolar;
	Cab=2;
	real K_diff per_second;
	K_diff=0.0005;
	real Ve 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.08;
	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;
	real dCaCalmoddt(time) millimolar_per_second;
	real dCaTropdt(time) millimolar_per_second;

	// <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)*(i_Stim+i_K1+i_b_Na+i_b_Ca+i_b_K+i_NaK+i_NaCa+i_Na+i_si)/C_m);

	// <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="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">
	E_K=(RTONF*ln(K_c/K_i));
	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="second_inward_calcium_current">
	i_si=(i_siCa+i_siK+i_siNa);
	i_siCa=(4*P_si*d*CaChon*(V-(50 millivolt))/RTONF/(1-exp((-1)*1*(V-(50 millivolt))*2/RTONF))*(Ca_i*exp((100 millivolt)/RTONF)-Ca_o*exp((-1)*2*(V-(50 millivolt))/RTONF)));
	i_siK=(.002*P_si*d*CaChon*(V-(50 millivolt))/RTONF/(1-exp((-1)*1*(V-(50 millivolt))/RTONF))*(K_i*exp((50 millivolt)/RTONF)-K_c*exp((-1)*1*(V-(50 millivolt))/RTONF)));
	i_siNa=(.01*P_si*d*CaChon*(V-(50 millivolt))/RTONF/(1-exp((-1)*1*(V-(50 millivolt))/RTONF))*(Na_i*exp((50 millivolt)/RTONF)-Na_o*exp((-1)*1*(V-(50 millivolt))/RTONF)));

	// <component name="second_inward_calcium_current_d_gate">
	E0_d=(V+(24 millivolt)-(5 millivolt));
	alpha_d=(if (abs(E0_d)<delta_d) (120 per_second) else (30 per_millivolt_second)*E0_d/(1-exp((-1)*1*E0_d/(4 millivolt))));
	beta_d=(if (abs(E0_d)<delta_d) (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="second_inward_calcium_current_f_Ca_gate">
	E0_f=(V+(34 millivolt));
	alpha_f_Ca=(if (abs(E0_f)<delta_f) (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-(.6 per_second)*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+i_NaK*3+i_NaCa*3+i_siNa));

	// <component name="extracellular_calcium_concentration">
	Ca_o:time=((Cab-Ca_o)*K_diff-(1 mA_nA)*(i_siCa+i_NaCa+i_b_Ca)/(2*(1 litre_micrometre3)*Ve*F));

	// <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_siK+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);
	Ve=(V_Cell*V_e_ratio);
	Ca_i:time=((-1)*(1 mA_nA)/(2*(1 litre_micrometre3)*V_i*F)*(i_siCa+i_b_Ca-2*i_NaCa)+i_rel*V_rel_ratio/V_i_ratio-dCaCalmoddt-dCaTropdt-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);
	dCaCalmoddt=(alpha_Calmod*Ca_i*(Calmod-Ca_Calmod)-beta_Calmod*Ca_Calmod);
	dCaTropdt=(alpha_Trop*Ca_i*(Trop-Ca_Trop)-beta_Trop*Ca_Trop);
}