Research

LRd MODEL C++ CODE

/* The Luo-Rudy Dynamic (LRd) Model of the Mammalian Ventricular Myocyte */
/* Gregory Faber */
/* This code requires a C++ compiler */
/* Detailed list of equations and model description are provided in */

/* Circ Res 1991;68:1501-1526 */
/* Circ Res 1994;74:1071-1096 */
/* Circ Res 1994;74:1097-1113 */
/* Circ Res 1995;77:140-152 */
/* Biophys J 1995;68:949-964 */
/* Cardiovasc Res 1997;35:256-272 */
/* Circulation 1999;99:2466-2474 */
/* Cardiovas Res 1999;42:530-542 */
/* Nature 1999;400:566-569 */
/* Circulation 2000;101:1192-1198 */
/* Biophy J 2000;78:2392-2404 */


/* 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 30 // Number of Beats

/* List of variables and parameters (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 kir [beats] ; // Intracellular K 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; // 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; // Equilibrium 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; // Equilibrium constant of buffering for CMDN (mM)
const double kmtrpn = 0.0005; // Equilibrium 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 pca = 0.00054; // Permeability 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; // Permeability 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; // Permeability 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-independent 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 Permeability Ratio

/* Potassium Current (time-independent) */
double iki; // Time-independent 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 availability of open channels
const double natp = 0.24; // K dependence of ATP-sensitive K current
const double nicholsarea = 0.00005; // Nichol's area (cm^2)
const double atpi = 3; // Intracellular ATP concentraion (mM)
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 independent 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 dependence of inaca

/* Sodium-Potassium Pump */
double inak; // NaK pump current (uA/uF)
double fnak; // Voltage-dependence parameter of inak
double sigma; // [Na] o dependence 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; // Permeability 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-Independent 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.01; // 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 = 12.236437; // 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; // Initial Extracellular K (mM)
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;

/* 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 = 50;}
else
{nxstep = 500;}

if(utsc>=nxstep || dvdt>5 || irelcicr>0.01 || (t>=(tstim-udt) && t<=(tstim+udt))
|| (stimtime>=0 && stimtime<=0.5))
{

/* List of functions called for each timestep, currents commented
out are only used when modeling pathological conditions */

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 ();

//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 */

stimtime = stimtime+dt;

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

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 occurred at time " << t << endl;
}

utsc = 0;
dt = 0;

}

if(dvdt>3 || irelcicr>.01)
{printval = 50;}
else
{printval = 750;}
if(printdata>=printval)
{prttofile();
printdata = 0;}
printdata = printdata+1;


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\t%g\n", i, vmax [i] , dvdtmax [i] , apd [i] ,
toneapd [i] , ttwoapd [i] , nair [i] , kir [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);
}

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

xrss = 1/(1+exp(-(v+21.5)/7.5));
tauxr = 1/(0.00138*(v+14.2)/(1-exp(-0.123*(v+14.2)))+0.00061*(v+38.9)/(exp(0.145*(v+38.9))-1));
xr = xrss-(xrss-xr)*exp(-dt/tauxr);
r = 1/(1+exp((v+9)/22.4));
ikr = gkr*xr*r*(v-ekr);
}

void comp_iks ()
{
gks = 0.433*(1+0.6/(1+pow((0.000038/cai),1.4)));
eks = ((R*temp)/frdy)*log((ko+prnak*nao)/(ki+prnak*nai));

xs1ss = 1/(1+exp(-(v-1.5)/16.7));
xs2ss = xs1ss;
tauxs1 = 1/(0.0000719*(v+30)/(1-exp(-0.148*(v+30)))+0.000131*(v+30)/(exp(0.0687*(v+30))-1));
tauxs2 = 4*tauxs1;
xs1 = xs1ss-(xs1ss-xs1)*exp(-dt/tauxs1);
xs2 = xs2ss-(xs2ss-xs2)*exp(-dt/tauxs2);
iks = gks*xs1*xs2*(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 = 0.000195/nicholsarea;
patp = 1/(1+(pow((atpi/katp),hatp)));
gkbaratp = gkatp*patp*(pow((ko/4),natp));

ikatp = gkbaratp*(v-ekatp);
}

void comp_ito ()
{
gitodv = 0.5;
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;
kir [i] = ki;
cair [i] = cai;
cout << t << endl;}

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;

if(stimtime>=0 && stimtime<=0.5)
{dko = dt*((kbm-ko)/taudiff+(kiont+st)*acap/(vcleft*frdy));}
else
{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(t>(0) && t<(bcl*beats))
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\n",
t, v, nai, ki, cai, jsr, nsr, ina, ikr, iks, iki, ilca, icat, inab, icab, ipca, inaca, inak);

}

 

Return to LRd Model

Copyright © 2005-2008 Yoram Rudy, Ph.D. All Rights Reserved.  •  Contact Webmaster