/*
 * The Ebihara-Johnson Sodium Current Model (1980)
 * 
 * Model Status
 * 
 * This model is valid CellML and contains consistent units but
 * is not fully parameterised.
 * 
 * Model Structure
 * 
 * In 1980, Lisa Ebihara and Edward A. Johnson published the first
 * electrophysiological model to specifically target a single ion
 * channel and attempt to quantify its parameters. They used the
 * Hodgkin-Huxley formulation to revise the parameters of the fast
 * sodium current in cardiac muscle. This model exhibits a faster
 * sodium current than the Beeler-Reuter model, and when the two
 * models are coupled, the Ebihara-Johnson can be used as a direct
 * replacement of the sodium kinetics of the Beeler-Reuter.
 * 
 * model diagram Schematic diagram of the Ebihara and Johnson model.
 * 
 * The complete original paper reference is cited below:
 * 
 * Fast Sodium Current In Cardiac Muscle, Lisa Ebihara and Edward
 * A. Johnson, 1980, Biophys. J. 32, 779-790. PubMed ID: 7260301
 */

import nsrunit;
// Warning: unit conversion turned off due to unit errors in 6 equation(s)
unit conversion off;
unit ms=.001 second^1;
unit per_ms=1E3 second^(-1);
unit mV=.001 kilogram^1*meter^2*second^(-3)*ampere^(-1);
unit per_mV=1E3 kilogram^(-1)*meter^(-2)*second^3*ampere^1;
unit per_mV_ms=1E6 kilogram^(-1)*meter^(-2)*second^2*ampere^1;
unit mS_per_mm2=1E3 kilogram^(-1)*meter^(-4)*second^3*ampere^2;
unit uF_per_mm2=1 kilogram^(-1)*meter^(-4)*second^4*ampere^2;
unit uA_per_mm2=1 meter^(-2)*ampere^1;
unit concentration_units=1 meter^(-3)*mole^1;
unit per_concentration_units=1 meter^3*mole^(-1);

math main {
	//Warning:  the following variables were set 'extern' or given
	//  an initial value of '0' because the model would otherwise be
	//  underdetermined:  Istim, time_dependent_outward_current.alpha_x1,
	//  time_dependent_outward_current.beta_x1, x1
	realDomain time ms;
	time.min=0;
	extern time.max;
	extern time.delta;
	real V(time) mV;
	when(time=time.min) V=-87.0;
	real C uF_per_mm2;
	C=0.013;
	real i_Na(time) uA_per_mm2;
	real i_s(time) uA_per_mm2;
	real i_x1(time) uA_per_mm2;
	real i_K1(time) uA_per_mm2;
	extern real Istim uA_per_mm2;
	real IStimC uA_per_mm2;
	real g_Na mS_per_mm2;
	g_Na=23.0e-2;
	real E_Na mV;
	E_Na=29.0;
	real m(time) dimensionless;
	when(time=time.min) m=0.0;
	real h(time) dimensionless;
	when(time=time.min) h=0.18;
	real alpha_m(time) per_ms;
	real beta_m(time) per_ms;
	real alpha_h(time) per_ms;
	real beta_h(time) per_ms;
	real g_s mS_per_mm2;
	g_s=9.0e-4;
	real E_s(time) mV;
	real Cai(time) concentration_units;
	when(time=time.min) Cai=0.000000177;
	real d(time) dimensionless;
	when(time=time.min) d=0.003;
	real f(time) dimensionless;
	when(time=time.min) f=0.994;
	real alpha_d(time) per_ms;
	real beta_d(time) per_ms;
	real alpha_f(time) per_ms;
	real beta_f(time) per_ms;
	extern real time_dependent_outward_current.alpha_x1 per_ms;
	extern real time_dependent_outward_current.beta_x1 per_ms;
	real x1(time) dimensionless;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) x1=0;
	real time_dependent_outward_current_x1_gate.alpha_x1(time) per_ms;
	real time_dependent_outward_current_x1_gate.beta_x1(time) per_ms;

	// <component name="environment">

	// <component name="membrane">
	V:time=((Istim-(i_Na+i_s+i_x1+i_K1))/C);
	IStimC=Istim;

	// <component name="fast_sodium_current">
	i_Na=(g_Na*m^3*h*(V-E_Na));

	// <component name="fast_sodium_current_m_gate">
	alpha_m=((.32 per_mV_ms)*((47.13 mV)+V)/(1-exp((-1)*V-(47.13 mV))));
	beta_m=((.08 per_ms)*exp((-1)*V/(11 mV)));
	m:time=(alpha_m*((1 per_mV)-m)-beta_m*m);

	// <component name="fast_sodium_current_h_gate">
	alpha_h=(if (V<((-40 mV))) (.135 per_ms)*exp((((-80 mV))-V)/(6.8 mV)) else (0 per_ms));
	beta_h=(if (V<((-40 mV))) (3.56 per_ms)*exp((.079 mV)*V)+(3.1E5 per_ms)*exp((.35 per_mV)*V) else 1/((.13 ms)*(1+exp((-1)*(V+(10.66 mV))/(11.1 mV)))));
	h:time=(alpha_h*(1-h)-beta_h*h);

	// <component name="slow_inward_current">
	E_s=(((-82.3 mV))-(13.0287 mV)*ln(Cai*(.001 per_concentration_units)));
	i_s=(g_s*d*f*(V-E_s));
	Cai:time=((-0.01)*i_s+.07*(1E-4-Cai));

	// <component name="slow_inward_current_d_gate">
	alpha_d=((.095 per_ms)*exp((-1)*((V-(5 mV))/(100 mV)))/(1+exp((-1)*((V-(5 mV))/(13.89 mV)))));
	beta_d=((.07 per_ms)*exp((-1)*((V+(44 mV))/(59 mV)))/(1+exp((V+(44 mV))/(20 mV))));
	d:time=(alpha_d*(1-d)-beta_d*d);

	// <component name="slow_inward_current_f_gate">
	alpha_f=((.012 per_ms)*exp((-1)*((V+(28 mV))/(125 mV)))/(1+exp((V+(28 mV))/(6.67 mV))));
	beta_f=((.0065 per_ms)*exp((-1)*((V+(30 mV))/(50 mV)))/(1+exp((-1)*((V+(30 mV))/(5 mV)))));
	f:time=(alpha_f*(1-f)-beta_f*f);

	// <component name="time_dependent_outward_current">
	i_x1=(x1*.008*((exp((.04 per_mV)*(V+(77 mV)))-1)/exp((.04 per_mV)*(V+(35 mV)))));

	// <component name="time_dependent_outward_current_x1_gate">
	time_dependent_outward_current_x1_gate.alpha_x1=((5E-4 per_ms)*(exp((V+(50 mV))/(12.1 mV))/(1+exp((V+(50 mV))/(17.5 mV)))));
	time_dependent_outward_current_x1_gate.beta_x1=((.0013 per_ms)*(exp((-1)*((V+(20 mV))/(16.67 mV)))/(1+exp((-1)*((V+(20 mV))/(25 mV))))));
	x1:time=(time_dependent_outward_current_x1_gate.alpha_x1*(1-x1)-time_dependent_outward_current_x1_gate.beta_x1*x1);

	// <component name="time_independent_outward_current">
	i_K1=(.0035*((4 per_ms)*((exp((.04 per_mV)*(V+(85 mV)))-1)/(exp((.08 per_mV)*(V+(53 mV)))+exp((.04 per_mV)*(V+(53 mV)))))+(.2 per_ms)*((V+(23 mV))/(1-exp(((-0.04 per_mV))*(V+(23 mV)))))));
}