ATS-AP.cpp

/* The modified Luo-Rudy Dynamic (LRd) Model of the mammalian ventricular myocyte */
/* The Markovian model of IKr and IKs channels was incorporated */
/* This code requires a C++ compiler */
 
/* Am J Physiol Heart Circ Physiol 2006 [Epub ahead of print] */ 
/* Implemented by Sheng-Nan Wu and Han-Dong Chang  2006/08/01 */
  
/* IMPORTANT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
USER SHOULD PACE THE MODEL UNTIL STEADY-STATE IS REACHED.  THIS REQUIRES APPROXIMATELY 5-20 MINUTES OF PACING DEPENDING ON THE RATE.*/
 
#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#include <fstream.h>
#include <stdlib.h>
#include <stdio.h>
 
#define bcl 1000   // Basic Cycle Length (ms)
#define beats 51  // Number of Beats

double cc3, cc2, cc1, oo, ii;

double cccc1, cccc2, cccc3, cccc4, cccc5, cccc6, cccc7, cccc8, cccc9, cccc10, cccc11, cccc12, cccc13, cccc14, cccc15, oooo1, oooo2;

double iko0, ikB10, ikB20, ikB30, ikB40, iko1, ikB11, ikB21, ikB31, ikB41;
 
/* List of variables and paramaters (this code uses all global variables) */
 
	/* Creation of Data File */
	FILE *ap;
	FILE *fmaxs;
	FILE *fpara;
	void prttofile ();	
	int printdata;	
	int printval;	
 
	/* Cell Geometry */
	const double l = 0.01;       // Length of the cell (cm)
	const double a = 0.0011;     // Radius of the cell (cm)
	const double pi = 3.141592;  // Pi
	double vcell;   // Cell volume (uL)
	double ageo;    // Geometric membrane area (cm^2)
	double acap;    // Capacitive membrane area (cm^2)
	double vmyo;    // Myoplasm volume (uL)
	double vmito;   // Mitochondria volume (uL)
	double vsr;     // SR volume (uL)
	double vnsr;    // NSR volume (uL)
	double vjsr;    // JSR volume (uL)
	double vcleft;  // Cleft volume (uL)
	
	/* Voltage */
	double v;       // Membrane voltage (mV)
	double vnew;    // New Voltage (mV)
	double dvdt;    // Change in Voltage / Change in Time (mV/ms)
	double dvdtnew; // New dv/dt (mV/ms)
	double flag; // Flag condition to test for dvdtmax
	
	/* Time Step */
	double dt;      // Time step (ms)
	double t;       // Time (ms)
	double udt;     // Universal Time Step
	int utsc;       // Universal Time Step Counter
	int nxstep;     // Interval Between Calculating Ion Currents
	int steps;      // Number of Steps
	int increment;  // Loop Control Variable

	/* Action Potential Duration and Max. Info */
	double vmax[beats];           // Max. Voltage (mV)
	double dvdtmax[beats];        // Max. dv/dt (mV/ms)
	double apd[beats];            // Action Potential Duration
	double toneapd[beats];        // Time of dv/dt Max.
	double ttwoapd[beats];        // Time of 90% Repolarization
	double rmbp[beats];           // Resting Membrane Potential
	double nair[beats];           // Intracellular Na At Rest
	double cair[beats];           // Intracellular Ca At Rest
	double caimax[beats];	      // Peak Intracellular Ca
	
	int i;                        // Stimulation Counter
	
	/* Total Current and Stimulus */
	double st;       // Constant Stimulus (uA/cm^2)
	double tstim;    // Time Stimulus is Applied (ms)
	double stimtime; // Time period during which stimulus is applied (ms)
	double it;       // Total current (uA/cm^2)
 
	/* Terms for Solution of Conductance and Reversal Potential */
	const double R = 8314;      // Universal Gas Constant (J/kmol*K)
	const double frdy = 96485;  // Faraday's Constant (C/mol)
	const double temp = 310;    // Temperature (K)
 
	/* Ion Valences */
	const double zna = 1;  // Na valence
	const double zk = 1;   // K valence
	const double zca = 2;  // Ca valence
 
	/* Ion Concentrations */
	double nai;    // Intracellular Na Concentration (mM)
	double nao;    // Extracellular Na Concentration (mM)
	double nabm;   // Bulk Medium Na Concentration (mM)
	double dnao;   // Change in Cleft Na Concentration (mM)
	double ki;     // Intracellular K Concentration (mM)
 	double ko;     // Extracellular K Concentration (mM)
	double kbm;    // Bulk Medium K Concentration (mM)
 	double dko;    // Change in Cleft K Concentration (mM)
	double cai;    // Intracellular Ca Concentration (mM)
	double cao;    // Extracellular Ca Concentration (mM)
	double cabm;   // Bulk Medium Ca Concentration (mM)
	double dcao;   // Change in Cleft Ca Concentration (mM)
	double cmdn;   // Calmodulin Buffered Ca Concentration (mM)
	double trpn;   // Troponin Buffered Ca Concentration (mM)
	double nsr;    // NSR Ca Concentration (mM)
	double jsr;    // JSR Ca Concentration (mM)
	double csqn;   // Calsequestrin Buffered Ca Concentration (mM)
	const double taudiff = 1000; // Diffusion Constant for Ion Movement from Bulk Medium to Cleft Space
 
	/* Myoplasmic Na Ion Concentration Changes */
	double naiont;  // Total Na Ion Flow (uA/uF)
	double dnai;    // Change in Intracellular Na Concentration (mM)
 
	/* Myoplasmic K Ion Concentration Changes */
	double kiont; // Total K Ion Flow (uA/uF)
	double dki;   // Change in Intracellular K Concentration (mM)
 
	/* NSR Ca Ion Concentration Changes */
	double dnsr;   // Change in [Ca] in the NSR (mM)
	double iup;    // Ca uptake from myo. to NSR (mM/ms)
	double ileak;  // Ca leakage from NSR to myo. (mM/ms)
	double kleak;  // Rate constant of Ca leakage from NSR (ms^-1)
	const double kmup = 0.00092;    // Half-saturation concentration of iup (mM)
	const double iupbar = 0.00875;  // Max. current through iup channel (mM/ms)
	const double nsrbar = 15;       // Max. [Ca] in NSR (mM) 
	
	/* JSR Ca Ion Concentration Changes */
	double djsr;                    // Change in [Ca] in the JSR (mM)
	const double tauon = 0.5;         // Time constant of activation of Ca release from JSR (ms)
	const double tauoff = 0.5;        // Time constant of deactivation of Ca release from JSR (ms)
	double tcicr;                   // t=0 at time of CICR (ms)
	double irelcicr;                // Ca release from JSR to myo. due to CICR (mM/ms)
	const double csqnth = 8.75;     // 8.75 Threshold for release of Ca from CSQN due to JSR overload (mM)
	const double gmaxrel = 150;     // Max. rate constant of Ca release from JSR due to overload (ms^-1)
	double grelbarjsrol;            // Rate constant of Ca release from JSR due to overload (ms^-1)
	double greljsrol;               // Rate constant of Ca release from JSR due to CICR (ms^-1)
	double tjsrol;                  // t=0 at time of JSR overload (ms)
	double ireljsrol;               // Ca release from JSR to myo. due to JSR overload (mM/ms)
	const double csqnbar = 10;      // Max. [Ca] buffered in CSQN (mM)
	const double kmcsqn = 0.8;      // Equalibrium constant of buffering for CSQN (mM)
	double bjsr;                    // b Variable for analytical computation of [Ca] in JSR (mM)
	double cjsr;                    // c Variable for analytical computation of [Ca] in JSR (mM)
	double on;                      // Time constant of activation of Ca release from JSR (ms)
	double off;                     // Time constant of deactivation of Ca release from JSR (ms)
	double magrel;                  // Magnitude of Ca release
	double dcaiont;                 // Rate of change of Ca entry
	double dcaiontnew;              // New rate of change of Ca entry
	double caiontold;               // Old rate of change of Ca entry
 
	/* Translocation of Ca Ions from NSR to JSR */
	double itr;                // Translocation current of Ca ions from NSR to JSR (mM/ms)
	const double tautr = 180;  // Time constant of Ca transfer from NSR to JSR (ms)
	
	/* Myoplasmic Ca Ion Concentration Changes */
	double caiont;  // Total Ca Ion Flow (uA/uF)
	double dcai;    // Change in myoplasmic Ca concentration (mM)
	double catotal; // Total myoplasmic Ca concentration (mM)
	double bmyo;    // b Variable for analytical computation of [Ca] in myoplasm (mM)
	double cmyo;    // c Variable for analytical computation of [Ca] in myoplasm (mM)
	double dmyo;    // d Variable for analytical computation of [Ca] in myoplasm (mM)
	double gpig;    // Tribute to all the guinea pigs killed for the advancement of knowledge
	const double cmdnbar = 0.050;   // Max. [Ca] buffered in CMDN (mM)
	const double trpnbar = 0.070;   // Max. [Ca] buffered in TRPN (mM)
	const double kmcmdn = 0.00238;  // Equalibrium constant of buffering for CMDN (mM)
	const double kmtrpn = 0.0005;   // Equalibrium constant of buffering for TRPN (mM)
 
	/* Fast Sodium Current (time dependant) */
	double ina;    // Fast Na Current (uA/uF)
	double gna;    // Max. Conductance of the Na Channel (mS/uF)
	double ena;    // Reversal Potential of Na (mV)
	double am;     // Na alpha-m rate constant (ms^-1)
	double bm;     // Na beta-m rate constant (ms^-1)
	double ah;     // Na alpha-h rate constant (ms^-1)
	double bh;     // Na beta-h rate constant (ms^-1)
	double aj;     // Na alpha-j rate constant (ms^-1)
	double bj;     // Na beta-j rate constant (ms^-1)
	double mtau;      // Na activation
	double htau;      // Na inactivation
	double jtau;      // Na inactivation
	double mss;      // Na activation
	double hss;      // Na inactivation
	double jss;      // Na inactivation
	double m;      // Na activation
	double h;      // Na inactivation
	double j;      // Na inactivation
 
	/* Current through L-type Ca Channel */
	double ilca;    // Ca current through L-type Ca channel (uA/uF)
	double ilcana;  // Na current through L-type Ca channel (uA/uF)
	double ilcak ;  // K current through L-type Ca channel (uA/uF)
	double ilcatot; // Total current through the L-type Ca channel (uA/uF)
	double ibarca;  // Max. Ca current through Ca channel (uA/uF)
	double ibarna;  // Max. Na current through Ca channel (uA/uF)
	double ibark;   // Max. K current through Ca channel (uA/uF)
	double d;       // Voltage dependant activation gate
	double dss;     // Steady-state value of activation gate d 
	double taud;    // Time constant of gate d (ms^-1)
	double f;       // Voltage dependant inactivation gate
	double fss;     // Steady-state value of inactivation gate f
	double tauf;    // Time constant of gate f (ms^-1)
	double fca;     // Ca dependant inactivation gate
	const double kmca = 0.0006;     // Half-saturation concentration of Ca channel (mM)
	const double atpi = 3; // Intracellular ATP concentration (mM)
	
	const double pca = 0.00054;     // Permiability of membrane to Ca (cm/s)
	const double gacai = 1;         // Activity coefficient of Ca
	const double gacao = 0.341;     // Activity coefficient of Ca
	const double pna = 0.000000675; // Permiability of membrane to Na (cm/s)
	const double ganai = 0.75;      // Activity coefficient of Na
	const double ganao = 0.75;      // Activity coefficient of Na
	const double pk = 0.000000193;  // Permiability of membrane to K (cm/s)
	const double gaki = 0.75;       // Activity coefficient of K
	const double gako = 0.75;       // Activity coefficient of K
 
	/* Current through T-type Ca Channel */
	double icat;    // Ca current through T-type Ca channel (uA/uF)
	double gcat;    // Max. Conductance of the T-type Ca channel (mS/uF)
	double eca;     // Reversal Potential of the T-type Ca channel (mV)
	double b;       // Voltage dependant activation gate
	double bss;     // Steady-state value of activation gate b 
	double taub;    // Time constant of gate b (ms^-1)
	double g;       // Voltage dependant inactivation gate
	double gss;     // Steady-state value of inactivation gate g
	double taug;    // Time constant of gate g (ms^-1)
 
	/* Rapidly Activating Potassium Current */
	double ikr;   // Rapidly Activating K Current (uA/uF)
	double gkr;   // Channel Conductance of Rapidly Activating K Current (mS/uF)
	double ekr;   // Reversal Potential of Rapidly Activating K Current (mV)
	double xr;    // Rapidly Activating K time-dependant activation
	double xrss;  // Steady-state value of inactivation gate xr
	double tauxr; // Time constant of gate xr (ms^-1)
	double r;     // K time-independant inactivation
	
	/* Slowly Activating Potassium Current */
	double iks;    // Slowly Activating K Current (uA/uF)
	double gks;    // Channel Conductance of Slowly Activating K Current (mS/uF)
	double eks;    // Reversal Potential of Slowly Activating K Current (mV)
	double xs1;    // Slowly Activating K time-dependant activation
	double xs1ss;  // Steady-state value of inactivation gate xs1
	double tauxs1; // Time constant of gate xs1 (ms^-1)
	double xs2;    // Slowly Activating K time-dependant activation
	double xs2ss;  // Steady-state value of inactivation gate xs2
	double tauxs2; // Time constant of gate xs2 (ms^-1)
	const double prnak = 0.01833;  // Na/K Permiability Ratio
	
	/* Potassium Current (time-independant) */
	double iki;    // Time-independant K current (uA/uF)
	double gki;    // Channel Conductance of Time Independant K Current (mS/uF)
	double eki;    // Reversal Potential of Time Independant K Current (mV)
	double aki;    // K alpha-ki rate constant (ms^-1)
	double bki;    // K beta-ki rate constant (ms^-1)
	double kin;    // K inactivation
 
	/* Plateau Potassium Current */
	double ikp;    // Plateau K current (uA/uF)
	double gkp;    // Channel Conductance of Plateau K Current (mS/uF)
	double ekp;    // Reversal Potential of Plateau K Current (mV)
	double kp;     // K plateau factor	
 
	/* Na-Activated K Channel */
	double ikna;   // Na activated K channel
	double pona;   // Open probability dependant on Nai
	double pov;    // Open probability dependant on Voltage
	double ekna;   // Reversal potential
	const double gkna = 0.12848;   // Maximum conductance (mS/uF)
	const double nkna = 2.8;       // Hill coefficient for Na dependance
	const double kdkna = 66;       // Dissociation constant for Na dependance(mM)
	
	/* ATP-Sensitive K Channel */
	double ikatp;    // ATP-sensitive K current (uA/uF)
	double ekatp;    // K reversal potential (mV)
	double gkbaratp; // Conductance of the ATP-sensitive K channel (mS/uF)
	double gkatp;    // Maximum conductance of the ATP-sensitive K channel (mS/uF)
	double patp;     // Percentage availibility of open channels
	const double natp = 0.24;          // K dependence of ATP-sensitive K current
	const double nicholsarea = 0.005;  // Nichol's ares (cm^2)
    const double hatp = 2;             // Hill coefficient
	const double katp = 0.250;         // Half-maximal saturation point of ATP-sensitive K current (mM)
	
	/* Ito Transient Outward Current (Dumaine et al. Circ Res 1999;85:803-809) */
	double ito;	      // Transient outward current
	double gitodv;	  // Maximum conductance of Ito
	double ekdv;	  // Reversal Potential of Ito
	double rvdv;      // Time independant voltage dependence of Ito
	double zdv;       // Ito activation
	double azdv;      // Ito alpha-z rate constant
	double bzdv;      // Ito beta-z rate constant
	double tauzdv;	  // Time constant of z gate
	double zssdv;     // Steady-state value of z gate
	double ydv;       // Ito inactivation
	double aydv;      // Ito alpha-y rate constant
	double bydv;      // Ito beta-y rate constant
	double tauydv;	  // Time constant of y gate
	double yssdv;     // Steady-state value of y gate
		
	/* Sodium-Calcium Exchanger V-S */
	double inaca;               // NaCa exchanger current (uA/uF)
	const double c1 = .00025;   // Scaling factor for inaca (uA/uF)
	const double c2 = 0.0001;   // Half-saturation concentration of NaCa exhanger (mM)
	const double gammas = .15;  // Position of energy barrier controlling voltage dependance of inaca
 
	/* Sodium-Potassium Pump */
	double inak;    // NaK pump current (uA/uF)
	double fnak;    // Voltage-dependance parameter of inak
	double sigma;   // [Na]o dependance factor of fnak
	const double ibarnak = 2.25;   // Max. current through Na-K pump (uA/uF)
	const double kmnai = 10;       // Half-saturation concentration of NaK pump (mM)
	const double kmko = 1.5;       // Half-saturation concentration of NaK pump (mM)
	
	/* Nonspecific Ca-activated Current */
	double insna;     // Non-specific Na current (uA/uF)
	double insk;      // Non-specific K current (uA/uF)
	double ibarnsna;  // Max. Na current through NSCa channel (uA/uF)
	double ibarnsk;   // Max. K current through NSCa channel (uA/uF)
	const double pnsca = 0.000000175;  // Permiability of channel to Na and K (cm/s)
	const double kmnsca = 0.0012;      // Half-saturation concentration of NSCa channel (mM)
 
	/* Sarcolemmal Ca Pump */
	double ipca;                 // Sarcolemmal Ca pump current (uA/uF)
	const double ibarpca = 1.15; // Max. Ca current through sarcolemmal Ca pump (uA/uF)
	const double kmpca = 0.0005; // Half-saturation concentration of sarcolemmal Ca pump (mM)
	
	/* Ca Background Current */
	double icab;  // Ca background current (uA/uF)
	double gcab;  // Max. conductance of Ca background (mS/uF)
	double ecan;  // Nernst potential for Ca (mV)
 
	/* Na Background Current */
	double inab;  // Na background current (uA/uF)
	double gnab;  // Max. conductance of Na background (mS/uF)
	double enan;  // Nernst potential for Na (mV)
 
	/* Ion Current Functions */	
	void comp_ina ();    // Calculates Fast Na Current
	void comp_ical ();   // Calculates Currents through L-Type Ca Channel
	void comp_icat ();   // Calculates Currents through T-Type Ca Channel
	void comp_ikr ();    // Calculates Rapidly Activating K Current
	void comp_iks ();    // Calculates Slowly Activating K Current
	void comp_iki ();    // Calculates Time-Independant K Current
	void comp_ikp ();    // Calculates Plateau K Current
	void comp_ikna ();   // Calculates Na-activated K Current
	void comp_ikatp ();  // Calculates ATP-Sensitive K Current
	void comp_ito ();    // Calculates Transient Outward Current
	void comp_inaca ();  // Calculates Na-Ca Exchanger Current
	void comp_inak ();   // Calculates Na-K Pump Current
	void comp_insca ();  // Calculates Non-Specific ca-Activated Current
	void comp_ipca ();   // Calculates Sarcolemmal Ca Pump Current
	void comp_icab ();   // Calculates Ca Background Current
	void comp_inab ();   // Calculates Na Background Current
	void comp_it ();     // Calculates Total Current
 
	/* Ion Concentration Functions */
	void conc_nai ();    // Calculates new myoplasmic Na ion concentration
	void conc_ki ();     // Calculates new myoplasmic K ion concentration
	void conc_nsr ();    // Calculates new NSR Ca ion concentration
	void conc_jsr ();    // Calculates new JSR Ca ion concentration
	void calc_itr ();    // Calculates Translocation of Ca from NSR to JSR
	void conc_cai ();    // Calculates new myoplasmic Ca ion concentration
	void conc_cleft ();  // Calculates new cleft ion concentrations
 
