/*
 * An Electromechanical Model of Excitable Tissue to Study Reentrant
 * Cardiac Arrhythmias
 * 
 * Model Status
 * 
 * This is the original unchecked version of the model imported
 * from the previous CellML model repository, 24-Jan-2006.
 * 
 * Model Structure
 * 
 * Sudden cardiac death is the leading cause of premature death
 * in the developed world. Underlying these potentially lethal
 * cardiac arrhythmias are reentrant electrical sources, or rotating
 * spiral waves of excitation, and it has been shown that each
 * class of cardiac arrhythmia is associated with a particular
 * type of spiral wave. Sudden cardiac death during ventricular
 * fibrillation results in a compromised mechanical pump function
 * of the heart. Disrupted patterns of electrical excitation during
 * ventricular fibrillation results in desynchronised cardiac contractions.
 * 
 * Mathematical modelling is a useful technique for investigating
 * the electrical activity of the heart. Over the past 40 years
 * models have also been used to study the spatio-temporal dynamics
 * of reentrant arrhythmias in 2D, 3D, and in anatomically accurate
 * geometric models of the heart. The mathematical models are based
 * on physiological data, and in parallel with the increasing sophistication
 * of experimental techniques, the mathematical models has also
 * evolved in their complexity and biological accuracy. Previously,
 * mathematical models have tended to focus on either the electrical
 * properties of the heart, or on cardiac mechanics. In the study
 * described here, Nash and Panfilov establish a mathematical model
 * which combines the effects of cardiac mechanics and electrical
 * activity during arrhythmia.
 * 
 * The model described here in CellML is based on the Nash and
 * Panfilov 2004 paper, which describes an electromechanical model
 * of excitable tissue, which is used to study reentrant cardiac
 * arrhythmias. This mathematical model is a version of the modified
 * FitzHugh-Nagumo model, published by Aliev and Panfilov in 1996
 * (for details on the original model please see The FitzHugh-Nagumo
 * Model, 1961). While the original two-variable model described
 * a non-dimensional activation variable (u) and a non-dimensional
 * recovery variable (v), here Nash and Panfilov formulate the
 * model in terms of the real action potential given by the time
 * course of the transmembrane potential (Vm). In so doing, the
 * time rate of change of the activation variable describes the
 * total ionic current through the membrane with the original model
 * parameters adjusted to give the correct dimensionality. This
 * 1996 version of the FitzHugh-Nagumo model has been further modified
 * by Nash and Panfilov in this 2004 study to include the simplest
 * model of active tension development.
 * 
 * The main equations in the paper are Eqn 22abc and Eqn 23. When
 * the first author of this paper converted the model into CellML,
 * he noticed two main errors in the original paper:
 * 
 * Eqn 22b: b should be a, and
 * 
 * Eqn 23: the factor 10 should be on first case line (V is less
 * than 0.05). The second case should have the factor 1.
 * 
 * A new set of model parameters is also used.
 * 
 * The complete original paper reference is cited below:
 * 
 * Electromechanical model of excitable tissue to study reentrant
 * cardiac arrhythmias, M. P. Nash and A. V. Panfilov, 2004, Progress
 * in Biophysics and Molecular Biology, 85, 501-522. PubMed ID:
 * 15142759
 * 
 * reaction diagram
 * 
 * [[Image file: nash_2004.png]]
 * 
 * Properties of the action potential described by the model.
 */

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

math main {
	//Warning:  the following variables were set 'extern' or given
	//  an initial value of '0' because the model would otherwise be
	//  underdetermined:  Istim
	realDomain t ms;
	t.min=0;
	extern t.max;
	extern t.delta;
	real Cm uFpmmsq;
	Cm=1.0;
	real Vr mV;
	Vr=-80.0;
	real Vth mV;
	Vth=-70.0;
	real Vp mV;
	Vp=20.0;
	real k uApmmsq;
	k=8.0;
	real epsilon pms;
	epsilon=0.01;
	real mu1 pms;
	mu1=0.2;
	real mu2 dimensionless;
	mu2=0.3;
	real e0 dimensionless;
	e0=1.0;
	real kTa kPa;
	kTa=47.9;
	extern real Istim uApmmsq;
	real Vm(t) mV;
	when(t=t.min) Vm=-85;
	real Iion(t) uApmmsq;
	real r(t) dimensionless;
	when(t=t.min) r=0.0;
	real Ta(t) kPa;
	when(t=t.min) Ta=0.0;
	real IStimC uApmmsq;
	real u(t) dimensionless;
	real a dimensionless;
	real eps(t) pms;
	real e(t) dimensionless;

	// <component name="interface">
	IStimC=Istim;

	// <component name="membrane_potential">
	Vm:t=((Istim-Iion)/Cm);
	u=((Vm-Vr)/(Vp-Vr));
	a=((Vth-Vr)/(Vp-Vr));

	// <component name="ionic_current">
	Iion=(k*u*(u-a)*(u-1)+r*u);

	// <component name="recovery_variable">
	r:t=(eps*((-1)*r-k*u*(u-(a+1))));
	eps=(epsilon+mu1*r/(mu2+u));

	// <component name="active_tension">
	Ta:t=(e*(kTa*u-Ta));
	e=(if (u<.05) 10*e0 else 1*e0);
}