/*
 * Di Francesco-Noble Purkinje Fibre Model 1985
 * 
 * Model Status
 * 
 * This model has been curated and unit-checked by Penny Noble
 * of Oxford University and is known to run in PCEnv and COR and
 * reproduce the results published in the paper it is based on.
 * 
 * Model Structure
 * 
 * During the years that followed the formulation of the McAllister-Noble-Tsien
 * Purkinje fibre model in 1975 and the Beeler-Reuter mammalian
 * ventricular model in 1977, many experiments were performed which
 * provided a greater insight into the working of the ion channels
 * in cardiac tissue. D. Di Francesco and D. Noble (1985) constructed
 * a new model of cardiac electrical activity which sought to incorporate
 * much of this new data (see the figure below).
 * 
 * The complete original paper reference is cited below:
 * 
 * A Model of the Cardiac Electrical Activity Incorporating Ionic
 * Pumps and Concentration Changes - Simulations of Ionic Currents
 * and Concentration Changes, Di Francesco, D. and Noble, D. Phil.
 * Trans. R. Soc. Lond. , B307, 353-398. PubMed ID: 2578676
 * 
 * cell diagram of the DFN model showing ionic currents, pumps
 * and exchangers within the sarcolemma and the sarcoplasmic reticulum
 * 
 * [[Image file: difrancesco_noble_1985.png]]
 * 
 * A schematic diagram describing the current flows across the
 * cell membrane that are captured in the DFN model.
 * 
 * the cellml rendering of the DFN model
 * 
 * [[Image file: cellml_rendering.gif]]
 * 
 * The network defined in the CellML description of the Di Francesco-Noble
 * model. A key describing the significance of the shapes of the
 * components and the colours of the connections between them is
 * in the notation guide. For simplicity, not all the variables
 * are shown.
 * 
 * The membrane physically contains the currents, exchangers and
 * pumps, as indicated by the blue arrows in . The currents act
 * independently and are not connected to each other. Several of
 * the channels encapsulate and contain further components which
 * represent activation and inactivation gates. The addition of
 * an encapsulation relationship informs modellers and processing
 * software that the gates are important parts of the current model.
 * It also prevents any other components that aren't also encapsulated
 * by the parent component from connecting to its gates, effectively
 * hiding them from the rest of the model.
 * 
 * The breakdown of the model into components and the definition
 * of encapsulation and containment relationships between them
 * is somewhat arbitrary. When considering how a model should be
 * broken into components, modellers are encouraged to consider
 * which parts of a model might be re-used and how the physiological
 * elements of the system being modelled are naturally bounded.
 * Containment relationships should be used to provide simple rendering
 * information for processing software (ideally, this will correspond
 * to the layout of the physical system), and encapsulation should
 * be used to group sets of components into sub-models.
 */

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 millivolt_second=.001 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 microS_per_millimolar=1E-6 kilogram^(-1)*meter^1*second^3*ampere^2*mole^(-1);
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 micrometre predefined
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=-87;
	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 microF;
	C=0.075;
	real i_pulse nanoA;
	i_pulse=0;
	real i_f(time) nanoA;
	real i_K(time) nanoA;
	real i_K1(time) nanoA;
	real i_to(time) nanoA;
	real i_Na_b(time) nanoA;
	real i_Ca_b(time) nanoA;
	real i_p(time) nanoA;
	real i_NaCa(time) nanoA;
	real i_Na(time) nanoA;
	real i_si(time) nanoA;
	real i_fNa(time) nanoA;
	real E_Na(time) millivolt;
	real E_K(time) millivolt;
	real i_fK(time) nanoA;
	real g_f_Na microS;
	g_f_Na=3;
	real g_f_K microS;
	g_f_K=3;
	real Km_f millimolar;
	Km_f=45;
	real Kc(time) millimolar;
	when(time=time.min) Kc=4;
	real Ki(time) millimolar;
	when(time=time.min) Ki=140;
	real Nai(time) millimolar;
	when(time=time.min) Nai=8;
	real Nao millimolar;
	Nao=140;
	real y(time) dimensionless;
	when(time=time.min) y=0.2;
	real alpha_y(time) per_second;
	real beta_y(time) per_second;
	real delta_y millivolt;
	delta_y=1e-5;
	real E0_y(time) millivolt;
	real I_K(time) nanoA;
	real i_K_max nanoA;
	i_K_max=180;
	real x(time) dimensionless;
	when(time=time.min) x=0.01;
	real alpha_x(time) per_second;
	real beta_x(time) per_second;
	real g_K1 microS;
	g_K1=920;
	real Km_K1 millimolar;
	Km_K1=210;
	real Km_to millimolar;
	Km_to=10;
	real Km_Ca millimolar;
	Km_Ca=0.0005;
	real g_to microS_per_millimolar;
	g_to=0.28;
	real Cai(time) millimolar;
	when(time=time.min) Cai=5e-5;
	real s(time) dimensionless;
	when(time=time.min) s=1;
	real alpha_s(time) per_second;
	real beta_s(time) per_second;
	real g_Nab microS;
	g_Nab=0.18;
	real E_Ca(time) millivolt;
	real g_Cab microS;
	g_Cab=0.02;
	real Cao millimolar;
	Cao=2;
	real I_p nanoA;
	I_p=125;
	real K_mK millimolar;
	K_mK=1;
	real K_mNa millimolar;
	K_mNa=40;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real n_NaCa dimensionless;