int main ()
{	
	/* Opening of Datafiles */
	ap = fopen("ap", "w");
	fpara = fopen("fpara", "w");
	fmaxs = fopen("fmaxs", "w");
	printdata = 60;
	
	/* Cell Geometry */
	vcell = 1000*pi*a*a*l;     //   3.801e-5 uL
	ageo = 2*pi*a*a+2*pi*a*l;  //   7.671e-5 cm^2
	acap = ageo*2;             //   1.534e-4 cm^2
	vmyo = vcell*0.68;
	vmito = vcell*0.26;
	vsr = vcell*0.06;
	vnsr = vcell*0.0552;
	vjsr = vcell*0.0048;
	vcleft = vcell*0.12/0.88;
	
	/* Time Loop Conditions */
	t = 0.0;           // Time (ms)
	udt = 0.002;       // Time step (ms)  
	steps = (bcl*beats)/udt; // Number of ms
	st = -80.0;        // Stimulus  
	tstim = 10.0;      // Time to begin stimulus
	stimtime = 10.0;   // Initial Condition for Stimulus
	v = -88.654973;    // Initial Voltage (mv)
 
	/* Beginning Ion Concentrations */
	nai = 15;       // Initial Intracellular Na (mM)
	nao = 140;      // Initial Extracellular Na (mM)
	nabm = 140;     // Initial Bulk Medium Na (mM)
	ki = 136.89149; // Initial Intracellular K (mM)
	ko = 4.5;       // 5.4 Initial Extracellular K (mM)//don
	kbm = 4.5;      // Initial Bulk Medium K (mM)
	cai = 0.000079; // Initial Intracellular Ca (mM)
	cao = 1.8;      // Initial Extracellular Ca (mM)
	cabm = 1.8;     // Initial Bulk Medium Ca (mM)
 
	/* Initial Gate Conditions */
	m = 0.000838;
	h = 0.993336;
	j = 0.995484;
	d = 0.000003;
	f = 0.999745;
	xs1 = 0.004503;
	xs2 = 0.004503;
	xr = 0.000129;
	b = 0.000994;
	g = 0.994041;
	zdv = 0.0120892;
	ydv = 0.999978;
	iko0=0, ikB10=0, ikB20=0, ikB30=0, ikB40=0.32, iko1=0, ikB11=0, ikB21, ikB31=0,ikB41=0.68;
	
	/* Initial Conditions */
	grelbarjsrol = 0;
	tjsrol = 1000;
	tcicr = 1000;
	jsr = 1.179991;
	nsr = 1.179991;
	trpn = 0.0143923;
	cmdn = 0.00257849;
	csqn = 6.97978;
	flag = 0;
	dt = udt;
	utsc = 50;
	dcaiont = 0;
	i=-1;
 
	/* Beginning of Time Loop */
	for (increment = 0; increment < steps; increment++)
	{
		if(abs(dvdt)<0.25 && v<0)
			{nxstep = 5;}
		else
			{nxstep = 5;}
 
	
		
		if(utsc>=nxstep || dvdt>5 || irelcicr>0.01 || (t>=(tstim-udt) && t<=(tstim+udt)) || (stimtime>=0 && stimtime<0.5))
		{
 
		comp_ina ();
		comp_ical ();
		comp_icat ();
		comp_ikr ();
		comp_iks ();
		comp_iki ();
		comp_ikp ();
		comp_ikna ();
		comp_ikatp ();
		comp_ito ();
		comp_inaca ();
		comp_inak ();
		comp_insca ();
		comp_ipca ();
		comp_icab ();
		comp_inab ();
		comp_it ();
		
		conc_nai ();
		conc_ki ();
		calc_itr ();
		conc_jsr ();
		conc_nsr ();
		conc_cai ();
 
		stimtime = stimtime+dt;
	
		//conc_cleft ();  /* Cleft Space disabled, if you want to use cleft space, make sure the initial conditions of ion concentrations in the bulk medium are the same as the extracellular concentrations */
		
		utsc = 0;
		dt = 0;
 
		}
 
			if(dvdt>3 || irelcicr>.01)
			{printval = 50;}
			else
			{printval = 750;}
			
			if(printdata>=printval) 
			{prttofile();	
			printdata = 0;}
			
			printdata = printdata+1;
 
	
	vnew = v-it*udt;
	dvdtnew = (vnew-v)/udt;
 
	
	if(i>=0)
		{
		if (vnew>vmax[i])
			vmax[i] = vnew;
		if (cai>caimax[i])
			caimax[i] = cai;
		if (dvdtnew>dvdtmax[i])
			{dvdtmax[i] = dvdtnew;
			toneapd[i] = t;}
		if (vnew>=(vmax[i]-0.9*(vmax[i]-rmbp[i])))
			ttwoapd[i] = t;
		}
 
	if(csqn>=csqnth && tjsrol>50)
		{grelbarjsrol = 4;
		tjsrol = 0;
		cout << "Spontaneous Release occured at time " << t << endl;
		}	
  
	v = vnew;
	dvdt = dvdtnew;
	caiontold = caiont;
	dcaiont = dcaiontnew;
	
	dt = dt+udt;
	utsc = utsc+1;
	
	t = t+udt;
 
	}
 
	fprintf(fpara,"%.3f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\t%.16f\n", t, v, nai, ki, cai, jsr, nsr, nao, ko, cao, m, h, j, d, f, xs1, xs2, xr, b, g, tcicr, flag);
	for(i=0;i<beats;i++)
	{apd[i] = ttwoapd[i]-toneapd[i];
       	fprintf(fmaxs, "%i\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", i, vmax[i], dvdtmax[i], apd[i], toneapd[i], ttwoapd[i], nair[i], cair[i], caimax[i], rmbp[i]);}
 
	cout << "\a\a";
 
	return (1);
}
 
