Cardiac Atrial Cell (Courtemanche et al 1998) (C++)

 Download zip file 
Help downloading and running models
Accession:79461
The mechanisms underlying many important properties of the human atrial action potential (AP) are poorly understood. Using specific formulations of the K+, Na+, and Ca2+ currents based on data recorded from human atrial myocytes, along with representations of pump, exchange, and background currents, we developed a mathematical model of the AP. The model AP resembles APs recorded from human atrial samples and responds to rate changes, L-type Ca2+ current blockade, Na+/Ca2+ exchanger inhibition, and variations in transient outward current amplitude in a fashion similar to experimental recordings. Rate-dependent adaptation of AP duration, an important determinant of susceptibility to atrial fibrillation, was attributable to incomplete L-type Ca2+ current recovery from inactivation and incomplete delayed rectifier current deactivation at rapid rates. Experimental observations of variable AP morphology could be accounted for by changes in transient outward current density, as suggested experimentally. We conclude that this mathematical model of the human atrial AP reproduces a variety of observed AP behaviors and provides insights into the mechanisms of clinically important AP properties.
Reference:
1 . Courtemanche M, Ramirez RJ, Nattel S (1998) Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. Am J Physiol 275:H301-21 [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):
Channel(s): I Na,t; I L high threshold; I K; I Sodium; I Calcium; I Potassium; Na/K pump;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s): Acetylcholine;
Simulation Environment: C or C++ program;
Model Concept(s): Ion Channel Kinetics; Action Potentials; Pathophysiology; Heart disease; Sodium pump;
Implementer(s): Wu, Sheng-Nan [snwu at mail.ncku.edu.tw]; Lai, Hsing-Jung ;
Search NeuronDB for information about:  I Na,t; I L high threshold; I K; I Sodium; I Calcium; I Potassium; Na/K pump; Acetylcholine;
/
HAM_AP
readme.html
Cour_HAM_con_1.cpp
ham_ap.jpg
                            
/* The Courtemanche Model of the Human Atrial Myocyte */

/* This code requires a C++ compiler */
/* Acetylcholine-activated K channel current (GIRK3.1) was incoporated */

/* Am. J. Physiol. 1998:275 H301-H321*/

/*  Modified to C by HJ LAI 2006/08/01 */
 
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define bcl 1000     /* Basic cycle length (ms) */
#define beats 50     /* Number of Beats */
#define S2 1000