//	n_NaCa=3;
	real K_NaCa nanoA;
	K_NaCa=0.02;
	real d_NaCa dimensionless;
	d_NaCa=0.001;
	real gamma dimensionless;
	gamma=0.5;
	real g_Na microS;
	g_Na=750;
	real E_mh(time) millivolt;
	real m(time) dimensionless;
	when(time=time.min) m=0.01;
	real h(time) dimensionless;
	when(time=time.min) h=0.8;
	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_siCa(time) nanoA;
	real i_siK(time) nanoA;
	real i_siNa(time) nanoA;
	real P_si nanoA_per_millimolar;
	P_si=15;
	real d(time) dimensionless;
	when(time=time.min) d=0.005;
	real f(time) dimensionless;
	when(time=time.min) f=1;
	real f2(time) dimensionless;
	when(time=time.min) f2=1;
	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(time) per_second;
	real beta_f(time) per_second;
	real delta_f millivolt;
	delta_f=0.0001;
	real E0_f(time) millivolt;
	real alpha_f2 per_second;
	alpha_f2=5;
	real beta_f2(time) per_second;
	real K_mf2 millimolar;
	K_mf2=0.001;
	real radius micrometre;
	radius=0.05;
	real length micrometre;
	length=2;
	real V_e_ratio dimensionless;
	V_e_ratio=0.1;
	real V_Cell micrometre3;
	real Vi micrometre3;
	real V_up micrometre3;
	real V_rel micrometre3;
	real i_up(time) nanoA;
	real i_tr(time) nanoA;
	real i_rel(time) nanoA;
	real Ca_up(time) millimolar;
	when(time=time.min) Ca_up=2;
	real Ca_rel(time) millimolar;
	when(time=time.min) Ca_rel=1;
	real Ca_up_max millimolar;
	Ca_up_max=5;
	real K_mCa millimolar;
	K_mCa=0.001;
	real p(time) dimensionless;
	when(time=time.min) p=1;
	real alpha_p(time) per_second;
	real beta_p(time) per_second;
	real tau_up second;
	tau_up=0.025;
	real tau_rep second;
	tau_rep=2;
	real tau_rel second;
	tau_rel=0.05;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real rCa dimensionless;
