/*
 * Hodgkin Huxley Squid Axon Model 1952
 * 
 * Model Status
 * 
 * This particular variant of the CellML model is based on the
 * original model in the 1952 Hodgkin-Huxley published paper. Previous
 * versions of the CellML model description have been modified
 * from the original model to to be consistent with the modern
 * convention of describing cardiac models. This particular model
 * has been tested in both PCEnv and COR. To run the model correctly
 * in COR you need to set the duration of the simulation to 50
 * ms and, to make the rendering of the results more accurate,
 * the output to 0.1 ms.
 * 
 * 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 CellML model itself is intended to represent the original
 * model from the published paper. To date, all the other versions
 * of the Hodgkin-Hulxley model have been slightly modified versions
 * of the original published model. In particular the current descriptions
 * were reversed to be consistent with the modern convention proposed
 * by Prof. Denis Noble, now commonly adopted for cardiac muscle
 * model descriptions.
 * 
 * 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.
 * 
 * The complete original paper reference 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
 */

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=0;
	real E_R millivolt;
	E_R=0;
	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))) (-1)*(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 per_millivolt_millisecond)*(V+(25 millivolt))/(exp((V+(25 millivolt))/(10 millivolt))-1));
	beta_m=((4 per_millisecond)*exp(V/(18 millivolt)));
	m:time=(alpha_m*(1-m)-beta_m*m);

	// <component name="sodium_channel_h_gate">
	alpha_h=((.07 per_millisecond)*exp(V/(20 millivolt)));
	beta_h=((1 per_millisecond)/(exp((V+(30 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=((.01 per_millivolt_millisecond)*(V+(10 millivolt))/(exp((V+(10 millivolt))/(10 millivolt))-1));
	beta_n=((.125 per_millisecond)*exp(V/(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));
}