/* 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.0008;     /* 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;    /* Capacity */
        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) */
        
        /* 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 boolien; /* Boolien 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+1];           /* Max. Voltage (mV) */
        double dvdtmax[beats+1];        /* Max. dv/dt (mV/ms) */
        double apd[beats+1];            /* Action Potential Duration */
        double toneapd[beats+1];        /* Time of dv/dt Max. */
        double ttwoapd[beats+1];        /* Time of 90% Repolarization */
        double trep[beats+1];           /* Time of Full Repolarization */
        double di[beats+1];             /* Diastolic Interval */
        double rmbp[beats+1];           /* Resting Membrane Potential */
        double nair[beats+1];           /* Intracellular Na At Rest */
        double cair[beats+1];           /* Intracellular Ca At Rest */
        double caimax[beats+1];         /* Peak Intracellular Ca */
        int i;                        /* Stimulus 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 ki;     /* Intracellular K Concentration (mM) */
        double ko;     /* Extracellular K Concentration (mM) */
        double cai;    /* Intracellular Ca Concentration (mM) */
        double cao;    /* Extracellular 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) */

        /* Myoplasmic Na Ion Concentration Changes */
        double naiont;  /* Total Na Ion Flow (mM/ms) */
        double dnai;    /* Change in Intracellular Na Concentration (mM) */

        /* Myoplasmic K Ion Concentration Changes */
        double kiont; /* Total K Ion Flow (mM/ms) */
        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.005;  /* 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) */
        double urel;                   /* Activation gate u of Ca release from jsr*/
		double urelss;                 /* Steady state of activation gate u*/
		double tauurel;                /* Time constant of activation gate u*/
		double vrel;                   /* Activation gate v of Ca release from jsr*/
		double vrelss;                 /* Steady state of activation gate v*/
		double tauvrel;                /* Time constant of activation gate v*/
		double wrel;                   /* Inactivation gate w of Ca release from jsr*/
		double wrelss;                 /* Steady state of inactivation gate w*/
		double tauwrel;                /* Time constant of inactivation gate w*/
		double fn;
		const double grelbarjsrol = 30; /* 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 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*/

        /* 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 (mM/ms) */
        double dcai;    /* Change in myoplasmic Ca concentration (mM) */
        double b1cai;
		double b2cai;
		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 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 am;     /* Na alpha-m rate constant (ms^-1) */
        double bm;     /* Na beta-m rate constant (ms^-1) */
        double h;      /* Na activation */
        double j;      /* Na inactivation */
        double m;      /* Na inactivation */
        double gB;
	
        /* Current through L-type Ca Channel */
        double ilca;    /* Ca 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 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 */
        double taufca;  /* Time constant of gate fca (ms^-1) */
		double fcass;   /* Steady-state value of activation gate fca  */

		const double gcalbar = 0.1238;
        
    	/* Acetylcholine-Activated Potassium Current */
		/* modified from Matsuoka et al., Jap J Physiol 2003;53:105-123 */
		double ikach; /* Acetylcholine-activated K current (uA/uF) */
		double gkach; /* Channel conductance of acetylcholine-activated K current (mS/uF) */
        double ekach; /* Reversal potential of acetylcholine-activated K current (mV) */
    	double alphayach; /* Alpha rate constant (ms^-1) */
	    double betayach; /* Beta rate constant (ms^-1) */
	    double tauyach; /* Time constant (ms) */
	    double yachss; /* Steady-state value */
	    double yach;
	    const double ach = 0.0; /* Acetylcholine concentration */

        /* Ultra-Rapidly Activating Potassium Current */
        double ikur;   /* Ultra-rapidly activating K current (uA/uF) */
        double gkur;   /* Channel conductance of ultra-rapidly activating K current (mS/uF) */
        double ekur;   /* Reversal potential of ultra-rapidly activating K current (mV) */
        double uakur;    /* Ultra-rapidly activating K activation gate ua */
        double uakurss;  /* Steady-state value of activation gate ua */
        double tauuakur; /* Time constant of gate ua (ms^-1) */
        double alphauakur; /* Alpha rate constant of activation gate ua (ms^-1) */
		double betauakur;  /* Beta rate constant of activation gate ua (ms^-1) */
		double uikur;    /* Ultra-rapidly activating K activation gate ui*/
        double uikurss;  /* Steady-state value of activation gate ui */
        double tauuikur; /* Time constant of gate ui (ms) */
		double alphauikur; /* Alpha rate constant of activation gate ui (ms^-1) */
		double betauikur; /* Beta rate constant of activation gate ui (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 xs;    /* Slowly activating potassium current activation gate*/
        double xsss;  /* Steady-state value of activation gate xs */
        double tauxs; /* Time constant of gate xs (ms^-1) */
		const double prnak = 0.01833;  /* Na/K Permiability Ratio */
        
        /* Time-Independent Potassium Current */
		/*Partly modified from Matsuoka, et al, Jap J Physiol,2003:53:105-123*/
        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 kin;    /* K inactivation */
		double iku ;	/*Attaching rate constant of Magnesium to iki*/
		double ikl ;    /*Detaching rate constant of Magnesium to iki*/
		double ikay ;   /*Attaching rate constant of spermine to iki*/
		double ikby ;   /*Detaching rate constant of spermine to iki*/
		double tauiky ; /*Time constant of spermine attachment*/
		double ikyss ;  /*Steady state of spermine attachment*/
		double iky ;    /*Spermine attachment*/
		double foiki ;  /*Fraction of channel free from attachment of Magnesium*/
		double fbiki ;  /*Fraction of channel with attachment of Magnesium*/
        
        /* Transient Outward Potassium Current */
        double ito;       /* Transient outward current */
        double gito;      /* Maximum conductance of Ito */
        double erevto;    /* Reversal potential of Ito */
        double ato;       /* Ito activation */
        double alphaato;  /* Ito alpha-a rate constant */
        double betaato;   /* Ito beta-a rate constant */
        double tauato;    /* Time constant of a gate */
        double atoss;     /* Steady-state value of a gate */
        double iito;      /* Ito inactivation */
        double alphaiito; /* Ito alpha-i rate constant */
        double betaiito;  /* Ito beta-i rate constant */
        double tauiito;   /* Time constant of i gate */
        double iitoss;    /* Steady-state value of i gate */
                
        /* Sodium-Calcium Exchanger */
        double inaca;               /* NaCa exchanger current (uA/uF) */
		const double kmnancx = 87.5;  /* Na saturation constant for NaCa exchanger */
		const double ksatncx = 0.1;   /* Saturation factor for NaCa exchanger */
		const double kmcancx = 1.38;  /* Ca saturation factor for NaCa exchanger */
        const double gammas = 0.35;  /* 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 = 1.0933;   /* 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) */
        
        /* Sarcolemmal Ca Pump */
        double ipca;                 /* Sarcolemmal Ca pump current (uA/uF) */
        const double ibarpca = 0.275; /* 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) */
      
        /* Total Ca current */
        double icatot;

        /* Ion Current Functions */     
        void comp_ina ();    /* Calculates Fast Na Current */
        void comp_ical ();   /* Calculates Currents through L-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_ikach();   /* Calculates Acetylcholine-sensitive potassium*/
		void comp_ikur ();   /* Calculates Ultra-Rapidly activation 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_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 main ()
{       
        /* Opening of Datafiles */

        ap = fopen("ap", "w");
