Simulation study of Andersen-Tawil syndrome (Sung et al 2006)

 Download zip file 
Help downloading and running models
Accession:79237
Patients with Andersen-Tawil syndrome (ATS) mostly have mutations on the KCNJ2 gene producing loss of function or dominant-negative suppression of the inward rectifier K(+) channel Kir2.1. However, clinical manifestations of ATS including dysmorphic features, periodic paralysis (hypo-, hyper-, or normokalemic), long QT, and ventricular arrhythmias (VA) are considerably variable. Using a modified dynamic Luo-Rudy simulation model of cardiac ventricular myocyte, we elucidate the mechanisms of VA in ATS. We adopted a kinetic model of KCNJ2 in which channel block by Mg(+2) and spermine was incorporated. In this study, we attempt to examine the effects of KCNJ2 mutations on the ventricular action potential (AP), single-channel Markovian models were reformulated and incorporated into the dynamic Luo-Rudy model for rapidly and slowly delayed rectifying K(+) currents and KCNJ2 channel. During pacing at 1.0 Hz with [K(+)]o at 5.4 mM, a stepwise 10% reduction of Kir2.1 channel conductance progressively prolonged the terminal repolarization phase of AP along with gradual depolarization of the resting membrane potential (RMP). At 90% reduction, early after- depolarizations (EADs) became inducible and RMP was depolarized to -55.0 mV (control: -90.1 mV) followed by emergence of spontaneous action potentials (SAP). Both EADs and SAP were facilitated by a decrease in [K(+)]o and suppressed by increase in [K(+)]o. beta-adrenergic stimulation enhanced delayed after-depolarizations (DADs) and could also facilitate EADs as well as SAP in the setting of low [K(+)]o and reduced Kir2.1 channel conductance. In conclusion, the spectrum of VA in ATS includes (1) triggered activity mediated by EADs and/or DADs, and (2) abnormal automaticity manifested as SAP. These VA can be aggravated by a decrease in [K(+)]o and beta-adrenergic stimulation, and may potentially induce torsades de pointes and cause sudden death. In patients with ATS, the hypokalemic form of periodic paralysis should have the highest propensity to VA especially during physical activities.
Reference:
1 . Sung RJ, Wu SN, Wu JS, Chang HD, Luo CH (2006) Electrophysiological mechanisms of ventricular arrhythmias in relation to Andersen-Tawil syndrome under conditions of reduced IK1: a simulation study. Am J Physiol Heart Circ Physiol 291:H2597-605 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell; Electrogenic pump;
Brain Region(s)/Organism:
Cell Type(s): Heart cell;
Channel(s): I Na,t; I L high threshold; I T low threshold; I K; I Mixed; I Potassium; Na/K pump;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s): Ions;
Simulation Environment: C or C++ program;
Model Concept(s): Activity Patterns; Ion Channel Kinetics; Action Potentials; Heart disease; Long-QT; Sodium pump;
Implementer(s): Wu, Sheng-Nan [snwu at mail.ncku.edu.tw]; Chang, Han-Dong;
Search NeuronDB for information about:  I Na,t; I L high threshold; I T low threshold; I K; I Mixed; I Potassium; Na/K pump; Ions;
/
ATS-file
readme.html
ATS-AP.cpp
Fig7.jpg
                            
/* 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);
 
} 

 


Loading data, please wait...