//	rCa=2;
	real Ve micrometre3;
	Ve=0.00157;
	real Kb millimolar;
	Kb=4;
	real i_mK(time) nanoA;
	real pf per_second;
	pf=0.7;

	// <component name="environment">

	// <component name="membrane">
	RTONF=(R*T/F);
	V:time=((-1)*(i_f+i_K+i_K1+i_to+i_Na_b+i_Ca_b+i_p+i_NaCa+i_Na+i_si+i_pulse)/C);

	// <component name="hyperpolarising_activated_current">
	E_Na=(RTONF*ln(Nao/Nai));
	E_K=(RTONF*ln(Kc/Ki));
	i_fNa=(y*Kc/(Kc+Km_f)*g_f_Na*(V-E_Na));
	i_fK=(y*Kc/(Kc+Km_f)*g_f_K*(V-E_K));
	i_f=(i_fNa+i_fK);

	// <component name="hyperpolarising_activated_current_y_gate">
	alpha_y=((.05 per_second)*exp((-1)*(.067 per_millivolt)*(V+(52 millivolt)-(10 millivolt))));
	E0_y=(V+(52 millivolt)-(10 millivolt));
	beta_y=(if (abs(E0_y)<delta_y) (2.5 per_second) else (1 per_millivolt_second)*E0_y/(1-exp((-1)*(.2 per_millivolt)*E0_y)));
	y:time=(alpha_y*(1-y)-beta_y*y);

	// <component name="time_dependent_potassium_current">
	I_K=(i_K_max*(Ki-Kc*exp((-1)*V/RTONF))/(140 millimolar));
	i_K=(x*I_K);

	// <component name="time_dependent_potassium_current_x_gate">
	alpha_x=((.5 per_second)*exp((.0826 per_millivolt)*(V+(50 millivolt)))/(1+exp((.057 per_millivolt)*(V+(50 millivolt)))));
	beta_x=((1.3 per_second)*exp((-1)*(.06 per_millivolt)*(V+(20 millivolt)))/(1+exp((-1)*(.04 per_millivolt)*(V+(20 millivolt)))));
	x:time=(alpha_x*(1-x)-beta_x*x);

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

	// <component name="transient_outward_current">
	i_to=(s*g_to*(.2+Kc/(Km_to+Kc))*Cai/(Km_Ca+Cai)*(V+(10 millivolt))/(1-exp((-1)*(.2 per_millivolt)*(V+(10 millivolt))))*(Ki*exp(.5*V/RTONF)-Kc*exp((-1)*.5*V/RTONF)));

	// <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)*(V+(10 millivolt))/(8 millivolt))));
	s:time=(alpha_s*(1-s)-beta_s*s);

	// <component name="sodium_background_current">
	i_Na_b=(g_Nab*(V-E_Na));

	// <component name="calcium_background_current">
	E_Ca=(.5*RTONF*ln(Cao/Cai));
	i_Ca_b=(g_Cab*(V-E_Ca));

	// <component name="sodium_potassium_pump">
	i_p=(I_p*Kc/(K_mK+Kc)*Nai/(K_mNa+Nai));

	// <component name="Na_Ca_exchanger">
	i_NaCa=(K_NaCa*(exp(gamma*(3-2)*V/RTONF)*Nai^3*Cao-exp((gamma-1)*(3-2)*V/RTONF)*Nao^3*Cai)/(((1 millimolar4)+d_NaCa*(Cai*Nao^3+Cao*Nai^3))*(1+Cai/(.0069 millimolar))));

	// <component name="fast_sodium_current">
	E_mh=(RTONF*ln((Nao+.12*Kc)/(Nai+.12*Ki)));
	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)/(320*exp((-1)*(.1 per_millivolt)*(V+(75 millivolt)))+1));
	h:time=(alpha_h*(1-h)-beta_h*h);

	// <component name="second_inward_current">
	i_siCa=(4*P_si*(V-(50 millivolt))/(RTONF*(1-exp((-1)*1*(V-(50 millivolt))*2/RTONF)))*(Cai*exp((100 millivolt)/RTONF)-Cao*exp((-1)*2*(V-(50 millivolt))/RTONF))*d*f*f2);
	i_siK=(.01*P_si*(V-(50 millivolt))/(RTONF*(1-exp((-1)*1*(V-(50 millivolt))/RTONF)))*(Ki*exp((50 millivolt)/RTONF)-Kc*exp((-1)*1*(V-(50 millivolt))/RTONF))*d*f*f2);
	i_siNa=(.01*P_si*(V-(50 millivolt))/(RTONF*(1-exp((-1)*1*(V-(50 millivolt))/RTONF)))*(Nai*exp((50 millivolt)/RTONF)-Nao*exp((-1)*1*(V-(50 millivolt))/RTONF))*d*f*f2);
	i_si=(i_siCa+i_siK+i_siNa);

	// <component name="second_inward_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_current_f_gate">
	E0_f=(V+(34 millivolt));
	alpha_f=(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=((50 per_second)/(1+exp((-1)*1*(V+(34 millivolt))/(4 millivolt))));
	f:time=(alpha_f*(1-f)-beta_f*f);

	// <component name="second_inward_current_f2_gate">
	f2:time=(alpha_f2-f2*(alpha_f2+beta_f2));
	beta_f2=(Cai*alpha_f2/K_mf2);

	// <component name="extracellular_sodium_concentration">

	// <component name="intracellular_sodium_concentration">
	V_Cell=(3.141592654*radius^2*length);
	Vi=(V_Cell*(1-V_e_ratio));
	Nai:time=((-1)*(1 mA_nA)*(i_Na+i_Na_b+i_fNa+i_siNa+i_p*3+i_NaCa*3/(3-2))/((1 litre_micrometre3)*Vi*F));

	// <component name="extracellular_calcium_concentration">

	// <component name="intracellular_calcium_concentration">
	i_up=(2*(1 litre_micrometre3)*Vi*F/((1 mA_nA)*tau_up*Ca_up_max)*Cai*(Ca_up_max-Ca_up));
	i_tr=(2*(1 litre_micrometre3)*V_rel*F/((1 mA_nA)*tau_rep)*p*(Ca_up-Ca_rel));
	alpha_p=((.625 per_millivolt_second)*(V+(34 millivolt))/(exp((V+(34 millivolt))/(4 millivolt))-1));
	beta_p=((5 per_second)/(1+exp((-1)*1*(V+(34 millivolt))/(4 millivolt))));
	i_rel=(2*(1 litre_micrometre3)*V_rel*F/((1 mA_nA)*tau_rel)*Ca_rel*Cai^2/(Cai^2+K_mCa^2));
	p:time=(alpha_p*(1-p)-beta_p*p);
	Ca_up:time=((1 mA_nA)*(i_up-i_tr)/(2*(1 litre_micrometre3)*V_up*F));
	V_up=(Vi*.05);
	Ca_rel:time=((1 mA_nA)*(i_tr-i_rel)/(2*(1 litre_micrometre3)*V_rel*F));
	V_rel=(Vi*.02);
	Cai:time=((-1)*(1 mA_nA)*(i_siCa+i_Ca_b-2*i_NaCa/(3-2)-i_rel+i_up)/(2*(1 litre_micrometre3)*Vi*F));

	// <component name="extracellular_potassium_concentration">
	i_mK=(i_K1+i_K+i_fK+i_siK+i_to-2*i_p);
	Kc:time=((-1)*pf*(Kc-Kb)+(1 mA_nA)*i_mK/((1 litre_micrometre3)*Ve*F));

	// <component name="intracellular_potassium_concentration">
	Ki:time=((-1)*(1 mA_nA)*i_mK/((1 litre_micrometre3)*Vi*F));
}