/*      
        fpara = fopen("fpara", "w");  
        fmaxs = fopen("fmaxs", "w");
*/

        /* Cell Geometry */
        vcell = 1000*pi*a*a*l;     /*   3.801e-5 uL */
        ageo = 2*pi*a*a+2*pi*a*l;
		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;
        
        /* Time Loop Conditions */
        t = 0;           /* Time (ms) */
        udt = 0.01;     /* Time step (ms) */
        steps = (S2 + bcl*beats)/udt; /* Number of ms */
        st = -200;        /* Stimulus */
        tstim = 10;       /* Time to begin stimulus */
        stimtime = 10;   /* Initial Condition for Stimulus */
        v = -81.2;       /* Initial Voltage (mv) */

        /* Beginning Ion Concentrations */
        nai = 11.2;       /* Initial Intracellular Na (mM) */
        nao = 140;      /* Initial Extracellular Na (mM) */
        ki = 139;       /* Initial Intracellular K (mM) */
        ko = 4.5;       /* Initial Extracellular K (mM) */
        cai = 0.000102;  /* Initial Intracellular Ca (mM) */
        cao = 1.8;      /* Initial Extracellular Ca (mM) */

        /* Initial Gate Conditions */
        m = 0.00291;
        h = 0.965;
        j = 0.978;
        d = 0.000137;
        f = 0.999837;
		xs = 0.0187;
        xr = 0.0000329;
		ato = 0.0304;
		iito = 0.999;
		uakur = 0.00496;
		uikur = 0.999;
		fca = 0.775;
		ireljsrol=0;
        
        /* Initial Conditions */
        jsr = 1.49;
        nsr = 1.49;
		trpn = 0.0118;
        cmdn = 0.00205;
        csqn = 6.51;
        boolien = 1;
        dt = udt;
        utsc = 50;
		urel = 0.00;
		vrel = 1.00;
		wrel = 0.999;
		yach = 2.54e-2;
        iky = 0.6;
		i=-1;

        /* Beginning of Time Loop */
        for (increment = 0; increment < steps; increment++)
        {               
               /* List of functions called for each timestep, currents commented out are only used when modeling pathological conditions */
                comp_ina ();
                comp_ical ();
                comp_ikr ();
                comp_iks ();
                comp_iki ();
				comp_ikach ();
                comp_ikur ();
				comp_ito (); 
                comp_inaca ();
                comp_inak ();
                comp_ipca ();
                comp_icab ();
                comp_inab ();
                comp_it ();
                
                conc_nai ();
                conc_ki ();
                conc_nsr ();
                conc_jsr ();
                calc_itr ();
                conc_cai ();
                
                utsc = 0;
                dt = udt;
		
	prttofile();      

        vnew = v-it*udt;
        dvdtnew = (vnew-v)/udt;

                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 (vnew>=(vmax[i]-0.98*(vmax[i]-rmbp[i])))
                        trep[i] = t;


        v = vnew;
        utsc = utsc+1;      
        t = t+udt;     
        }

        printf ("Main loop passed...\n"); 

        if (beats > 0) {       
          apd[beats] = ttwoapd[beats]-toneapd[beats];
          di[i] = toneapd[beats]-ttwoapd[beats-1]; 
          printf("DI = %g, APD = %g\n", di[beats], apd[beats]);
        }