/********************************************************/
 
/* Functions that describe the currents begin here */
  
void comp_ina ()
{
	gna = 16;
	ena = ((R*temp)/frdy)*log(nao/nai);
 
	am = 0.32*(v+47.13)/(1-exp(-0.1*(v+47.13)));
	bm = 0.08*exp(-v/11);
	
	if (v < -40)
	{ah = 0.135*exp((80+v)/-6.8);
	bh = 3.56*exp(0.079*v)+310000*exp(0.35*v);
	aj = (-127140*exp(0.2444*v)-0.00003474*exp(-0.04391*v))*((v+37.78)/(1+exp(0.311*(v+79.23))));
	bj = (0.1212*exp(-0.01052*v))/(1+exp(-0.1378*(v+40.14)));}
 
	else
	{ah = 0;
	bh = 1/(0.13*(1+exp((v+10.66)/-11.1)));
	aj = 0;
	bj = (0.3*exp(-0.0000002535*v))/(1+exp(-0.1*(v+32)));}
	
	mtau = 1/(am+bm);
	htau = 1/(ah+bh);
	jtau = 1/(aj+bj);
 
	mss = am*mtau;
	hss = ah*htau;
	jss = aj*jtau;
 
       	m = mss-(mss-m)*exp(-dt/mtau);
       	h = hss-(hss-h)*exp(-dt/htau);
       	j = jss-(jss-j)*exp(-dt/jtau);
 	
	ina = gna*m*m*m*h*j*(v-ena);
}
 
