/* 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 */ /* Biophy J 2007;In Press */ #include #include #include #include #include #include #include #define bcl 1000 // Basic Cycle Length (ms) #define beats 10 // Number of Beats #define betaad 0 // Set this Value to 1 to Introduce Beta Adrenergic Effects #define csqnmut 0 // Set this Value to 1 to Introdue CASQ2 Mutation Effects /* List of variables and paramaters (this code uses all global variables) */ /* Creation of Data File */ FILE *ap; FILE *fmaxs; FILE *fpara; FILE *statedata; FILE *ryrdata; 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 volrss; // Restricted Subspace (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]; // 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; // 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 /* Myoplasmic 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 buffersmyo; // Ratio of Buffered [Ca] in the Myoplasm 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) /* SR Ion Concentrations */ 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 (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) /* Myoplasmic Ca Ion Concentration Changes */ double caiont; // Total Ca Ion Flow (uA/uF) double dcai; // Change in myoplasmic Ca 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) double kmup; // Half-saturation concentration of iup (mM) double iupbar; // Max. current through iup channel (mM/ms) double nsrbar; // Max. [Ca] in NSR (mM) /* JSR Ca Ion Concentration Changes */ double djsr; // Change in [Ca] in the JSR (mM) double csqnbar; // Max. [Ca] buffered in CSQN (mM) const double kmcsqn = 0.8; // Equalibrium constant of buffering for CSQN (mM) double caiontold; // Old rate of change of Ca entry double buffersjsr; // Ratio of Buffered [Ca] in the Junctional SR /* Subspace Ion Concentrations */ double caiss; // [Ca] in the Subspace (mM) double naiss; // [Na] in the Subspace (mM) double kiss; // [K] in the Subspace (mM) double dcaresspace; // Change of [Ca] in the Subspace (mM) double dnaresspace; // Change of [Na] in the Subspace (mM) double dkresspace; // Change of [K] in the Subspace (mM) double idiffca; // [Ca] Transfer from Subspace to Myoplasm (mM/ms) double idiffna; // [Na] Transfer from Subspace to Myoplasm (mM/ms) double idiffk; // [K] Transfer from Subspace to Myoplasm (mM/ms) double taudiffss; // Time constant of ion transfer from Subspace to Myoplasm (ms) double bsr; // Anionic SR binding sites for Calcium in the Subspace (mM) const double bsrbar = 0.047; // Max. [Ca] buffered in Anionic SR binding sites (mM) const double kmbsr = 0.00087; // Equalibrium constant of buffering for BSR (mM) double bsl; // Anionic Sarcolemmal binding sites for Calcium in the Subspace (mM) const double bslbar = 2.124; // Max. [Ca] buffered in Anionic Sarcolemmal binding sites (mM) const double kmbsl = 0.127; // Equalibrium constant of buffering for BSL (mM) double buffersss; // Ratio of Buffered [Ca] in the Subspace /* Translocation of Ca Ions from NSR to JSR */ double itr; // Translocation current of Ca ions from NSR to JSR (mM/ms) const double tautr = 120; // Time constant of Ca transfer from NSR to JSR (ms) /* 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 m; // Na activation double h; // Na inactivation double j; // Na inactivation /* 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 /* Sodium-Calcium Exchanger V-S */ double inaca; // NaCa exchanger current (uA/uF) double inacass; // NaCa exchanger current in the Subspace (uA/uF) double c1; // NaCa Scaling factor double c2; // Half-saturation concentration of NaCa exhanger (mM) double c3; // NaCa exhanger Extracellular Na Dependence double gammanaca; // 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 double ibarnak; // 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) double ibarpca; // Max. Ca current through sarcolemmal Ca pump (uA/uF) /* 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) /* L-type Calcium Channel Markov Formulation */ void comp_ltypesc (); // Function That Will Compute L-type Current void comp_ltypeca (); // Calculates Currents through Markov L-Type Ca Channel void comp_rates (); // Calculates the Rate Transitions of L-type Markov Model double ical; // Ca current through L-type Ca channel (uA/uF) 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) const double pca = 0.006615; // Permiability of membrane to Ca (cm/s) const double gacai = 0.01; // Activity coefficient of Ca const double gacao = 0.341; // Activity coefficient of Ca const double pna = 0.000008265; // 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.000002363; // 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 double ltypeCzero; double dltypeCzero; double Czero1; double Czero2; double ltypeCone; double dltypeCone; double Cone1; double Cone2; double ltypeCtwo; double dltypeCtwo; double Ctwo1; double Ctwo2; double ltypeCthree; double dltypeCthree; double Cthree1; double Cthree2; double ltypeO; double dltypeO; double O1; double O2; double ltypeIVf; double dltypeIVf; double IVf1; double IVf2; double ltypeIVs; double dltypeIVs; double IVs1; double IVs2; double icaCzero; double icaCone; double icaCtwo; double icaCthree; double icaO; double icaIVf; double icaIVs; double icaCzerocdi; double icaConecdi; double icaCtwocdi; double icaCthreecdi; double icaIcdi; double icaIVfcdi; double icaIVscdi; double alphaequation; double betaequation; double alpha0; double alpha1; double alpha2; double alpha3; double beta0; double beta1; double beta2; double beta3; double gammmaf; double gammmas; double phif; double phis; double lambdaf; double lambdas; double omegaf; double omegas; double omegasf; double omegafs; double fcasc; double caon; double caoff; double fmode0; double mode0on; double mode0off; /* RyR Channal Markov Formulation */ void comp_ryrrates(); // Computes RyR Rates void comp_ryr(); // Computes RyR State Residencies double sponrel; // Spontaneous SR Ca Release Function double sponrelss; // JSR dependence of Spontaneous SR Ca Release double spontau; // Time constant of Spontaneous Release double gradedrel; // Controls Graded Release double vgainofrel; // Controls Variable Dependent Gain of Relase double vg; double grel; // Irel Scalar double irel; // SR Ca Release double ryrCone; double dryrCone; double ryrCone1; double ryrCone11; double ryrCone2; double ryrCtwo; double dryrCtwo; double ryrCtwo1; double ryrCtwo11; double ryrCtwo2; double ryrCthree; double dryrCthree; double ryrCthree1; double ryrCthree11; double ryrCthree2; double ryrCfour; double dryrCfour; double ryrCfour1; double ryrCfour11; double ryrCfour2; double ryrOone; double dryrOone; double ryrOone1; double ryrOone11; double ryrOone2; double ryrIone; double dryrIone; double ryrIone1; double ryrIone11; double ryrIone2; double ryrItwo; double dryrItwo; double ryrItwo1; double ryrItwo11; double ryrItwo2; double ryrIthree; double dryrIthree; double ryrIthree1; double ryrIthree11; double ryrIthree2; double ryrIfour; double dryrIfour; double ryrIfour1; double ryrIfour11; double ryrIfour2; double ryrIfive; double dryrIfive; double ryrIfive1; double ryrIfive11; double ryrIfive2; double C1_C2; double C2_C3; double C3_C4; double C4_O1; double I1_I2; double I2_I3; double I3_I4; double I4_I5; double O1_C4; double C4_C3; double C3_C2; double C2_C1; double I5_I4; double I4_I3; double I3_I2; double I2_I1; double C1_I1; double C2_I2; double C3_I3; double C4_I4; double O1_I5; double I1_C1; double I2_C2; double I3_C3; double I4_C4; double I5_O1; /* 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_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 comp_Cacon (); // Calculates Ca Concentration in Restricted Space int main () { /* Opening of Datafiles */ ap = fopen("ap", "w"); fpara = fopen("fpara", "w"); fmaxs = fopen("fmaxs", "w"); statedata = fopen("statedata", "w"); ryrdata = fopen("ryrdata", "w"); printdata = 60; /* Cell Geometry */ vcell = 1000*pi*a*a*l; // vell = 3.801e-5 uL ageo = 2*pi*a*a+2*pi*a*l; // ageo = 7.671e-5 cm^2 acap = ageo*2; // acap = 1.534e-4 cm^2 volrss = vcell*0.02; // volrss = 7.602e-7 uL vmyo = vcell*0.68; // vmyo = 2.584e-5 uL vmito = vcell*0.24; // vmito = 9.122e-6 uL vsr = vcell*0.06; // vsr = 2.281e-6 uL vnsr = vcell*0.0552; // vnsr = 2.098e-6 uL vjsr = vcell*0.0048; // vjsr = 1.824e-7 uL /* Time Loop Conditions */ t = 0; // Time (ms) udt = 0.001; // Time step (ms) steps = (bcl*beats)/udt; // Number of ms st = -100; // Stimulus tstim = 10; // Time to begin stimulus stimtime = 10; // Initial Condition for Stimulus v = -91.34; // Initial Voltage (mv) /* Beginning Ion Concentrations */ nai = 10.71; // Initial Intracellular Na (mM) naiss = 10.71; // Initial Subspace Na (mM) nao = 130; // Initial Extracellular Na (mM) ki = 137.6; // Initial Intracellular K (mM) kiss = 137.6; // Initial Subspace K (mM) ko = 4.5; // Initial Extracellular K (mM) cai = 0.0000821; // Initial Intracellular Ca (mM) caiss = 0.0000672; // Initial Subspace Ca (mM) cao = 1.8; // Initial Extracellular Ca (mM) /* Initial Gate Conditions */ m = 0.000642; h = 0.996271; j = 0.997286; xs1 = 0.004761; xs2 = 0.03291; xr = 0.0001096; b = 0.000880231; g = 0.994305; fmode0 = .22; fcasc = 1; /* Initial Conditions */ jsr = 2.479; nsr = 2.553; trpn = 0.00987; cmdn = 0.00166; csqn = 7.56; boolien = 1; dt = udt; utsc = 50; i=-1; ryrCone = 0.89; ryrCtwo = 0.021; ryrCthree = 0.003; ryrCfour = 0.001; ryrOone = 0; ryrIone = 0.082; ryrItwo = 0.002; ryrIthree = 0.001; ryrIfour = 0; ryrIfive = 1-(ryrCone+ryrCtwo+ryrCthree+ryrCfour+ryrOone+ryrIone+ryrItwo+ryrIthree+ryrIfour); ltypeCzero = 0.948; ltypeCone = 0.052; ltypeCtwo = 0; ltypeCthree = 0.; ltypeIVf = 0; ltypeIVs = 0; ltypeO = 1-(ltypeCzero+ltypeCone+ltypeCtwo+ltypeCthree+ltypeIVf+ltypeIVs); /* Beginning of Time Loop */ for (increment = 0; increment < steps; increment++) { if(v<-50) {nxstep = 10;} else {nxstep = 1;} if(utsc>=nxstep || dvdt>5 || irel>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_ltypesc (); comp_icat (); comp_ikr (); comp_iks (); comp_iki (); comp_ikp (); comp_inaca (); comp_inak (); comp_ipca (); comp_icab (); comp_inab (); comp_it (); conc_nai (); conc_ki (); conc_nsr (); comp_ryrrates (); comp_ryr (); conc_jsr (); calc_itr (); conc_cai (); stimtime = stimtime+dt; vnew = v-it*dt; dvdtnew = (vnew-v)/dt; utsc = 0; dt = 0; } if(v>-85 || (stimtime>=0 && stimtime<0.5)) {printval = 500;} else {printval = 10000;} if(printdata>=printval) {prttofile(); printdata = 0;} printdata = printdata+1; 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; } v = vnew; dvdt = dvdtnew; caiontold = caiont; dt = dt+udt; utsc = utsc+1; t = t+udt; } 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%g\n", t, v, nai, ki, cai, naiss, kiss, caiss, jsr, nsr, nao, ko, cao, m, h, j, xs1, xs2, xr, b, g); for(i=0;i=tstim && t<(tstim+dt)) {stimtime = 0; tstim = tstim+bcl; boolien = 0; i = i+1; rmbp[i]=v; nair[i] = nai; cair[i] = cai; kir[i] = ki;} if(stimtime>=0 && stimtime<0.5) {it = st+naiont+kiont+caiont; cout << t << endl;cout << v << endl;cout << nair[i] << endl;} 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-3*inacass-ilcana)*acap)/(vmyo*zna*frdy)-idiffna*(volrss/vmyo)); nai = dnai + nai; } void conc_ki () { if(stimtime>=0 && stimtime<0.5) {dki = -dt*(((kiont-ilcak+st)*acap)/(vmyo*zk*frdy)-idiffk*(volrss/vmyo));} else {dki = -dt*(((kiont-ilcak)*acap)/(vmyo*zk*frdy)-idiffk*(volrss/vmyo));} ki = dki + ki; } void conc_nsr () { kmup = 0.00092; iupbar = 0.017325; nsrbar = 15; kleak = 0.00875/nsrbar; ileak = kleak*nsr; if(betaad==0) { iup = iupbar*cai/(cai+kmup); } if(betaad==1) { iup = 1.2*iupbar*cai/(cai+kmup); } dnsr = dt*(iup-ileak-itr*vjsr/vnsr); nsr = nsr+dnsr; } void comp_ryrrates() { C1_C2 = 1750*caiss; C2_C3 = 5600*caiss; C3_C4 = 5600*caiss; C4_O1 = 5600*caiss; I1_I2 = 1750*caiss; I2_I3 = 5600*caiss; I3_I4 = 5600*caiss; I4_I5 = 5600*caiss; C2_C1 = 5; C3_C2 = 2.625; C4_C3 = 1; O1_C4 = 6.25; I2_I1 = 5; I3_I2 = 2.625; I4_I3 = 1; I5_I4 = 6.25; C1_I1 = 0.4*caiss; C2_I2 = 1.2*caiss; C3_I3 = 2.8*caiss; C4_I4 = 5.2*caiss; O1_I5 = 8.4*caiss; I1_C1 = 0.01/(1+pow((csqnbar*0.75/csqn),9)); I2_C2 = 0.001/(1+pow((csqnbar*0.75/csqn),9)); I3_C3 = 0.0001/(1+pow((csqnbar*0.75/csqn),9)); I4_C4 = 0.00001/(1+pow((csqnbar*0.75/csqn),9)); I5_O1 = 0.000001/(1+pow((csqnbar*0.75/csqn),9)); } void comp_ryr() { ryrCone1 = ryrCone; ryrCtwo1 = ryrCtwo; ryrCthree1 = ryrCthree; ryrCfour1 = ryrCfour; ryrOone1 = ryrOone; ryrIone1 = ryrIone; ryrItwo1 = ryrItwo; ryrIthree1 = ryrIthree; ryrIfour1 = ryrIfour; ryrIfive1 = 1-(ryrCone1+ryrCtwo1+ryrCthree1+ryrCfour1+ryrOone1+ryrIone1+ryrItwo1+ryrIthree1+ryrIfour1); dryrCone = dt*(C2_C1*ryrCtwo1+I1_C1*ryrIone1-ryrCone1*(C1_C2+C1_I1)); dryrCtwo = dt*(C1_C2*ryrCone1+C3_C2*ryrCthree1+I2_C2*ryrItwo1-ryrCtwo1*(C2_C1+C2_C3+C2_I2)); dryrCthree = dt*(C2_C3*ryrCtwo1+C4_C3*ryrCfour1+I3_C3*ryrIthree1-ryrCthree1*(C3_C2+C3_C4+C3_I3)); dryrCfour = dt*(O1_C4*ryrOone1+C3_C4*ryrCthree1+I4_C4*ryrIfour1-ryrCfour1*(C4_C3+C4_O1+C4_I4)); dryrOone = dt*(C4_O1*ryrCfour1+I5_O1*ryrIfive1-ryrOone1*(O1_I5+O1_C4)); dryrIone = dt*(C1_I1*ryrCone1+I2_I1*ryrItwo1-ryrIone1*(I1_C1+I1_I2)); dryrItwo = dt*(C2_I2*ryrCtwo1+I1_I2*ryrIone1+I3_I2*ryrIthree1-ryrItwo1*(I2_C2+I2_I1+I2_I3)); dryrIthree = dt*(C3_I3*ryrCthree1+I2_I3*ryrItwo1+I4_I3*ryrIfour1-ryrIthree1*(I3_C3+I3_I2+I3_I4)); dryrIfour = dt*(C4_I4*ryrCfour1+I3_I4*ryrIthree1+I5_I4*ryrIfive1-ryrIfour1*(I4_C4+I4_I3+I4_I5)); dryrIfive = dt*(O1_I5*ryrOone1+I4_I5*ryrIfour1-ryrIfive1*(I5_O1+I5_I4)); ryrCone2 = ryrCone+dryrCone; ryrCtwo2 = ryrCtwo+dryrCtwo; ryrCthree2 = ryrCthree+dryrCthree; ryrCfour2 = ryrCfour+dryrCfour; ryrOone2 = ryrOone+dryrOone; ryrIone2 = ryrIone+dryrIone; ryrItwo2 = ryrItwo+dryrItwo; ryrIthree2 = ryrIthree+dryrIthree; ryrIfour2 = ryrIfour+dryrIfour; ryrIfive2 = ryrIfive+dryrIfive; ryrCone = ryrCone2; ryrCtwo = ryrCtwo2; ryrCthree = ryrCthree2; ryrCfour = ryrCfour2; ryrOone = ryrOone2; ryrIone = ryrIone2; ryrItwo = ryrItwo2; ryrIthree = ryrIthree2; ryrIfour = ryrIfour2; ryrIfive = 1-(ryrCone2+ryrCtwo2+ryrCthree2+ryrCfour2+ryrOone2+ryrIone2+ryrItwo2+ryrIthree2+ryrIfour2); } void conc_jsr () { sponrelss = 25/(1+exp((5.3-jsr)/.001)); spontau = 125; sponrel = sponrelss-(sponrelss-sponrel)*exp(-dt/spontau); gradedrel = (1/(1+exp((20+ical)/6))-0.034445); vgainofrel = (1+exp(((0.015*ibarca)+1.25)/0.75)); vg = gradedrel/vgainofrel; grel = 250*(vg+sponrel); irel = grel*(ryrOone)*(jsr-caiss); if (csqnmut == 0) {csqnbar = 10;} if (csqnmut == 1) {csqnbar = 4;} csqn = csqnbar*(jsr/(jsr+kmcsqn)); buffersjsr = 1/(1+(csqnbar*kmcsqn/(pow(kmcsqn+jsr,2)))); djsr = dt*buffersjsr*(itr-irel); jsr = jsr+djsr; } void calc_itr () { itr = (nsr-jsr)/tautr; } void conc_cai () { trpn = trpnbar*(cai/(cai+kmtrpn)); cmdn = cmdnbar*(cai/(cai+kmcmdn)); buffersmyo = 1/(1+(trpnbar*kmtrpn)/((kmtrpn+cai)*(kmtrpn+cai))+(cmdnbar*kmcmdn)/((kmcmdn+cai)*(kmcmdn+cai))); dcai = -dt*buffersmyo*((((caiont-ilca+2*inacass)*acap)/(vmyo*zca*frdy))+((iup-ileak)*(vnsr/vmyo))-idiffca*(volrss/vmyo)); cai = cai+dcai; } void comp_ltypesc () { comp_rates (); comp_ltypeca (); comp_Cacon (); ibarca = pca*zca*zca*((v*frdy*frdy)/(R*temp))*((gacai*caiss*exp((zca*(v-0)*frdy)/(R*temp))-gacao*cao)/(exp((zca*(v-0)*frdy)/(R*temp))-1)); ibarna = pna*zna*zna*((v*frdy*frdy)/(R*temp))*((ganai*naiss*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*kiss*exp((zk*v*frdy)/(R*temp))-gako*ko)/(exp((zk*v*frdy)/(R*temp))-1)); fcasc = caoff/(caon+caoff)-((caoff/(caon+caoff))-fcasc)*exp(-dt/(1/(caon+caoff))); fmode0 = mode0off/(mode0on+mode0off)-((mode0off/(mode0on+mode0off))-fmode0)*exp(-dt/(1/(mode0on+mode0off))); ical = ibarca*fcasc*fmode0*ltypeO; ilcana = ibarna*fcasc*fmode0*ltypeO; ilcak = ibark*fcasc*fmode0*ltypeO; ilca = ical; } void comp_rates () { if(betaad==0) { alphaequation = 0.925*exp(v/30); betaequation = 0.39*exp(v/-40); alpha0 = 4*alphaequation; alpha1 = 3*alphaequation; alpha2 = 2*alphaequation; alpha3 = 1*alphaequation; beta0 = 1*betaequation; beta1 = 2*betaequation; beta2 = 3*betaequation; beta3 = 4*betaequation; gammmaf = 0.245*exp(v/10); gammmas = 0.005*exp(v/-40); phif = 0.02*exp(v/500); phis = 0.03*exp(v/-280); lambdaf = 0.035*exp(v/-300); lambdas = 0.0011*exp(v/500); omegaf = (lambdaf*beta3*gammmaf)/(phif*alpha3); omegas = (lambdas*beta3*gammmas)/(phis*alpha3); omegasf = (lambdas*phif)/(lambdaf); omegafs = phis; caon = 4/(1+1/caiss); caoff = 0.01; mode0on = 0.0008; mode0off = 0.00018; } if(betaad==1) { alphaequation = 0.925*exp(v/30); betaequation = 0.39*exp(v/-40); alpha0 = 4*alphaequation; alpha1 = 3*alphaequation; alpha2 = 2*alphaequation; alpha3 = 1*alphaequation; beta0 = 1*betaequation; beta1 = 2*betaequation; beta2 = 3*betaequation; beta3 = 4*betaequation; gammmaf = 0.245*exp(v/10); gammmas = 0.005*exp(v/-40); phif = 0.02*exp(v/500); phis = 0.014*exp(v/-20); lambdaf = 0.015*exp(v/-300); lambdas = 0.0022*exp(v/500); omegaf = (lambdaf*beta3*gammmaf)/(phif*alpha3); omegas = (lambdas*beta3*gammmas)/(phis*alpha3); omegasf = (lambdas*phif)/(lambdaf); omegafs = phis; caon = 5/(1+1/caiss); caoff = 0.001; mode0on = 0.0004; mode0off = 0.0004545; } } void comp_ltypeca() { Czero1 = ltypeCzero; Cone1 = ltypeCone; Ctwo1 = ltypeCtwo; Cthree1 = ltypeCthree; IVf1 = ltypeIVf; IVs1 = ltypeIVs; O1 = 1-(Czero1+Cone1+Ctwo1+Cthree1+IVf1+IVs1); dltypeCzero = dt*(Cone1*beta0-Czero1*(alpha0)); dltypeCone = dt*(Ctwo1*beta1+Czero1*alpha0-Cone1*(alpha1+beta0)); dltypeCtwo = dt*(Cone1*alpha1+Cthree1*beta2-Ctwo1*(beta1+alpha2)); dltypeCthree = dt*(O1*beta3+IVf1*omegaf+IVs1*omegas+Ctwo1*alpha2-Cthree1*(gammmaf+gammmas+alpha3+beta2)); dltypeO = dt*(IVf1*lambdaf+IVs1*lambdas+Cthree1*alpha3-O1*(phif+phis+beta3)); dltypeIVf = dt*(O1*phif+Cthree1*gammmaf+IVs1*omegasf-IVf1*(omegaf+lambdaf+omegafs)); dltypeIVs = dt*(O1*phis+Cthree1*gammmas+IVf1*omegafs-IVs1*(omegas+lambdas+omegasf)); Czero2 = ltypeCzero+dltypeCzero; Cone2 = ltypeCone+dltypeCone; Ctwo2 = ltypeCtwo+dltypeCtwo; Cthree2 = ltypeCthree+dltypeCthree; O2 = ltypeO+dltypeO; IVf2 = ltypeIVf+dltypeIVf; IVs2 = ltypeIVs+dltypeIVs; ltypeCzero = Czero2; ltypeCone = Cone2; ltypeCtwo = Ctwo2; ltypeCthree = Cthree2; ltypeIVf = IVf2; ltypeIVs = IVs2; ltypeO = 1-(Czero2+Cone2+Ctwo2+Cthree2+IVf2+IVs2); } void comp_Cacon () { taudiffss = 0.1; bsr = bsrbar*(caiss/(caiss+kmbsr)); bsl = bslbar*(caiss/(caiss+kmbsl)); buffersss = 1/(1+(bsrbar*kmbsr)/((kmbsr+caiss)*(kmbsr+caiss))+(bslbar*kmbsl)/((kmbsl+caiss)*(kmbsl+caiss))); idiffca = (caiss-cai)/taudiffss; dcaresspace = dt*buffersss*(-(ical-2*inacass)*acap/(volrss*frdy*2)+(irel*(vjsr/volrss))-idiffca); caiss = caiss+dcaresspace; idiffna = (naiss-nai)/taudiffss; dnaresspace = dt*(-(3*inacass+ilcana)*acap/(volrss*frdy)-idiffna); naiss = naiss+dnaresspace; idiffk = (kiss-ki)/taudiffss; dkresspace = dt*(-ilcak*acap/(volrss*frdy)-idiffk); kiss = kiss+dkresspace; } /* Values are printed to a file called ap. The voltage and currents can be plotted versus time using graphing software. */ void prttofile() { /* These equations compute the state residencies for the 14 state model; the first 7 equations are residencies in ModeV and the latter 7 equations are residencies in ModeCa */ icaCzero = ltypeCzero*fcasc; icaCone = ltypeCone*fcasc; icaCtwo = ltypeCtwo*fcasc; icaCthree = ltypeCthree*fcasc; icaO = ltypeO*fcasc; icaIVf = ltypeIVf*fcasc; icaIVs = ltypeIVs*fcasc; icaCzerocdi = ltypeCzero*(1-fcasc); icaConecdi = ltypeCone*(1-fcasc); icaCtwocdi = ltypeCtwo*(1-fcasc); icaCthreecdi = ltypeCthree*(1-fcasc); icaIcdi = ltypeO*(1-fcasc); icaIVfcdi = ltypeIVf*(1-fcasc); icaIVscdi = ltypeIVs*(1-fcasc); /* Prints to File */ if(t>(0) && t<(300000)) {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\n", t, v, 1000*cai, ilca, irel, inaca, inacass, inak, ikr, iks, iki, ikp, ina, inab, icab, icat, nsr, jsr, csqn);} if(t>(0) && t<(300000)) {fprintf(statedata,"%.3f\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", t, ilca, icaO, icaIcdi, icaIVf, icaIVfcdi, icaIVs, icaIVscdi, (icaCzero+icaCone+icaCtwo+icaCthree), (icaCzerocdi+icaConecdi+icaCtwocdi+icaCthreecdi));} if(t>(0) && t<(300000)) {fprintf(ryrdata,"%.3f\t%g\t%g\t%g\t%g\t%g\n", t, irel, (ryrCone+ryrCtwo+ryrCthree+ryrCfour), ryrOone, (ryrIone+ryrItwo+ryrIthree), (ryrIfour+ryrIfive));} }