/*
 * Mathematical modeling and qualitative analysis of insulin therapies
 * 
 * Model Status
 * 
 * This CellML version of the model has been checked in COR and
 * OpenCell and is known to be valid CellML. The units are consistent
 * and the model runs to recreate the published results.
 * 
 * Model Structure
 * 
 * ABSTRACT: A model of the rabbit sinoatrial action potential
 * is introduced, based on a model by Morris and Lecar. One cell
 * is described by two nonlinear first-order ordinary differential
 * equations, with ten constant parameters. The model is much simpler
 * than most other models in use, but can reproduce perfectly experimentally
 * recorded action potentials. The dynamics of two coupled cells,
 * with and without the presence of periodic acetylcholine pulses,
 * shows examples of bifurcations and strange attractors, mathematical
 * phenomena characterizing chaotic motion. It remains to be clarified
 * whether such dynamics is actually observed, for example in the
 * small irregular variations of the normal heart rate.
 * 
 * model diagram
 * 
 * [[Image file: endresen_1996.png]]
 * 
 * Schematic diagram of the cell model
 * 
 * The original paper reference is cited below:
 * 
 * Chaos in weakly-coupled pacemaker cells, Lars Petter Endresen,
 * 1997, Journal of Theoretical Biology, 184, 41-50. PubMed ID:
 * 9039399
 */

import nsrunit;
unit conversion on;
unit per_second=1 second^(-1);
// unit millivolt predefined
unit per_millivolt=1E3 kilogram^(-1)*meter^(-2)*second^3*ampere^1;
unit picoS=1E-12 kilogram^(-1)*meter^(-2)*second^3*ampere^2;
unit picoF=1E-12 kilogram^(-1)*meter^(-2)*second^4*ampere^2;
unit femtoA=1E-15 ampere^1;
// unit molar predefined

math main {
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real V(time) millivolt;
	when(time=time.min) V=-52.07606;
	real Cm picoF;
	Cm=60;
	real i_s(time) femtoA;
	real i_K(time) femtoA;
	real i_K_ACh(time) femtoA;
	real i_j(time) femtoA;
	real g_s picoS;
	g_s=382.9118;
	real V_s millivolt;
	V_s=214.1429;
	real V_1 millivolt;
	V_1=-35.9358;
	real V_2 millivolt;
	V_2=7.8589;
	real g_K picoS;
	g_K=536.1093;
	real V_K millivolt;
	V_K=-259.0783;
	real w(time) dimensionless;
	when(time=time.min) w=0.0008971;
	real lambda_w per_second;
	lambda_w=20.7796;
	real V_3 millivolt;
	V_3=-27.9375;
	real V_4 millivolt;
	V_4=6.321;
	real u(time) dimensionless;
	when(time=time.min) u=0.2344555;
	real alpha per_second;
	real beta(time) per_second;
	real ACh molar;
	ACh=1e-6;
	real g_j picoS;
	g_j=0;
	real V_B millivolt;
	V_B=-50;

	// <component name="environment">

	// <component name="membrane">
	V:time=((-1)*(i_s+i_K+i_K_ACh+i_j)/Cm);

	// <component name="calcium_channel">
	i_s=(1/2*g_s*(1+tanh((V-V_1)/V_2))*(V-V_s));

	// <component name="potassium_channel">
	i_K=(g_K*w*(V-V_K));

	// <component name="potassium_channel_w_gate">
	w:time=(lambda_w*cosh((V-V_3)/(2*V_4))*(1/2*(1+tanh((V-V_3)/V_4))-w));

	// <component name="acetyl_choline_activated_potassium_channel">
	i_K_ACh=(1*(.27 picoS)*u*(V+(90 millivolt)));

	// <component name="acetyl_choline_activated_potassium_channel_u_gate">
	alpha=((.012332 per_second)/(1+(4.2E-6 molar)/ACh));
	beta=((.01 per_second)*exp((.0133 per_millivolt)*(V+(40 millivolt))));
	u:time=(alpha*(1-u)-beta*u);

	// <component name="coupling_current">
	i_j=(g_j*(V-V_B));
}