void comp_ical ()
{
	dss = 1/(1+exp(-(v+10)/6.24));
	taud = dss*(1-exp(-(v+10)/6.24))/(0.035*(v+10));
 
	fss = (1/(1+exp((v+32)/8)))+(0.6/(1+exp((50-v)/20)));
	tauf = 1/(0.0197*exp(-pow(0.0337*(v+10),2))+0.02);
 
	d = dss-(dss-d)*exp(-dt/taud);
	f = fss-(fss-f)*exp(-dt/tauf);
	
	ibarca = pca*zca*zca*((v*frdy*frdy)/(R*temp))*((gacai*cai*exp((zca*v*frdy)/(R*temp))-gacao*cao)/(exp((zca*v*frdy)/(R*temp))-1));
	ibarna = pna*zna*zna*((v*frdy*frdy)/(R*temp))*((ganai*nai*exp((zna*v*frdy)/(R*temp))-ganao*nao)/(exp((zna*v*frdy)/(R*temp))-1));
	ibark = pk*zk*zk*((v*frdy*frdy)/(R*temp))*((gaki*ki*exp((zk*v*frdy)/(R*temp))-gako*ko)/(exp((zk*v*frdy)/(R*temp))-1));
		
	fca = 1/(1+cai/kmca);
	
	ilca = d*f*fca*ibarca;
	ilcana = d*f*fca*ibarna;
	ilcak = d*f*fca*ibark;
	
	ilcatot = ilca+ilcana+ilcak;
}
 