/*
        fprintf(fpara,"%.3f\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n%g\n", t, v, nai, ki, cai, jsr, nsr, nao, ko, cao, m, h, j, d, f, xs1, xs2, xr, b, g);

        for(i=1;i<beats;i++)
        {apd[i] = ttwoapd[i]-toneapd[i];
         di[i] = toneapd[i]-ttwoapd[i-1]; 
         fprintf(fmaxs, "%i\t%g\t%g\t%g\t%g\n", i, di[i], apd[i], toneapd[i], trep[i]);} 
*/
}

/********************************************************/

/* Functions that describe the currents begin here */


void comp_ina ()
{
        gna = 7.8;
        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.0;
        bh = 1/(0.13*(1+exp((v+10.66)/-11.1)));
        aj = 0.0;
        bj = (0.3*exp(-0.0000002535*v))/(1+exp(-0.1*(v+32)));}

        h = ah/(ah+bh)-((ah/(ah+bh))-h)*exp(-dt/(1/(ah+bh)));
        j = aj/(aj+bj)-((aj/(aj+bj))-j)*exp(-dt/(1/(aj+bj)));
        m = am/(am+bm)-((am/(am+bm))-m)*exp(-dt/(1/(am+bm)));
        
        ina = gna*m*m*m*h*j*(v-ena);
}

void comp_ical ()
{
        dss = 1/(1+exp(-(v+10)/8));
        taud = (1-exp((v+10)/-6.24))/(0.035*(v+10)*(1+exp((v+10)/-6.24)));

        fss = 1/(1+exp((v+28)/6.9));
        tauf = 9/(0.0197*exp(-pow((0.0337*(v+10)),2))+0.02);

		fcass = 1/(1+cai/0.00035);
		taufca = 2;

        d = dss-(dss-d)*exp(-dt/taud);
        f = fss-(fss-f)*exp(-dt/tauf);
		fca = fcass-(fcass-fca)*exp(-dt/tauf);
        
		ibarca = gcalbar*(v-65);                
        
        ilca = d*f*fca*ibarca;
        
        ilcatot = ilca;
}

void comp_ikr ()
{
        gkr = 0.0294*sqrt(ko/5.4);
        ekr = ((R*temp)/frdy)*log(ko/ki);

        xrss = 1/(1+exp(-(v+14.1)/6.5));
        tauxr = 1/(0.0003*(v+14.1)/(1-exp(-(v+14.1)/5))+0.000073898*(v-3.3328)/(exp((v-3.3328)/5.1237)-1));
                
        xr = xrss-(xrss-xr)*exp(-dt/tauxr);
        
        r = 1/(1+exp((v+15)/22.4));
                
        ikr = gkr*xr*r*(v-ekr);
}

void comp_iks ()
{
        gks = 0.129;
        eks = ((R*temp)/frdy)*log(ko/ki);
		tauxs = 0.5/(0.00004*(v-19.9)/(1-exp(-(v-19.9)/17))+0.000035*(v-19.9)/(exp((v-19.9)/9)-1));
		xsss = 1/pow((1+exp(-(v-19.9)/12.7)),0.5);
		xs = xsss-(xsss-xs)*exp(-dt/tauxs);

		iks = gks*xs*xs*(v-eks);
}

