/*
 * The Drouhard-Roberge Sodium Current Model (1987)
 * 
 * Model Status
 * 
 * This model is valid CellML but can not be integrated.
 * 
 * Model Structure
 * 
 * In 1987, Jean-Pierre Drouhard and Fernand A. Roberge published
 * a revised formulation of the Hodgkin-Huxley representation of
 * the sodium current in ventricular myocardial cells. Their improvements
 * have ensured that the action potential upstroke is much faster,
 * and in agreement with experimental observations the peak depolarisation
 * is close to the sodium potential. Like the Ebihara-Johnson model
 * (1980), the Drouhard-Roberge model can be used to replace the
 * sodium kinetics of the Beeler-Reuter model (1977).
 * 
 * The complete original paper reference is cited below:
 * 
 * Revised Formulation of the Hodgkin-Huxley Representation of
 * the Sodium Current in Cardiac Cells, Jean-Pierre Drouhard and
 * Fernand A. Roberge, 1987, Computers and Biomedical Research
 * , 20, 333-350. PubMed ID: 3621918
 */

import nsrunit;
// Warning: unit conversion turned off due to unit errors in 8 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:  V, 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;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) V=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=0.15;
	real E_Na mV;
	E_Na=40.0;
	real m(time) dimensionless;
	when(time=time.min) m=0.0;
	real h(time) dimensionless;
	when(time=time.min) h=0.99;
	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=((.9 per_mV_ms)*(V+(42.65 mV))/(1-exp((-0.22)*(V+(42.65 mV)))));
	beta_m=((1.437 per_ms)*exp((-0.085)*(V+(39.75 mV))));
	m:time=(alpha_m*((1 per_mV)-m)-beta_m*m);

	// <component name="fast_sodium_current_h_gate">
	alpha_h=((.1 per_ms)*exp((-0.193)*(V+(79.65 mV))));
	beta_h=(1.7/(1+exp((-0.095)*(V+(20.5 mV)))));
	h:time=(alpha_h*(1-h)-beta_h*h);

	// <component name="slow_inward_current">
	E_s=(((-82.3 mV))-(13.0287 mV)*ln((.001 per_concentration_units)*Cai));
	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=(.0034999999999999996*((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)))))));
}