void comp_icat ()
{
	bss = 1/(1+exp(-(v+14.0)/10.8));
	taub = 3.7+6.1/(1+exp((v+25.0)/4.5));
 
	gss = 1/(1+exp((v+60.0)/5.6));
	if (v<=0)
	taug = -0.875*v+12.0;
	else
	taug = 12.0;
 
	b = bss-(bss-b)*exp(-dt/taub);
	g = gss-(gss-g)*exp(-dt/taug);
	
	gcat = 0.05;
	eca = (R*temp/(2*frdy))*log(cao/cai);
	
	icat = gcat*b*b*g*(v-eca);
}

/* Markovian model for IKr */ 
void comp_ikr ()
{
if (t==0)
{cc3=1, cc2=0, cc1=0, oo=0, ii=0;}

gkr = 1.35e-2*pow((ko),0.59); // 0.2 WU

double aa = 1.31e-2*exp(1.48*v*frdy/(R*temp)) ; //aa=alpha
double bb = 2.9e-3*exp(-5.77e-1*v*frdy/(R*temp)); //bb=beta
double a1 = 2.17; // alpha1
double b1 = 1.08; //beta1
double a2 = 3.02e-2*exp(1.48*v*frdy/(R*temp));// alpha2
double b2 = 2.90e-3*exp(-9.78e-1*v*frdy/(R*temp)); //beta2
double bi = 8.20e-1*exp(5.04e-1*v*frdy/(R*temp))*pow((4.5/ko),3); //betai
double ai = 5.45e-1*exp(-8.04e-1*v*frdy/(R*temp))*(4.5/ko); //alphai
double u = (ai*b2)/(bi);
	
double dcc3=(bb*cc2-aa*cc3)*dt;
double dcc2=(b1*cc1+aa*cc3-(a1+bb)*cc2)*dt;
double dcc1=(a1*cc2+b2*oo+u*ii-(b1+2*a2)*cc1)*dt;
double doo=(ai*ii+a2*cc1-(bi+b2)*oo)*dt;
double dii=(a2*cc1+bi*oo-(u+ai)*ii)*dt;

 cc3=cc3+dcc3;
 cc2=cc2+dcc2;
 cc1=cc1+dcc1;
 oo=oo+doo;
 ii=ii+dii;

double Okr=(1-cc1-cc2-cc3-ii);

	ekr = ((R*temp)/frdy)*log(ko/ki);
 	ikr = gkr*Okr*(v-ekr);
}
 
