/*
 * Beeler-Reuter Mammalian Ventricular Model 1977
 * 
 * Model Status
 * 
 * This model has been curated by Penny Noble using Flavio Fenton's
 * Java code as a reference (See http://thevirtualheart.org/ for
 * Java applet rendering of model - Java code is available from
 * Dr Fenton.) An artificial stimulus component has been added
 * this model to allow it to reproduce the action potential simulation
 * shown in Figure 4 of the publication. The model is known to
 * run and integrate in the PCEnv and COR CellML environments.
 * A PCEnv session file is also associated with this model.
 * 
 * ValidateCellML detects unit inconsistency within this model.
 * 
 * 
 * 
 * Model Structure
 * 
 * In contrast to the earlier Purkinje fibre ionic current models
 * of D. Noble (1962) and R.E. McAllister, D. Noble and R.W. Tsien
 * (1975), the G.W. Beeler and H. Reuter 1977 model was developed
 * to describe the mammalian ventricular action potential. Not
 * all the ionic currents of the Purkinje fibre model are present
 * in ventricular tissue; therefore, this model is simpler than
 * the MNT model. The total ionic flux is divided into only four
 * discrete, individual ionic currents (see below). The main additional
 * feature of the Beeler-Reuter ionic current model is a representation
 * of the intracellular calcium ion concentration.
 * 
 * The complete original paper reference is cited below:
 * 
 * Reconstruction of the action potential of ventricular myocardial
 * fibres, Beeler, G.W. and Reuter, H. 1977 Journal of Physiology
 * , 268, 177-210. PubMed ID: 874889
 * 
 * cell diagram of the Beeler-Reuter model showing ionic currents
 * across the cell surface membrane
 * 
 * [[Image file: beeler_reuter_1977.png]]
 * 
 * A schematic diagram describing the current flows across the
 * cell membrane that are captured in the BR model.
 * 
 * the cellml rendering of the Beeler-Reuter model
 * 
 * [[Image file: cellml_rendering.gif]]
 * 
 * The network defined in the CellML description of the Beeler-Reuter
 * model. A key describing the significance of the shapes of the
 * components and the colours of the connections between them is
 * in the notation guide. For simplicity, not all the variables
 * are shown.
 * 
 * The membrane physically contains the currents as indicated by
 * the blue arrows in . The currents act independently and are
 * not connected to each other. Several of the channels encapsulate
 * and contain further components which represent activation and
 * inactivation gates. The addition of an encapsulation relationship
 * informs modellers and processing software that the gates are
 * important parts of the current model. It also prevents any other
 * components that aren't also encapsulated by the parent component
 * from connecting to its gates, effectively hiding them from the
 * rest of the model.
 * 
 * The breakdown of the model into components and the definition
 * of encapsulation and containment relationships between them
 * is somewhat arbitrary. When considering how a model should be
 * broken into components, modellers are encouraged to consider
 * which parts of a model might be re-used and how the physiological
 * elements of the system being modelled are naturally bounded.
 * Containment relationships should be used to provide simple rendering
 * information for processing software (ideally, this will correspond
 * to the layout of the physical system), and encapsulation should
 * be used to group sets of components into sub-models.
 */

import nsrunit;
unit conversion on;
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);
unit coulomb_per_mole=1 second^1*ampere^1*mole^(-1);
unit per_mm=1E3 meter^(-1);

math main {
	realDomain time ms;
	time.min=0;
	extern time.max;
	extern time.delta;
	real V(time) mV;
	when(time=time.min) V=-84.624;
	real C uF_per_mm2;
	C=0.01;
	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;
	real Istim(time) uA_per_mm2;
	real g_Na mS_per_mm2;
	g_Na=4e-2;
	real E_Na mV;
	E_Na=50;
	real g_Nac mS_per_mm2;
	g_Nac=3e-5;
	real m(time) dimensionless;
	when(time=time.min) m=0.011;
	real h(time) dimensionless;
	when(time=time.min) h=0.988;
	real j(time) dimensionless;
	when(time=time.min) j=0.975;
	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 alpha_j(time) per_ms;
	real beta_j(time) per_ms;
	real g_s mS_per_mm2;
	g_s=9e-4;
	real E_s(time) mV;
	real Cai(time) concentration_units;
	when(time=time.min) Cai=1e-4;
	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;
	real x1(time) dimensionless;
	when(time=time.min) x1=0.0001;
	real alpha_x1(time) per_ms;
	real beta_x1(time) per_ms;
	real IstimStart ms;
	IstimStart=10;
	real IstimEnd ms;
	IstimEnd=50000;
	real IstimAmplitude uA_per_mm2;
	IstimAmplitude=0.5;
	real IstimPeriod ms;
	IstimPeriod=1000;
	real IstimPulseDuration ms;
	IstimPulseDuration=1;

	// <component name="environment">

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

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

	// <component name="sodium_current_m_gate">
	alpha_m=((-1)*(1 per_mV_ms)*(V+(47 mV))/(exp((-1)*(.1 per_mV)*(V+(47 mV)))-1));
	beta_m=((40 per_ms)*exp((-1)*(.056 per_mV)*(V+(72 mV))));
	m:time=(alpha_m*(1-m)-beta_m*m);

	// <component name="sodium_current_h_gate">
	alpha_h=((.126 per_ms)*exp((-1)*(.25 per_mV)*(V+(77 mV))));
	beta_h=((1.7 per_ms)/(exp((-1)*(.082 per_mV)*(V+(22.5 mV)))+1));
	h:time=(alpha_h*(1-h)-beta_h*h);

	// <component name="sodium_current_j_gate">
	alpha_j=((.055 per_ms)*exp((-1)*(.25 per_mV)*(V+(78 mV)))/(exp((-1)*(.2 per_mV)*(V+(78 mV)))+1));
	beta_j=((.3 per_ms)/(exp((-1)*(.1 per_mV)*(V+(32 mV)))+1));
	j:time=(alpha_j*(1-j)-beta_j*j);

	// <component name="slow_inward_current">
	E_s=((-1)*(82.3 mV)-(13.0287 mV)*ln(Cai*(.001 per_concentration_units)));
	i_s=(g_s*d*f*(V-E_s));
	Cai:time=((-1)*(.01 per_mm)*i_s/(1 coulomb_per_mole)+(.07 per_ms)*((1E-4 concentration_units)-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 uA_per_mm2)*(exp((.04 per_mV)*(V+(77 mV)))-1)/exp((.04 per_mV)*(V+(35 mV))));

	// <component name="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))));
	beta_x1=((.0013 per_ms)*exp((-1)*(V+(20 mV))/(16.67 mV))/(1+exp((-1)*(V+(20 mV))/(25 mV))));
	x1:time=(alpha_x1*(1-x1)-beta_x1*x1);

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

	// <component name="stimulus_protocol">
	Istim=(if (((time>=IstimStart) and (time<=IstimEnd)) and ((time-IstimStart-floor((time-IstimStart)/IstimPeriod)*IstimPeriod)<=IstimPulseDuration)) IstimAmplitude else (0 uA_per_mm2));
}