/*
 * Model Status
 * 
 * This version has had a stimulus protocol component added to
 * allow the model to simulate multiple action potentials, has
 * been unit checked and curated and is known to run in COR and
 * PCEnv.
 * 
 * Model Structure
 * 
 * In a series of papers published in 1952, A.L. Hodgkin and A.F.
 * Huxley presented the results of a series of experiments in which
 * they investigated the flow of electric current through the surface
 * membrane of the giant nerve fibre of a squid. In the summary
 * paper of the Hodgkin and Huxley model, the authors developed
 * a mathematical description of the behaviour of the membrane
 * based upon these experiments, which accounts for the conduction
 * and excitation of the fibre. The form of this description has
 * been used as the basis for almost all other ionic current models
 * of excitable tissues, including Purkinje fibres and cardiac
 * atrial and ventricular muscle.
 * 
 * The summary paper is cited below:
 * 
 * A Quantitative Description of Membrane Current and its Application
 * to Conduction and Excitation in Nerve, A. L. Hodgkin and A.
 * F. Huxley, 1952, The Journal of Physiology, 117, 500-544. PubMed
 * ID: 12991237
 * 
 * Electrical circuit describing the current across the cell membrane
 * 
 * [[Image file: hodgkin_1952.png]]
 * 
 * A schematic cell diagram describing the current flows across
 * the cell membrane that are captured in the Hodgkin Huxley model.
 */

import nsrunit;
unit conversion on;
// unit millisecond predefined
unit per_millisecond=1E3 second^(-1);
// unit millivolt predefined
unit per_millivolt_millisecond=1E6 kilogram^(-1)*meter^(-2)*second^2*ampere^1;
unit milliS_per_cm2=10 kilogram^(-1)*meter^(-4)*second^3*ampere^2;
unit microF_per_cm2=.01 kilogram^(-1)*meter^(-4)*second^4*ampere^2;
unit microA_per_cm2=.01 meter^(-2)*ampere^1;

math main {
	realDomain time millisecond;
	time.min=0;
	extern time.max;
	extern time.delta;
	real V(time) millivolt;
	when(time=time.min) V=-75;
	real E_R millivolt;
	E_R=-75;
	real Cm microF_per_cm2;
	Cm=1;
	real i_Na(time) microA_per_cm2;
	real i_K(time) microA_per_cm2;
	real i_L(time) microA_per_cm2;
	real i_Stim(time) microA_per_cm2;
	real g_Na milliS_per_cm2;
	g_Na=120;
	real E_Na millivolt;
	real m(time) dimensionless;
	when(time=time.min) m=0.05;
	real h(time) dimensionless;
	when(time=time.min) h=0.6;
	real alpha_m(time) per_millisecond;
	real beta_m(time) per_millisecond;
	real alpha_h(time) per_millisecond;
	real beta_h(time) per_millisecond;
	real g_K milliS_per_cm2;
	g_K=36;
	real E_K millivolt;
	real n(time) dimensionless;
	when(time=time.min) n=0.325;
	real alpha_n(time) per_millisecond;
	real beta_n(time) per_millisecond;
	real g_L milliS_per_cm2;
	g_L=0.3;
	real E_L millivolt;

	// <component name="environment">

	// <component name="membrane">
	i_Stim=(if ((time>=(10 millisecond)) and (time<=(10.5 millisecond))) (20 microA_per_cm2) else (0 microA_per_cm2));
	V:time=((-1)*((-1)*i_Stim+i_Na+i_K+i_L)/Cm);

	// <component name="sodium_channel">
	E_Na=(E_R+(115 millivolt));
	i_Na=(g_Na*m^3*h*(V-E_Na));

	// <component name="sodium_channel_m_gate">
	alpha_m=((-1)*(.1 per_millivolt_millisecond)*(V+(50 millivolt))/(exp((-1)*(V+(50 millivolt))/(10 millivolt))-1));
	beta_m=((4 per_millisecond)*exp((-1)*(V+(75 millivolt))/(18 millivolt)));
	m:time=(alpha_m*(1-m)-beta_m*m);

	// <component name="sodium_channel_h_gate">
	alpha_h=((.07 per_millisecond)*exp((-1)*(V+(75 millivolt))/(20 millivolt)));
	beta_h=((1 per_millisecond)/(exp((-1)*(V+(45 millivolt))/(10 millivolt))+1));
	h:time=(alpha_h*(1-h)-beta_h*h);

	// <component name="potassium_channel">
	E_K=(E_R-(12 millivolt));
	i_K=(g_K*n^4*(V-E_K));

	// <component name="potassium_channel_n_gate">
	alpha_n=((-1)*(.01 per_millivolt_millisecond)*(V+(65 millivolt))/(exp((-1)*(V+(65 millivolt))/(10 millivolt))-1));
	beta_n=((.125 per_millisecond)*exp((V+(75 millivolt))/(80 millivolt)));
	n:time=(alpha_n*(1-n)-beta_n*n);

	// <component name="leakage_current">
	E_L=(E_R+(10.613 millivolt));
	i_L=(g_L*(V-E_L));
}