/* Markovian model for IKs */ 
void comp_iks ()
{
if (t==0)
{cccc1=1.0, cccc2=0, cccc3=0, cccc4=0, cccc5=0, cccc6=0, cccc7=0, cccc8=0, cccc9=0, cccc10=0, cccc11=0, cccc12=0, cccc13=0, cccc14=0, cccc15=0, oooo1=0, oooo2=0;}

gks = 0.779*(1+0.6/(1+pow((3.8e-5/cai),1.4))); //7/23 for M cell Cai: mM

double sa = 3.98e-4*exp( 3.61e-1*v*frdy/(R*temp));
double sb = 5.74e-5*exp(-9.23e-2*v*frdy/(R*temp));
double sr = 3.41e-3*exp( 8.68e-1*v*frdy/(R*temp));
double sd = 1.2e-3*exp(  -3.3e-1*v*frdy/(R*temp));
double stheta = 6.47e-3;
double seta = 1.25e-2*exp(-4.81e-1*v*frdy/(R*temp));
double spsi = 6.33e-3*exp(    1.27*v*frdy/(R*temp));
double somega = 4.91e-3*exp(-6.79e-1*frdy/(R*temp));

double dcccc1 =  (cccc2*sb-cccc1*4*sa)*dt;
double dcccc2 =  (cccc1*4*sa+cccc3*2*sb+cccc6*sd-cccc2*(sb+3*sa+sr))*dt;
double dcccc3 =  (cccc2*3*sa+cccc4*3*sb+cccc7*sd-cccc3*(2*sb+2*sa+2*sr))*dt;
double dcccc4 =  (cccc3*2*sa+cccc5*4*sb+cccc8*sd-cccc4*(3*sb+sa+3*sr))*dt;
double dcccc5 =  (cccc4*sa+cccc9*sd-cccc5*(4*sb+4*sr))*dt;
double dcccc6 =  (cccc2*sr+cccc7*sb-cccc6*(sd+3*sa))*dt;
double dcccc7 =  (cccc6*3*sa+cccc8*2*sb+cccc3*2*sr+cccc10*2*sd-cccc7*(sb+2*sa+sd+sr))*dt;
double dcccc8 =  (cccc7*2*sa+cccc9*3*sb+cccc4*3*sr+cccc11*2*sd-cccc8*(2*sb+sa+sd+2*sr))*dt; 
double dcccc9 =  (cccc8*sa+cccc5*4*sr+cccc12*2*sd-cccc9*(3*sb+sd+3*sr))*dt;
double dcccc10 = (cccc11*sb+cccc7*sr-cccc10*(2*sa+2*sd))*dt; 
double dcccc11 = (cccc10*2*sa+cccc12*2*sb+cccc8*2*sr+cccc13*3*sd-cccc11*(1*sb+sa+2*sd+sr))*dt;
double dcccc12 = (cccc11*sa+cccc9*3*sr+cccc14*3*sd-cccc12*(2*sb+2*sd+2*sr))*dt;
double dcccc13 = (cccc11*sr+cccc14*sb-cccc13*(3*sd+sa))*dt;
double dcccc14 = (cccc13*sa+cccc12*2*sr+cccc15*4*sd-cccc14*(sb+3*sd+sr))*dt;
double dcccc15 = (cccc14*sr+oooo1*seta-cccc15*(4*sd+stheta))*dt;
double doooo1 =  (cccc15*stheta+oooo2*somega-oooo1*(seta+spsi))*dt;
double doooo2 =  (oooo1*spsi-oooo2*somega)*dt;

cccc1=cccc1+dcccc1;
cccc2=cccc2+dcccc2;
cccc3=cccc3+dcccc3;
cccc4=cccc4+dcccc4;
cccc5=cccc5+dcccc5;
cccc6=cccc6+dcccc6;
cccc7=cccc7+dcccc7;
cccc8=cccc8+dcccc8;
cccc9=cccc9+dcccc9;
cccc10=cccc10+dcccc10;
cccc11=cccc11+dcccc11;
cccc12=cccc12+dcccc12;
cccc13=cccc13+dcccc13;
cccc14=cccc14+dcccc14;
cccc15=cccc15+dcccc15;
oooo1=oooo1+doooo1;
oooo2=oooo2+doooo2;

double Oks=oooo1+oooo2;

	eks = ((R*temp)/frdy)*log((ko+prnak*nao)/(ki+prnak*nai));
 	
	iks = gks*Oks*(v-eks);
}
 
