/*
 * A Simplified Ventricular Myocyte Model
 * 
 * Model Status
 * 
 * This version of this model is known to run in both OpenCell
 * and COR. The units have been checked and they are consistent.
 * A generic stimulus protocol has been added to allow the model
 * to simulate trains of action potentials. Although the model
 * does run, the simulation output is still not quite the same
 * as the original published model. The original model authors
 * have been contacted and we will continue to curate the CellML
 * model.
 * 
 * Model Structure
 * 
 * ABSTRACT: In this paper we introduce and study a model for electrical
 * activity of cardiac membrane which incorporates only an inward
 * and an outward current. This model is useful for three reasons:
 * (1) Its simplicity, comparable to the FitzHugh-Nagumo model,
 * makes it useful in numerical simulations, especially in two
 * or three spatial dimensions where numerical efficiency is so
 * important. (2) It can be understood analytically without recourse
 * to numerical simulations. This allows us to determine rather
 * completely how the parameters in the model affect its behavior
 * which in turn provides insight into the effects of the many
 * parameters in more realistic models. (3) It naturally gives
 * rise to a one-dimensional map which specifies the action potential
 * duration as a function of the previous diastolic interval. For
 * certain parameter values, this map exhibits a new phenomenon--subcritical
 * alternans--that does not occur for the commonly used exponential
 * map.
 * 
 * The original paper reference is cited below:
 * 
 * A two-current model for the dynamics of cardiac membrane, Colleen
 * C. Mitchell, David G. Schaeffer, 2003, Bulletin of Mathematical
 * Biology, 65, (5), 767-793. PubMed ID: 12909250
 * 
 * cell diagram
 * 
 * [[Image file: mitchell_2003.png]]
 * 
 * A schematic diagram of the two ionic currents described by the
 * Mitchell-Schaeffer model of a ventricular myocyte.
 */

import nsrunit;
unit conversion on;
unit ms=.001 second^1;
unit per_ms=1E3 second^(-1);

math main {
	realDomain time ms;
	time.min=0;
	extern time.max;
	extern time.delta;
	real J_stim(time) per_ms;
	real IstimStart ms;
	IstimStart=0;
	real IstimEnd ms;
	IstimEnd=50000;
	real IstimAmplitude per_ms;
	IstimAmplitude=0.2;
	real IstimPeriod ms;
	IstimPeriod=500;
	real IstimPulseDuration ms;
	IstimPulseDuration=1;
	real Vm(time) dimensionless;
	when(time=time.min) Vm=0.00000820413566106744;
	real J_in(time) per_ms;
	real J_out(time) per_ms;
	real tau_in ms;
	tau_in=0.3;
	real h(time) dimensionless;
	when(time=time.min) h=0.8789655121804799;
	real tau_open ms;
	tau_open=120.0;
	real tau_close ms;
	tau_close=150.0;
	real V_gate dimensionless;
	V_gate=0.13;
	real tau_out ms;
	tau_out=6.0;

	// <component name="environment">

	// <component name="J_stim">
	J_stim=(if (((time>=IstimStart) and (time<=IstimEnd)) and ((time-IstimStart-floor((time-IstimStart)/IstimPeriod)*IstimPeriod)<=IstimPulseDuration)) IstimAmplitude else (0 per_ms));

	// <component name="membrane">
	Vm:time=(J_in+J_out+J_stim);

	// <component name="J_in">
	J_in=(h*(Vm^2*(1-Vm))/tau_in);

	// <component name="J_in_h_gate">
	h:time=(if (Vm<V_gate) (1-h)/tau_open else (-1)*h/tau_close);

	// <component name="J_out">
	J_out=((-1)*(Vm/tau_out));
}