void comp_iki ()
{
        gki = 0.09*pow(ko/5.4,0.4);
        eki = ((R*temp)/frdy)*log(ko/ki);

        kin = 1/(1+exp(0.07*(v+80)));
        
        iki = gki*kin*(v-eki);
/*

  // modified from Matsuoka, et al Jap J Physiol 2003:53:105-123
	iku = 0.75*exp(0.035*(v-eki-10))/(1+exp(0.015*(v-eki-140)));
	ikl = 3*exp(-0.048*(v-eki-10))*(1+exp(0.064*(v-eki-38)))/(1+exp(0.03*(v-eki-70)));
	ikay =1/(8000*exp((v-eki-97)/8.5)+7*exp((v-eki-97)/300));
	ikby =1/(0.00014*exp(-(v-eki-97)/9.1)+0.2*exp(-(v-eki-97)/500));
	tauiky = 1/(ikay+ikby);
	ikyss = ikay/(ikay+ikby);
	iky = ikyss - (ikyss-iky)*exp(-dt/tauiky);
	foiki = ikl/(iku+ikl);
	fbiki = iku/(iku+ikl);

	
	iki = gki*(pow(foiki,4)+8*pow(foiki,3)*fbiki/3+2*foiki*foiki*fbiki*fbiki)*iky*(v-eki); 
*/
  }

void comp_ikach ()
{
		gkach = 0.135;
		ekach = ((R*temp)/frdy)*log(ko/ki);
		alphayach= 1.232e-2/(1+0.0042/ach)+0.0002475;
		betayach = 0.01*exp(0.0133*(v+40));
		tauyach = 1/(alphayach+betayach);
		yachss = alphayach/(alphayach+betayach);

		yach = yachss-(yachss-yach)*exp(-dt/tauyach);
		ikach = gkach*yach*(v-ekach)/(1+exp((v+20)/20));
}

void comp_ikur ()
{
		gkur = 0.005+0.05/(1+exp(-(v-15)/13));
        ekur = ((R*temp)/frdy)*log(ko/ki);
		alphauakur = 0.65/(exp(-(v+10)/8.5)+exp(-(v-30)/59.0));
		betauakur = 0.65/(2.5+exp((v+82)/17.0));
		tauuakur = 1/(3*(alphauakur+betauakur));
		uakurss = 1/(1+exp(-(v+30.3)/9.6));
		alphauikur = 1/(21+exp(-(v-185)/28));
		betauikur = exp((v-158)/16);
		tauuikur = 1/(3*(alphauikur+betauikur));
		uikurss = 1/(1+exp((v-99.45)/27.48));

		uakur = uakurss-(uakurss-uakur)*exp(-dt/tauuakur);
		uikur = uikurss-(uikurss-uikur)*exp(-dt/tauuikur);

		ikur = gkur*uakur*uakur*uakur*uikur*(v-ekur);
}

void comp_ito ()
{
        gito = 0.1652;
        erevto = ((R*temp)/frdy)*log(ko/ki);
        
        alphaato = 0.65/(exp(-(v+10)/8.5)+exp(-(v-30)/59));
        betaato = 0.65/(2.5+exp((v+82)/17));
        tauato = 1/(3*(alphaato+betaato));
        atoss = 1/(1+exp(-(v+20.47)/17.54));
        ato = atoss-(atoss-ato)*exp(-dt/tauato);

        alphaiito = 1/(18.53+exp((v+113.7)/10.95));
        betaiito = 1/(35.56+exp(-(v+1.26)/7.44));
        tauiito = 1/(3*(alphaiito+betaiito));
        iitoss = 1/(1+exp((v+43.1)/5.3));
        iito = iitoss-(iitoss-iito)*exp(-dt/tauiito);
        
        ito = gito*ato*ato*ato*iito*(v-erevto);
}

void comp_inaca ()
{
	inaca = 1750*(exp(gammas*frdy*v/(R*temp))*nai*nai*nai*cao-exp((gammas-1)*frdy*v/(R*temp))*nao*nao*nao*cai)/((pow(kmnancx,3)+pow(nao,3))*(kmcancx+cao)*(1+ksatncx*exp((gammas-1)*frdy*v/(R*temp))));
}

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)));
		fnak=(v+150)/(v+200);
        inak = ibarnak*fnak*(1/(1+pow((kmnai/nai),1.5)))*(ko/(ko+kmko));
}