void comp_iki ()
{
	gki = 0.75*(sqrt(ko/5.4));
	eki = ((R*temp)/frdy)*log(ko/ki);
 
	aki = 1.02/(1+exp(0.2385*(v-eki-59.215)));
	bki = (0.49124*exp(0.08032*(v-eki+5.476))+exp(0.06175*(v-eki-594.31)))/(1+exp(-0.5143*(v-eki+4.753)));
	
	kin = aki/(aki+bki);
	
	iki = gki*kin*(v-eki);
}
 
void comp_ikp ()
{
	gkp = 0.00552;
	ekp = eki;
 
	kp = 1/(1+exp((7.488-v)/5.98));	
 
	ikp = gkp*kp*(v-ekp);
}
 
void comp_ikna ()
{
	ekna = ((R*temp)/frdy)*log(ko/ki);
	pona = 0.85/(1+pow((kdkna/nai),2.8));
	pov = 0.8-0.65/(1+exp((v+125)/15));
	ikna = gkna*pona*pov*(v-ekna);	
}
 
void comp_ikatp ()
{
/* Note: If you wish to use this current in your simulations, there are additional       */
/* changes which must be made to the code as detailed in Cardiovasc Res 1997;35:256-272  */
 
	ekatp = ((R*temp)/frdy)*log(ko/ki);
	gkatp = 1*0.000195/nicholsarea;//20fold
	patp = 1/(1+(pow((atpi/katp),hatp)));
	gkbaratp = gkatp*patp*(pow((ko/4),natp));
 
	ikatp = gkbaratp*(v-ekatp);	
}
 
void comp_ito ()
{
	gitodv = 0.5*0.01;
	ekdv = ((R*temp)/frdy)*log((ko)/(ki));
	rvdv = exp(v/100);
	
	azdv = (10*exp((v-40)/25))/(1+exp((v-40)/25));
	bzdv = (10*exp(-(v+90)/25))/(1+exp(-(v+90)/25));
	tauzdv = 1/(azdv+bzdv);
	zssdv = azdv/(azdv+bzdv);
	zdv = zssdv-(zssdv-zdv)*exp(-dt/tauzdv);
 
	aydv = 0.015/(1+exp((v+60)/5));
	bydv = (0.1*exp((v+25)/5))/(1+exp((v+25)/5));
	tauydv = 1/(aydv+bydv);
	yssdv = aydv/(aydv+bydv);
	ydv = yssdv-(yssdv-ydv)*exp(-dt/tauydv);
	
	ito = gitodv*zdv*zdv*zdv*ydv*rvdv*(v-ekdv);
}
 
void comp_inaca ()
{
	inaca = c1*exp((gammas-1)*v*frdy/(R*temp))*((exp(v*frdy/(R*temp))*nai*nai*nai*cao-nao*nao*nao*cai)/(1+c2*exp((gammas-1)*v*frdy/(R*temp))*(exp(v*frdy/(R*temp))*nai*nai*nai*cao+nao*nao*nao*cai)));
}
 
void comp_inak ()
{
	sigma = (exp(nao/67.3)-1)/7;
 
	fnak = 1/(1+0.1245*exp((-0.1*v*frdy)/(R*temp))+0.0365*sigma*exp((-v*frdy)/(R*temp)));
 
	inak = ibarnak*fnak*(1/(1+pow(kmnai/nai,2)))*(ko/(ko+kmko));
}
 
void comp_insca ()
{
	ibarnsna = pnsca*zna*zna*((v*frdy*frdy)/(R*temp))*((ganai*nai*exp((zna*v*frdy)/(R*temp))-ganao*nao)/(exp((zna*v*frdy)/(R*temp))-1));
	ibarnsk = pnsca*zk*zk*((v*frdy*frdy)/(R*temp))*((gaki*ki*exp((zk*v*frdy)/(R*temp))-gako*ko)/(exp((zk*v*frdy)/(R*temp))-1));
 
	insna = ibarnsna/(1+pow(kmnsca/cai,3));
	insk = ibarnsk/(1+pow(kmnsca/cai,3));
} 
 
