/*
 * Voltage Clamp Experiments in Striated Muscle Fibres
 * 
 * Model Status
 * 
 * This model has had a stimulus protocol added to it to allow
 * simulation of action potentials. Unfortunately, the details
 * of the stimulation defined in the original paper are not known
 * and as such, the parameters of the stimulus (magnitude, duration,
 * dimensions,) may not be appropriate. Currently, however, the
 * model runs in PCEnv and is able to simulate a train of action
 * potentials. COR will not currently run this model.
 * 
 * ValidateCellML detects unit inconsistencies within this model.
 * 
 * Model Structure
 * 
 * In this 1970 publication, Adrian, Chandler and Hodgkin developed
 * a mathematical model of an action potential in frog striated
 * muscle. Their model was based on data from voltage-clamp experiments,
 * and they include descriptions of three currents:
 * 
 * INa, a fast inward sodium current;
 * 
 * IK, a slow outward potassium current; and
 * 
 * IL, a leak current.
 * 
 * The style of the model equations is based on the The Hodgkin-Huxley
 * Squid Axon Model, 1952. However, the authors acknowledge that
 * the situation in striated muscle is complicated by tubular resistance
 * and capacity. The transverse tubular system which exists in
 * striated myocytes is represented in the model by a linear resistance
 * and capacity in series.
 * 
 * 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 model has been described here in CellML (the raw CellML
 * description of the Adrian et al. 1970 model can be downloaded
 * in various formats as described in ).
 * 
 * As the paper was published in 1970, there is no online version.
 * However, the complete reference is cited below:
 * 
 * Voltage Clamp Experiments in Striated Muscle Fibres. R. H. Adrian,
 * W. K. Chandler, and A. L. Hodgkin. Journal of Physiology, (1970),
 * 208, pp 607-644. PubMed ID: 5499787
 * 
 * In order to complete the model and run simulations, some parameters
 * were also taken from the following paper:
 * 
 * Reconstruction of the Action Potential of Frog Sartorius Muscle.
 * R. H. Adrian and L. D. Peachey. Journal of Physiology, (1973),
 * 235, pp 103-131. PubMed ID: 4778131
 */

import nsrunit;
// Warning: unit conversion turned off due to unit errors in 4 equation(s)
unit conversion off;
unit mV=.001 kilogram^1*meter^2*second^(-3)*ampere^(-1);
unit uApmmsq=1 meter^(-2)*ampere^1;
unit uApmmcu=1E3 meter^(-3)*ampere^1;
unit uFpmmsq=.001 kilogram^(-1)*meter^(-3)*second^4*ampere^2;
unit mSpmmsq=1E3 kilogram^(-1)*meter^(-4)*second^3*ampere^2;
unit mmsqpmS=.001 kilogram^1*meter^4*second^(-3)*ampere^(-2);
unit ms=.001 second^1;
unit pms=1E3 second^(-1);
unit pmm=1E3 meter^(-1);

math main {
	realDomain t ms;
	t.min=0;
	extern t.max;
	extern t.delta;
	real Cm uFpmmsq;
	Cm=0.009;
	real Am pmm;
	Am=200.0;
	real Istim(t) uApmmcu;
	real Vm(t) mV;
	when(t=t.min) Vm=-95.0;
	real Vt(t) mV;
	when(t=t.min) Vt=-95.0;
	real m(t) dimensionless;
	when(t=t.min) m=0.0;
	real h(t) dimensionless;
	when(t=t.min) h=1.0;
	real n(t) dimensionless;
	when(t=t.min) n=0.0;
	real INa(t) uApmmsq;
	real IK(t) uApmmsq;
	real IL(t) uApmmsq;
	real IT(t) uApmmsq;
	real IStimC(t) uApmmcu;
	real AmC pmm;
	real IstimStart ms;
	IstimStart=10;
	real IstimEnd ms;
	IstimEnd=50000;
	real IstimAmplitude uApmmcu;
	IstimAmplitude=0.5;
	real IstimPeriod ms;
	IstimPeriod=1000;
	real IstimPulseDuration ms;
	IstimPulseDuration=1;
	real gNa_max mSpmmsq;
	gNa_max=1.8;
	real ENa mV;
	ENa=50.0;
	real alpha_m(t) pms;
	real beta_m(t) pms;
	real alpha_m_max pms;
	alpha_m_max=0.208;
	real beta_m_max pms;
	beta_m_max=2.081;
	real Em mV;
	Em=-42.0;
	real v_alpha_m dimensionless;
	v_alpha_m=10.0;
	real v_beta_m mV;
	v_beta_m=18.0;
	real alpha_h(t) pms;
	real beta_h(t) pms;
	real alpha_h_max pms;
	alpha_h_max=0.0156;
	real beta_h_max pms;
	beta_h_max=3.382;
	real Eh mV;
	Eh=-41.0;
	real v_alpha_h mV;
	v_alpha_h=14.7;
	real v_beta_h mV;
	v_beta_h=7.6;
	real gK_max mSpmmsq;
	gK_max=0.415;
	real EK mV;
	EK=-70.0;
	real alpha_n(t) pms;
	real beta_n(t) pms;
	real alpha_n_max pms;
	alpha_n_max=0.0229;
	real beta_n_max pms;
	beta_n_max=0.09616;
	real En mV;
	En=-40.0;
	real v_alpha_n dimensionless;
	v_alpha_n=7.0;
	real v_beta_n mV;
	v_beta_n=40.0;
	real EL mV;
	EL=-95.0;
	real gL_max mSpmmsq;
	gL_max=0.0024;
	real Rs mmsqpmS;
	Rs=15.0;
	real Ct uFpmmsq;
	Ct=0.04;

	// <component name="interface">
	IStimC=Istim;
	AmC=Am;
	Istim=(if (((t>=IstimStart) and (t<=IstimEnd)) and ((t-IstimStart-floor((t-IstimStart)/IstimPeriod)*IstimPeriod)<=IstimPulseDuration)) IstimAmplitude else (0 uApmmcu));

	// <component name="membrane">
	Vm:t=((Istim-(INa+IK+IL+IT))/Cm);

	// <component name="sodium_current">
	INa=(gNa_max*m*m*m*h*(Vm-ENa));

	// <component name="sodium_current_m_gate">
	alpha_m=(alpha_m_max*(Vm-Em)/((1 mV)-exp((Em-Vm)/v_alpha_m)));
	beta_m=(beta_m_max*exp((Em-Vm)/v_beta_m));
	m:t=(alpha_m*(1-m)-beta_m*m);

	// <component name="sodium_current_h_gate">
	alpha_h=(alpha_h_max*exp((Eh-Vm)/v_alpha_h));
	beta_h=(beta_h_max/(1+exp((Eh-Vm)/v_beta_h)));
	h:t=(alpha_h*(1-h)-beta_h*h);

	// <component name="potassium_current">
	IK=(gK_max*n*n*n*n*(Vm-EK));

	// <component name="potassium_current_n_gate">
	alpha_n=(alpha_n_max*(Vm-En)/(1-exp((En-Vm)/v_alpha_n)));
	beta_n=(beta_n_max*exp((En-Vm)/v_beta_n));
	n:t=(alpha_n*(1-n)-beta_n*n);

	// <component name="leak_current">
	IL=(gL_max*(Vm-EL));

	// <component name="Ttubular_current">
	IT=((Vm-Vt)/Rs);

	// <component name="Ttubular_current_Vt_var">
	Vt:t=((Vm-Vt)/(Rs*Ct));
}