void comp_ipca ()
{
        ipca = (ibarpca*cai)/(kmpca+cai);       
}

void comp_icab ()
{
        gcab = 0.00113;
        ecan = ((R*temp)/frdy)*log(cao/cai);
                
        icab = gcab*(v-ecan);
}

void comp_inab ()
{
        gnab = 0.000674;
        enan = ((R*temp)/frdy)*log(nao/nai);
        
        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+3*inak+3*inaca+1.5e-2;
        kiont = ikr+iks+iki-2*inak+ito+ikur+ikach+1.5e-2;
        caiont = ilca+icab+ipca-2*inaca;
        
        if (t>=tstim && t<(tstim+dt))     
        {stimtime = 0;
        i = i+1;
        if (i < beats-1) 
            tstim = tstim+bcl; 
        else tstim = tstim+S2; 
        printf ("Stimulus %d applied, ", i+1);
        printf ("Time = %f\n", t);
        boolien = 0;
        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;}

        stimtime = stimtime+dt;
}

/* Functions that calculate intracellular ion concentrations begins here */

void conc_nai ()
{
        dnai = -dt*naiont*acap/(vmyo*zna*frdy);
        nai = dnai + nai;
}

void conc_ki ()
{
        dki = -dt*kiont*acap/(vmyo*zk*frdy);
        ki = dki + ki;
}

void conc_nsr ()
{
        kleak = iupbar/nsrbar;
        ileak = kleak*nsr;

        iup = iupbar*cai/(cai+kmup);
		csqn = csqnbar*(jsr/(jsr+kmcsqn));
		
		dnsr = dt*(iup-ileak-itr*vjsr/vnsr);
        nsr = dnsr+nsr;
}

void conc_jsr ()
{

//		fn = vjsr*(1e-12)*ireljsrol-(5e-13)*(ilca/2+inaca/5)*acap/frdy; 
		fn = vjsr*(1e-12)*ireljsrol-(1e-12)*caiont*acap/(2*frdy); 

		tauurel = 8.0;
		urelss = 1/(1+exp(-(fn-3.4175e-13)/13.67e-16));
		tauvrel = 1.91+2.09/(1+exp(-(fn-3.4175e-13)/13.67e-16));
		vrelss = 1-1/(1+exp(-(fn-6.835e-14)/13.67e-16));
		tauwrel = 6.0*(1-exp(-(v-7.9)/5))/((1+0.3*exp(-(v-7.9)/5))*(v-7.9));
		wrelss = 1-1/(1+exp(-(v-40)/17));

		urel = urelss-(urelss-urel)*exp(-dt/tauurel);
		vrel = vrelss-(vrelss-vrel)*exp(-dt/tauvrel);
		wrel = wrelss-(wrelss-wrel)*exp(-dt/tauwrel);
		
        greljsrol = grelbarjsrol*urel*urel*vrel*wrel;
        ireljsrol = greljsrol*(jsr-cai); 

        
		djsr = dt*(itr-0.5*ireljsrol)/(1+csqnbar*kmcsqn/pow((jsr+kmcsqn),2)); //LAI

        jsr = djsr+jsr;
}

void calc_itr ()
{
        itr = (nsr-jsr)/tautr;
}

void conc_cai ()
{
        trpn = trpnbar*(cai/(cai+kmtrpn));
        cmdn = cmdnbar*(cai/(cai+kmcmdn));
		
		b1cai = -caiont*acap/(2*frdy*vmyo)+(vnsr*(ileak-iup)+0.5*ireljsrol*vjsr)/vmyo; //LAI
		b2cai = 1+trpnbar*kmtrpn/pow((cai+kmtrpn),2)+cmdn*kmcmdn/pow((cai+kmcmdn),2);
		dcai = dt*b1cai/b2cai;


        cai = dcai+cai;
}


/* Values are printed to a file called ap. The voltage and currents can be plotted versus time using graphing software. */

void prttofile ()
{
if(increment%100 == 0)
   fprintf (ap, "%.3f\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", t, v, ikr, ilca, iks, ito, cai, iki, ina);
}

Loading data, please wait...