void comp_ipca ()
{
	ipca = (ibarpca*cai)/(kmpca+cai);	
}
 
void comp_icab ()
{
	gcab = 0.003016;
	ecan = ((R*temp)/(2*frdy))*log(cao/cai);
		
	icab = gcab*(v-ecan);
}
 
void comp_inab ()
{
	gnab = 0.004;
	enan = ena;
	
	inab = gnab*(v-enan);
}
 
/* Total sum of currents is calculated here, if the time is between stimtime = 0 and stimtime = 0.5, a stimulus is applied */
 
 
void comp_it ()
{
	naiont = ina+inab+ilcana+insna+3*inak+3*inaca;
	kiont = ikr+iks+iki+ikp+ilcak+insk-2*inak+ito+ikna+ikatp;
	caiont = ilca+icab+ipca-2*inaca+icat;
	
	if (dvdtnew > 10 && tcicr > 10 && flag == 1)
	  {flag = 0;}
 
	if (t>=tstim && t<(tstim+dt))
	{stimtime = 0;
	i = i+1;

	tstim = tstim+bcl;
	rmbp[i]=v;
	nair[i] = nai;
	cair[i] = cai;}
	
	if(stimtime>=0 && stimtime<=0.5)
	{it = st+naiont+kiont+caiont;}
	else
	{it = naiont+kiont+caiont;}
}
 
/* Functions that calculate intracellular ion concentrations begins here */
 
void conc_nai ()
{
// The units of dnai is in mM.  Note that naiont should be multiplied by the
// cell capacitance to get the correct units.  Since cell capacitance = 1 uF/cm^2,
// it doesn't explicitly appear in the equation below.
// This holds true for the calculation of dki and dcai. */
 
	dnai = -dt*(naiont*acap)/(vmyo*zna*frdy);
	nai = dnai + nai;
}
 
void conc_ki ()
{
	if(stimtime>=0 && stimtime<=0.5)
	{dki = -dt*((kiont+st)*acap)/(vmyo*zk*frdy);}
	else
	{dki = -dt*(kiont*acap)/(vmyo*zk*frdy);}
 
	ki = dki + ki;
}
 
void calc_itr ()
{
	itr = (nsr-jsr)/tautr;
}
 
void conc_jsr ()
{
	kleak = iupbar/nsrbar;
	ileak = kleak*nsr;
 
	iup = iupbar*cai/(cai+kmup);
 
	dcaiontnew = (caiont-caiontold)/dt;
 
	if(v>-35 && dcaiontnew>dcaiont && flag==0)
		{flag = 1;
		tcicr = 0;}
	
	on = 1/(1+exp((-tcicr+4)/tauon));
	off = (1-1/(1+exp((-tcicr+4)/tauoff)));
	magrel = 1/(1+exp(((ilca+icab+ipca-2*inaca+icat)+5)/0.9));
	
	irelcicr = gmaxrel*on*off*magrel*(jsr-cai);
 
	tcicr = tcicr+dt;
 
	greljsrol = grelbarjsrol*(1-exp(-tjsrol/tauon))*exp(-tjsrol/tauoff);
	ireljsrol = greljsrol*(jsr-cai);
	tjsrol = tjsrol+dt;
 
	csqn = csqnbar*(jsr/(jsr+kmcsqn));
	djsr = dt*(itr-irelcicr-ireljsrol);
	bjsr = csqnbar-csqn-djsr-jsr+kmcsqn;
	cjsr = kmcsqn*(csqn+djsr+jsr);
 
	jsr = (sqrt(bjsr*bjsr+4*cjsr)-bjsr)/2;
}
 
void conc_nsr ()
{
	dnsr = dt*(iup-ileak-itr*vjsr/vnsr);
	nsr = nsr+dnsr;
}
 
 
void conc_cai ()
{
	dcai = -dt*(((caiont*acap)/(vmyo*zca*frdy))+((iup-ileak)*vnsr/vmyo)-(irelcicr*vjsr/vmyo)-(ireljsrol*vjsr/vmyo));
	trpn = trpnbar*(cai/(cai+kmtrpn));
	cmdn = cmdnbar*(cai/(cai+kmcmdn));
 
	catotal = trpn+cmdn+dcai+cai;
	bmyo = cmdnbar+trpnbar-catotal+kmtrpn+kmcmdn;
	cmyo = (kmcmdn*kmtrpn)-(catotal*(kmtrpn+kmcmdn))+(trpnbar*kmcmdn)+(cmdnbar*kmtrpn);
	dmyo = -kmtrpn*kmcmdn*catotal;
	
	gpig = sqrt(bmyo*bmyo-3*cmyo);
 
	cai = (2*gpig/3)*cos(acos((9*bmyo*cmyo-2*bmyo*bmyo*bmyo-27*dmyo)/(2*pow((bmyo*bmyo-3*cmyo),1.5)))/3)-(bmyo/3);
}
 
void conc_cleft()
{
	dnao = dt*((nabm-nao)/taudiff+naiont*acap/(vcleft*frdy));
	nao = dnao+nao;
	dko = dt*((kbm-ko)/taudiff+kiont*acap/(vcleft*frdy));
	ko = dko+ko;
	dcao = dt*((cabm-cao)/taudiff+caiont*acap/(vcleft*frdy*2));
	cao = dcao+cao;
}
 
/* Values are printed to a file called ap. The voltage and currents can be plotted versus time using graphing software. */
 
void prttofile()
{
	if (i >= 0)
	fprintf(ap,"%.3f\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", t-10, v, nai, ki, cai*1000, jsr, nsr, ina, ikr, iks, iki, ilca, icat, inab, icab, ipca, inaca, inak, ikatp,oo,cc3);
 
}