Location: Aslanidi_etal_2010 @ f9e156454803 / M2010.cpp

Author:
Dougal Cowan <devnull@localhost>
Date:
2012-06-19 14:50:18+12:00
Desc:
edited model links
Permanent Source URI:
https://models.physiomeproject.org/w/dcowan/Aslanidi_etal_2010/rawfile/f9e1564548039eb6ee890cb9463d897ae8a732da/M2010.cpp

/*************************************************************************
Rabbit Midmyocardial cell model, developed by Rakan Sleiman & Oleg Aslanidi
Results are published in Aslanidi et al, Biophys J. 2010; 98(11): 2420-31
Created: 25-02-2009, Updated: 14-06-2012. Email: oleg.aslanidi@kcl.ac.uk
*************************************************************************/

#include <cmath>
#include <cstdio>
#include <iostream>
#include <fstream>

using namespace std;

/***************Simulation constants**********************/
//#define bcl	350
double bcl;
#define beats	20
#define S2	bcl
#define HT	0.005         /* ms */

/*****************Basic physical constants****************/

#define RTF 26.7081
const float Faraday=96487;
#define zCa 2.0
#define zNa 1.0
#define zK 1.0

#define PNaK 0.01833    /* IKs - for determining EKs*/
#define fCaKd 7.0  	/* uM */
/***************Declarations**************************/

double INa,INal;
double IKs, IKr, IK1, IKp,IKp_Kp,ITO_slow, ITO_fast, ITO;
double ICa, ICl;
double IClB,ICaB, INaB;
double INaK, ISLCaP, INaCaX;
double Iion;
double Istim;

double EK, ECl,ENa, ECa, EKs;

double V, m, h, j, x, X2, d, f,xr,a, i, xto_slow, yto_slow, r_slow, xto_fast, yto_fast,fCa,f2,ml,hl;
double NaIono ,NaIoni,CaIono ,KIono ,KIoni,ClIono ,ClIoni;

double m_inf, h_inf, j_inf, alpha_m, beta_m, alpha_h, beta_h, alpha_j, beta_j, tau_m, tau_h, tau_j;
double alpha_ml,beta_ml,ml_inf,tau_ml,tau_hl,hl_inf;
double d_inf, f_inf, tau_d, tau_f,tau_f1,tau_f2;
double x_inf, tau_x;
double pCa;                                                                          /* IKs */
double xr_inf, tau_xr, GKr, r_inf, tau_r;
double alpha_k1, beta_k1, k1_inf;
double xto_slow_inf, yto_slow_inf, r_slow_inf, tau_xto_slow, tau_yto_slow, tau_r_slow;
double xto_fast_inf, yto_fast_inf, tau_xto_fast, tau_yto_fast;
double a_inf, i_inf, tau_i;

double GNa;
double GNal;
double GCaL;
double SLGKs;
double GKs;
double GKr_factor;
double GK1;
double GTO_slow;
double GTO_fast;
double Gtos_factor;
double Gtof_factor;
double Gto_factor;
double GCl;
double GClB;
double GCaB;
double GNaB;


double INaK_Kmr;
double INaK_Hill;
double sigma, fNaK, INaK_Vmf, INaK_Kmf;

double INaCaXCa_Kmf;        /* mM */
double INaCaXCa_Kmr;            /* mM */
double INaCaXNa_Kmf;          /* mM */
double INaCaXNa_Kmr;            /* mM */
double Nue;
double KSat;
double KdAct;                 /* uM */
double INaCaX_Hill;
double INaCaX_Vmf, KmCao, KmNai, KmCai, KmNao, SLKa;    /* INaCaX */

/* parameter(s) for ISLCaPump */
double ISLCaP_Hill;
double ISLCaP_Vmf,ISLCaP_Kmf;              /* ISLCaP, the same for both junction and SL */

/*****************************Ca Handling parameters******************************/
//------------------------------------------------------------------------------
// Constants
//------------------------------------------------------------------------------
int ii;

double Ca_buffer___Bmax_Calsequestrin;   // millimolar
double Ca_buffer___Bmax_SLB_SL;   // millimolar
double Ca_buffer___Bmax_SLB_jct;   // millimolar
double Ca_buffer___Bmax_SLHigh_SL;   // millimolar
double Ca_buffer___Bmax_SLHigh_jct;   // millimolar
double Ca_buffer___koff_Calsequestrin;   // per_millisecond
double Ca_buffer___koff_SLB;   // per_millisecond
double Ca_buffer___koff_SLHigh;   // per_millisecond
double Ca_buffer___kon_Calsequestrin;   // per_millimolar_per_millisecond
double Ca_buffer___kon_SL;   // per_millimolar_per_millisecond

double Jleak_SR___KSRleak;   // per_millisecond
double Jpump_SR___H;   // dimensionless
double Jpump_SR___Kmf;   // millimolar
double Jpump_SR___Kmr;   // millimolar
double Jpump_SR___Q10_SRCaP;   // dimensionless
double Jpump_SR___V_max;   // millimolar_per_millisecond
double Jrel_SR___EC50_SR;   // millimolar
double Jrel_SR___HSR;   // dimensionless
double Jrel_SR___Max_SR;   // dimensionless
double Jrel_SR___Min_SR;   // dimensionless
double Jrel_SR___kiCa;   // per_millimolar_per_millisecond
double Jrel_SR___kim;   // per_millisecond
double Jrel_SR___koCa;   // per_millimolar2_per_millisecond
double Jrel_SR___kom;   // per_millisecond
double Jrel_SR___ks;   // per_millisecond

double cytosolic_Ca_buffer___Bmax_Calmodulin;   // millimolar
double cytosolic_Ca_buffer___Bmax_Myosin_Ca;   // millimolar
double cytosolic_Ca_buffer___Bmax_Myosin_Mg;   // millimolar
double cytosolic_Ca_buffer___Bmax_SRB;   // millimolar
double cytosolic_Ca_buffer___Bmax_TroponinC;   // millimolar
double cytosolic_Ca_buffer___Bmax_TroponinC_Ca_Mg_Ca;   // millimolar
double cytosolic_Ca_buffer___Bmax_TroponinC_Ca_Mg_Mg;   // millimolar
double cytosolic_Ca_buffer___koff_Calmodulin;   // per_millisecond
double cytosolic_Ca_buffer___koff_Myosin_Ca;   // per_millisecond
double cytosolic_Ca_buffer___koff_Myosin_Mg;   // per_millisecond
double cytosolic_Ca_buffer___koff_SRB;   // per_millisecond
double cytosolic_Ca_buffer___koff_TroponinC;   // per_millisecond
double cytosolic_Ca_buffer___koff_TroponinC_Ca_Mg_Ca;   // per_millisecond
double cytosolic_Ca_buffer___koff_TroponinC_Ca_Mg_Mg;   // per_millisecond
double cytosolic_Ca_buffer___kon_Calmodulin;   // per_millimolar_per_millisecond
double cytosolic_Ca_buffer___kon_Myosin_Ca;   // per_millimolar_per_millisecond
double cytosolic_Ca_buffer___kon_Myosin_Mg;   // per_millimolar_per_millisecond
double cytosolic_Ca_buffer___kon_SRB;   // per_millimolar_per_millisecond
double cytosolic_Ca_buffer___kon_TroponinC;   // per_millimolar_per_millisecond
double cytosolic_Ca_buffer___kon_TroponinC_Ca_Mg_Ca;   // per_millimolar_per_millisecond
double cytosolic_Ca_buffer___kon_TroponinC_Ca_Mg_Mg;   // per_millimolar_per_millisecond
double ion_diffusion___A_SL_cytosol;   // cm2
double ion_diffusion___A_jct_SL;   // cm2
double ion_diffusion___D_Ca_SL_cytosol;   // dm2_per_second
double ion_diffusion___D_Ca_jct_SL;   // dm2_per_second
double ion_diffusion___D_Na_SL_cytosol;   // dm2_per_second
double ion_diffusion___D_Na_jct_SL;   // dm2_per_second
double ion_diffusion___x_SL_cytosol;   // micrometre
double ion_diffusion___x_jct_SL;   // micrometre
double model_parameters___Cao;   // millimolar
double model_parameters___Cli;   // millimolar
double model_parameters___Clo;   // millimolar
double model_parameters___Cm_per_area;   // farad_per_cm2
double model_parameters___F;   // coulomb_per_mole
double model_parameters___Ki;   // millimolar
double model_parameters___Ko;   // millimolar
double model_parameters___Mgi;   // millimolar
double model_parameters___Nao;   // millimolar
double model_parameters___R;   // joule_per_kilomole_kelvin
double model_parameters___T;   // kelvin
double model_parameters___cell_length;   // micrometre
double model_parameters___cell_radius;   // micrometre

//------------------------------------------------------------------------------
// Computed variables
//------------------------------------------------------------------------------

double Ca_buffer___dCa_SLB_SL;   // millimolar_per_millisecond
double Ca_buffer___dCa_SLB_jct;   // millimolar_per_millisecond
double Ca_buffer___dCa_SLHigh_SL;   // millimolar_per_millisecond
double Ca_buffer___dCa_SLHigh_jct;   // millimolar_per_millisecond
double Ca_buffer___dCa_SL_tot_bound;   // millimolar_per_millisecond
double Ca_buffer___dCa_jct_tot_bound;   // millimolar_per_millisecond
double Ca_buffer___dCalsequestrin;   // millimolar_per_millisecond
double Ca_buffer___i_Ca_SL_tot;   // microA_per_microF
double Ca_buffer___i_Ca_jct_tot;   // microA_per_microF

double Jleak_SR___j_leak_SR;   // millimolar_per_millisecond
double Jpump_SR___j_pump_SR;   // millimolar_per_millisecond
double Jrel_SR___RI;   // dimensionless
double Jrel_SR___j_rel_SR;   // millimolar_per_millisecond
double Jrel_SR___kCaSR;   // dimensionless
double Jrel_SR___kiSRCa;   // per_millimolar_per_millisecond
double Jrel_SR___koSRCa;   // per_millimolar2_per_millisecond
double cell___i_Stim;   // microA_per_microF
double cytosolic_Ca_buffer___dCa_Calmodulin;   // millimolar_per_millisecond
double cytosolic_Ca_buffer___dCa_Myosin;   // millimolar_per_millisecond
double cytosolic_Ca_buffer___dCa_SRB;   // millimolar_per_millisecond
double cytosolic_Ca_buffer___dCa_TroponinC;   // millimolar_per_millisecond
double cytosolic_Ca_buffer___dCa_TroponinC_Ca_Mg;   // millimolar_per_millisecond
double cytosolic_Ca_buffer___dCa_cytosol_tot_bound;   // millimolar_per_millisecond
double cytosolic_Ca_buffer___dMg_Myosin;   // millimolar_per_millisecond
double cytosolic_Ca_buffer___dMg_TroponinC_Ca_Mg;   // millimolar_per_millisecond
double ion_diffusion___J_Ca_SL_cytosol;   // millimole_per_millisecond
double ion_diffusion___J_Ca_SL_cytosol_old;   // millimole_per_millisecond
double ion_diffusion___J_Ca_jct_SL;   // millimole_per_millisecond
double ion_diffusion___J_Ca_jct_SL_old;   // millimole_per_millisecond
double ion_diffusion___J_Na_SL_cytosol;   // millimole_per_millisecond
double ion_diffusion___J_Na_SL_cytosol_old;   // millimole_per_millisecond
double ion_diffusion___J_Na_jct_SL;   // millimole_per_millisecond
double ion_diffusion___J_Na_jct_SL_old;   // millimole_per_millisecond
double model_parameters___Cm;   // farad
double model_parameters___Vol_Cell;   // litre
double model_parameters___Vol_SL;   // litre
double model_parameters___Vol_SR;   // litre
double model_parameters___Vol_cytosol;   // litre
double model_parameters___Vol_jct;   // litre
double model_parameters___Vol_mito;   // litre

/*************State Variables*****************/
double Y[39];
double dY[39];
/************************************************/

//------------------------------------------------------------------------------
// Initialisation
//------------------------------------------------------------------------------
void init()
{
   //---------------------------------------------------------------------------
   // State variables
   //---------------------------------------------------------------------------
fCa=0.510436;

Y[0] = 1.28293;// Ca_buffer___Ca_Calsequestrin (millimolar)
Y[1] = 0.000131111;// Ca_buffer___Ca_SL (millimolar)
Y[2] = 0.0121366;// Ca_buffer___Ca_SLB_SL (millimolar)
Y[3] = 0.0116031;// Ca_buffer___Ca_SLB_jct (millimolar)
Y[4] = 0.132481; // Ca_buffer___Ca_SLHigh_SL (millimolar)
Y[5] = 0.0981466;// Ca_buffer___Ca_SLHigh_jct (millimolar)
Y[6] = 0.633156;// Ca_buffer___Ca_SR (millimolar)
Y[7] = 0.00026248;// Ca_buffer___Ca_jct (millimolar )
Y[8] = 9.93905e-05;// Ca_buffer___Cai (millimolar)
Y[11] = 0.04;
Y[23] = 2.67486e-07; // Jrel_SR___I (dimensionless)
Y[24] = 1.94911e-06;// Jrel_SR___O (dimensionless)
Y[25] = 0.879232;// Jrel_SR___R (dimensionless)
Y[26] = 8.87408;// Na_buffer___Na_SL (millimolar)
Y[27] = 0.776121;// Na_buffer___Na_SL_buf (millimolar)
Y[28] = 8.87282;// Na_buffer___Na_jct (millimolar)
Y[29] = 3.55706;// Na_buffer___Na_jct_buf (millimolar)
Y[30] = 8.87446;// Na_buffer___Nai (millimolar)
Y[31] = -85.7197; // cell___V (millivolt)
Y[32] = 0.000336009;// cytosolic_Ca_buffer___Ca_Calmodulin (millimolar)
Y[33] = 0.00244706;// cytosolic_Ca_buffer___Ca_Myosin (millimolar)
Y[34] = 0.00243042;// cytosolic_Ca_buffer___Ca_SRB (millimolar)
Y[35] = 0.00996049;// cytosolic_Ca_buffer___Ca_TroponinC (millimolar)
Y[36] = 0.12297;// cytosolic_Ca_buffer___Ca_TroponinC_Ca_Mg (millimolar)
Y[37] = 0.136961;// cytosolic_Ca_buffer___Mg_Myosin (millimolar)
Y[38] = 0.0079682; // cytosolic_Ca_buffer___Mg_TroponinC_Ca_Mg (millimolar)

NaIono = 140.0;
CaIono = 1.8;
KIono = 5.4;
KIoni = 135.0;
ClIono = 150.0;
ClIoni = 30.0;
NaIoni = 8.8;

/* initialization */
V = -84.9854;
m = 0.000662073;
h = 0.981596;
j = 0.968992;
x = 0.181464;
d = 6.49179e-07;
f = 1.00024;
f2=0.99;
xr = 8.86124e-05;
a = 0.000335844;
i = 0.999;
xto_slow = 0.0042119;
xto_fast = 0.00421165;
yto_fast = 0.994211;
yto_slow = 0.201896;
r_slow = 0.24237;

//---------------------------------------------------------------------------
// Constants
//---------------------------------------------------------------------------

GNa=8.0;
GNal=1.62e-3;
GCaL=0.44;
GKs=0.7;
GKr_factor=1.5;
GK1=0.42;
GTO_slow=0.02;
GTO_fast=0.06;
Gtos_factor=0.425;
Gtof_factor=0.85;
Gto_factor=1.0;
GCl=0.1;
GClB=2.25e-3;
GCaB=0.0002513;
GNaB=(0.297e-3*5);

INaK_Kmr= 1.5;
INaK_Hill=4.0;
INaK_Vmf=1.90719;
INaK_Kmf=11.0;
INaCaXCa_Kmf= 3.59e-3;        /* mM */
INaCaXCa_Kmr =1.3;            /* mM */
INaCaXNa_Kmf =12.29;          /* mM */
INaCaXNa_Kmr =1.0;            /* mM */
Nue= 0.35;
KSat= 0.27;
KdAct =0.256;                 /* uM */
INaCaX_Hill =3.0;
INaCaX_Vmf=9.0;
ISLCaP_Hill = 1.6;
ISLCaP_Vmf=0.0538*1.25, ISLCaP_Kmf=0.5;

Ca_buffer___Bmax_Calsequestrin = 0.14;   // millimolar
Ca_buffer___Bmax_SLB_SL = 0.0374;   // millimolar
Ca_buffer___Bmax_SLB_jct = 0.0046;   // millimolar
Ca_buffer___Bmax_SLHigh_SL = 0.0134;   // millimolar
Ca_buffer___Bmax_SLHigh_jct = 0.00165;   // millimolar
Ca_buffer___koff_Calsequestrin = 65.0;   // per_millisecond
Ca_buffer___koff_SLB = 1.3;   // per_millisecond
Ca_buffer___koff_SLHigh = 30.0e-3;   // per_millisecond
Ca_buffer___kon_Calsequestrin = 100.0;   // per_millimolar_per_millisecond
Ca_buffer___kon_SL = 100.0;   // per_millimolar_per_millisecond



Jleak_SR___KSRleak = 5.348e-6;   // per_millisecond
Jpump_SR___H = 1.787;   // dimensionless
Jpump_SR___Kmf = 0.000246;   // millimolar
Jpump_SR___Kmr = 1.7;   // millimolar
Jpump_SR___Q10_SRCaP = 2.6;   // dimensionless
Jpump_SR___V_max = 286.0e-6;   // millimolar_per_millisecond
Jrel_SR___EC50_SR = 0.45;   // millimolar
Jrel_SR___HSR = 2.5;   // dimensionless
Jrel_SR___Max_SR = 15.0;   // dimensionless
Jrel_SR___Min_SR = 1.0;   // dimensionless
Jrel_SR___kiCa = 0.5;   // per_millimolar_per_millisecond
Jrel_SR___kim = 0.005;   // per_millisecond
Jrel_SR___koCa = 10.0;   // per_millimolar2_per_millisecond
Jrel_SR___kom = 0.06;   // per_millisecond
Jrel_SR___ks = 25.0;   // per_millisecond

cytosolic_Ca_buffer___Bmax_Calmodulin = 0.024;   // millimolar
cytosolic_Ca_buffer___Bmax_Myosin_Ca = 0.14;   // millimolar
cytosolic_Ca_buffer___Bmax_Myosin_Mg = 0.14;   // millimolar
cytosolic_Ca_buffer___Bmax_SRB = 0.0171;   // millimolar
cytosolic_Ca_buffer___Bmax_TroponinC = 0.07;   // millimolar
cytosolic_Ca_buffer___Bmax_TroponinC_Ca_Mg_Ca = 0.14;   // millimolar
cytosolic_Ca_buffer___Bmax_TroponinC_Ca_Mg_Mg = 0.14;   // millimolar
cytosolic_Ca_buffer___koff_Calmodulin = 238.0e-3;   // per_millisecond
cytosolic_Ca_buffer___koff_Myosin_Ca = 0.46e-3;   // per_millisecond
cytosolic_Ca_buffer___koff_Myosin_Mg = 0.057e-3;   // per_millisecond
cytosolic_Ca_buffer___koff_SRB = 60.0e-3;   // per_millisecond
cytosolic_Ca_buffer___koff_TroponinC = 19.6e-3;   // per_millisecond
cytosolic_Ca_buffer___koff_TroponinC_Ca_Mg_Ca = 0.032e-3;   // per_millisecond
cytosolic_Ca_buffer___koff_TroponinC_Ca_Mg_Mg = 3.33e-3;   // per_millisecond
cytosolic_Ca_buffer___kon_Calmodulin = 34.0;   // per_millimolar_per_millisecond
cytosolic_Ca_buffer___kon_Myosin_Ca = 13.8;   // per_millimolar_per_millisecond
cytosolic_Ca_buffer___kon_Myosin_Mg = 15.7e-3;   // per_millimolar_per_millisecond
cytosolic_Ca_buffer___kon_SRB = 100.0;   // per_millimolar_per_millisecond
cytosolic_Ca_buffer___kon_TroponinC = 32.7;   // per_millimolar_per_millisecond
cytosolic_Ca_buffer___kon_TroponinC_Ca_Mg_Ca = 2.37;   // per_millimolar_per_millisecond
cytosolic_Ca_buffer___kon_TroponinC_Ca_Mg_Mg = 3.0e-3;   // per_millimolar_per_millisecond
ion_diffusion___A_SL_cytosol = 1.3e-4;   // cm2
ion_diffusion___A_jct_SL = 3.01e-6;   // cm2
ion_diffusion___D_Ca_SL_cytosol = 1.22e-6;   // dm2_per_second
ion_diffusion___D_Ca_jct_SL = 1.64e-6;   // dm2_per_second
ion_diffusion___D_Na_SL_cytosol = 1.79e-5;   // dm2_per_second
ion_diffusion___D_Na_jct_SL = 1.09e-5;   // dm2_per_second
ion_diffusion___x_SL_cytosol = 0.45;   // micrometre
ion_diffusion___x_jct_SL = 0.5;   // micrometre
model_parameters___Cao = 1.8;   // millimolar
model_parameters___Cli = 15.0;   // millimolar
model_parameters___Clo = 150.0;   // millimolar
model_parameters___Cm_per_area = 2.0e-6;   // farad_per_cm2
model_parameters___F = 96486.7;   // coulomb_per_mole
model_parameters___Ki = 135.0;   // millimolar
model_parameters___Ko = 5.4;   // millimolar
model_parameters___Mgi = 1.0;   // millimolar
model_parameters___Nao = 140.0;   // millimolar
model_parameters___R = 8314.3;   // joule_per_kilomole_kelvin
model_parameters___T = 310.0;   // kelvin
model_parameters___cell_length = 100.0;   // micrometre
model_parameters___cell_radius = 10.25;   // micrometre

//---------------------------------------------------------------------------
// Computed variables
//---------------------------------------------------------------------------

model_parameters___Vol_Cell = 3.141592654*pow(model_parameters___cell_radius/1000.0, 2.0)*model_parameters___cell_length/pow(1000.0, 3.0);
model_parameters___Vol_cytosol = 0.65*model_parameters___Vol_Cell;
model_parameters___Vol_SR = 0.035*model_parameters___Vol_Cell;
model_parameters___Vol_SL = 0.02*model_parameters___Vol_Cell;
model_parameters___Vol_jct = 0.00051*model_parameters___Vol_Cell;
model_parameters___Cm = model_parameters___Cm_per_area*2.0*model_parameters___cell_radius/10000.0*3.14159265358979*model_parameters___cell_length/10000.0;

}

/*****************Current calculatations********************/

void comp_ina(){
//labheart

alpha_m = 0.32 * (V+47.13)/(1- exp (-0.1* (V+47.13)));
beta_m = 0.08 * exp (-V/11);
if ( V >= -40 ){
alpha_h = 0;
beta_h = 1/(0.13 * (1+exp ((V+10.66)/-11.1)));
alpha_j = 0;
beta_j = 0.3 * exp (-2.535 * 10e-7 * V)/(1+ exp (-0.1* (V+32)));
}
else if ( V < -40){
alpha_h = 0.135 * exp ((80+V)/-6.8);
beta_h = 3.56 * exp (0.079 * V) + 3.1* 10e5 * exp (0.35*V);
alpha_j = ((-1.2714 * 10e5 * exp (0.2444* V) - (3.474 * 10e-5 * exp (-0.0439 * V)) * (V+37.78)))/(1+exp (0.311 * (V+79.23)));
beta_j = 0.1212 * exp (-0.01052* V)/(1+exp (-0.1378 * (V+40.14)));
}
tau_m	= 1.0/(alpha_m+beta_m);
tau_h	= 1.0/(alpha_h+beta_h);
tau_j	= 1.0/(alpha_j+beta_j);		
m_inf	= alpha_m*tau_m;
h_inf	= alpha_h*tau_h;
j_inf	= alpha_j*tau_j;

m += HT*(m_inf - m)/tau_m;
h += HT*(h_inf - h)/tau_h;
j += HT*(j_inf - j)/tau_j;

INa =GNa * pow(m, 3.0) *h*j* (V-ENa);

//Slow component

alpha_ml= 0.32*(V+47.13)/(1-exp(-0.1*(V+47.13)));
beta_ml	= 0.08*exp(-V/11);
tau_ml	= 1/(alpha_ml + beta_ml);
ml_inf	= alpha_ml/(alpha_ml+beta_ml);

hl_inf=1/(1+exp((V+80.5)/6.3));
tau_hl=132.4+ 112.8*exp(0.02325*V);

ml += HT*(ml_inf - ml)/tau_ml;
hl += HT*(hl_inf - hl)/tau_hl;

INal=GNal*ml*hl*(V-ENa);
}

void comp_ica_ipca(){

d_inf = 1.0/(1.0 + exp(-1.0*(V-19 + 27.5)/4.0));
tau_d=0.4*(1.0/(1.0+exp(-1.0*(V-19 + 27.5)/4.0)))*(1-exp(-1.0*(V-19 + 27.5)/4.0))/(0.035*(V-19 + 27.5));
f_inf = 1.0/(1.0 + exp((V -14+ 42.06)/6.0));
tau_f2 = (5+ 30.0/(1.0 + exp(-1.0*(V-30)/5)))+55;
tau_f1 = 8+ (20.0/(1.0 + exp(-1.0*(V-20)/5))-20.0/(1.0 + exp(-1.0*(V-40)/5)));

d += HT*(d_inf - d)/tau_d;
f += HT*(f_inf - f)/tau_f1;
f2 += HT*(f_inf - f2)/tau_f2;

ICa = GCaL*d*(0.8*f+0.2*f2)*(1-Y[11])*(V-60);
}

void comp_iks(){

pCa = -1.0*log10(Y[8])+3.0;
SLGKs = 2.0*0.07*(0.057+0.19/(1.0+exp((-7.2+pCa)/0.6)));

x_inf = 1.0/(1.0+exp(-1.0*(V-1.5)/13));
tau_x = 0.5*(600.0/(1.0+exp((V-20)/15.0)) + 250.0);
x += HT*(x_inf - x)/tau_x;

IKs = GKs* SLGKs * x * x * (V-EKs);

}

void comp_ikr(){
                                                          /* IKr */
V -= 35.0; 
xr_inf = 1.0/(1.0 + exp(-1.0 * (V + 50.0)/7.5));
tau_xr = 1.0/(0.00138 * (V + 7.0)/(1.0 - exp(-0.123 * (V + 7.0))) + 0.00061 * (V + 10.0)/(exp(0.145 * (V + 10.0)) - 1.0));
V += 35.0;
r_inf=1/(1+6*exp(0.05*V));

V -= 10.0;
tau_r = 200.0/(1.0 + exp((V + 40.0)/20.0)) + 200.0;
GKr = 0.03 * sqrt(KIono/5.4);  
V += 10.0;

xr += HT*(xr_inf - xr)/tau_xr;

IKr =  GKr_factor*GKr * xr * r_inf * (V - EK);

}

void comp_ik1(){

V -= 5.0; 
alpha_k1 = 1.02/(1.0 + exp(0.2385 * (V - EK - 59.215)));
beta_k1 = (0.49124 * exp(0.08032 * (V - EK + 5.476)) + exp(0.06175 * (V - EK - 594.31)))/(1.0 + exp(-0.5143 * (V - EK + 4.753)));
k1_inf = alpha_k1/(alpha_k1 + beta_k1);
V += 5.0;

IK1 = GK1 * k1_inf * (V - EK);

}
void comp_ikp(){
IKp_Kp = 1.0/(1.0 + exp((7.488 - V)/5.98));
IKp = 0.001 * IKp_Kp * (V - EK);
}

void comp_ito(){

xto_slow_inf = 1.0 / (1.0 + exp((V + 3.0)/-15.0));
yto_slow_inf = 1.0 / (1.0 + exp((V + 33.5)/10.0));
r_slow_inf = 1.0 / (1.0 + exp((V + 33.5)/10.0));
tau_xto_slow = 9.0 / (1.0 + exp((V + 3.0)/15.0)) + 0.5;
tau_yto_slow = 3000.0 / (1.0 + exp((V + 60.0)/10.0)) + 30.0;
tau_r_slow = 2800.0 / (1.0 + exp((V + 60.0)/10.0)) + 220.0;

xto_fast_inf = 1.0 / (1.0 + exp((V + 3.0)/-15.0));
yto_fast_inf = 1.0 / (1.0 + exp((V + 33.5)/10.0));
tau_xto_fast = 3.5 * exp(-1.0*(V / 30.0) * (V / 30.0)) + 1.5;
tau_yto_fast = 20.0 / (1.0 + exp((V + 33.5)/10.0)) + 20.0;

xto_slow += HT*(xto_slow_inf - xto_slow)/tau_xto_slow;
yto_slow += HT*(yto_slow_inf - yto_slow)/tau_yto_slow;
r_slow += HT*(r_slow_inf - r_slow)/tau_r_slow;
xto_fast += HT*(xto_fast_inf - xto_fast)/tau_xto_fast;
yto_fast += HT*(yto_fast_inf - yto_fast)/tau_yto_fast;


ITO_slow = Gtos_factor*  GTO_slow * xto_slow * (yto_slow + 0.5 * r_slow_inf) * (V - EK);
ITO_fast = Gtof_factor * GTO_fast * xto_fast * yto_fast * (V - EK);  

ITO = Gto_factor*(ITO_slow + ITO_fast);

}

void comp_icl(){

a_inf = 1.0/(1.0 + exp(-(V + 5.0)/10.0));
i_inf = 1.0/(1.0 + exp((V + 75.0)/10.0));
tau_i = 10.0 / (1.0 + exp((V + 33.5)/10.0)) + 10.0;

a += HT*(a_inf - a)/1.0;
i += HT*(i_inf - i)/tau_i;

ICl = GCl * a_inf * i / (1.0 + 0.1/(Y[8]*1000)) * (V - ECl);

}

void comp_ib(){

IClB = GClB * (V - ECl);
ICaB = GCaB * (V - ECa);
INaB = GNaB * (V - ENa);
}

void comp_inak(){

sigma = (exp(NaIono/67.3) - 1.0)/7.0;
fNaK = 1.0/(1.0 + 0.1245 * exp(-0.1 * V/(RTF)) + 0.0365 * sigma * exp(-1.0*V/(RTF)));
INaK = INaK_Vmf * fNaK * KIono /(1.0 + pow((INaK_Kmf/NaIoni), INaK_Hill))/(KIono + INaK_Kmr);
}

void comp_inaca(){

SLKa = 1.0 / (1.0 + (KdAct / (Y[8]*1000)) * (KdAct / (Y[8]*1000)) * (KdAct / (Y[8]*1000)));
KmCao = INaCaXCa_Kmr;
KmNai = INaCaXNa_Kmf;
KmCai = INaCaXCa_Kmf;
KmNao = 87.5;              

INaCaX = SLKa * INaCaX_Vmf * ((pow(NaIoni, INaCaX_Hill) * CaIono * exp(Nue * V/(RTF))) - pow(NaIono, INaCaX_Hill) * Y[8] * exp((Nue-1.0) * V/(RTF)))/((KmCao*pow(NaIoni, INaCaX_Hill)+pow(KmNao,INaCaX_Hill)*Y[8] +pow(KmNai,INaCaX_Hill)* CaIono *(1.0+Y[8] /KmCai)+KmCai*pow(NaIono,INaCaX_Hill)*(1.0+pow((NaIoni/KmNai),INaCaX_Hill))+pow(NaIoni,INaCaX_Hill)* CaIono +pow(NaIono,INaCaX_Hill)*Y[8])*(1.0+KSat*exp((Nue-1.0) * V/(RTF))));//taken from shannon bers 04 code

}

void comp_ip(){

ISLCaP = ISLCaP_Vmf/(1.0 + pow((ISLCaP_Kmf/(Y[8]*1000)), ISLCaP_Hill));
}

void comp_itot()
{
Iion = (-(IKs+IKr+IKp+IK1+ITO+ICl+IClB+INa+INal+ICa+INaB+ICaB+INaK+INaCaX+ISLCaP));
}

void comp_cai ()
{

   Ca_buffer___dCalsequestrin = Ca_buffer___kon_Calsequestrin*Y[6]*(Ca_buffer___Bmax_Calsequestrin*model_parameters___Vol_cytosol/model_parameters___Vol_SR-Y[0])-	Ca_buffer___koff_Calsequestrin*Y[0];

   dY[0] = Ca_buffer___dCalsequestrin;

   Ca_buffer___dCa_SLB_SL = Ca_buffer___kon_SL*Y[1]*(Ca_buffer___Bmax_SLB_SL*model_parameters___Vol_cytosol/model_parameters___Vol_SL-Y[2])-Ca_buffer___koff_SLB*Y[2];
   Ca_buffer___dCa_SLB_jct = Ca_buffer___kon_SL*Y[7]*(Ca_buffer___Bmax_SLB_jct*0.1*model_parameters___Vol_cytosol/model_parameters___Vol_jct-Y[3])-Ca_buffer___koff_SLB*Y[3];
   Ca_buffer___dCa_SLHigh_SL = Ca_buffer___kon_SL*Y[1]*(Ca_buffer___Bmax_SLHigh_SL*model_parameters___Vol_cytosol/model_parameters___Vol_SL-Y[4])-Ca_buffer___koff_SLHigh*Y[4];
   Ca_buffer___dCa_SLHigh_jct = Ca_buffer___kon_SL*Y[7]*(Ca_buffer___Bmax_SLHigh_jct*0.1*model_parameters___Vol_cytosol/model_parameters___Vol_jct-Y[5])-Ca_buffer___koff_SLHigh*Y[5];

   dY[2] = Ca_buffer___dCa_SLB_SL;
   dY[3] = Ca_buffer___dCa_SLB_jct;
   dY[4] = Ca_buffer___dCa_SLHigh_SL;
   dY[5] = Ca_buffer___dCa_SLHigh_jct;

   Ca_buffer___dCa_jct_tot_bound = Ca_buffer___dCa_SLB_jct+Ca_buffer___dCa_SLHigh_jct;
   Ca_buffer___dCa_SL_tot_bound = Ca_buffer___dCa_SLB_SL+Ca_buffer___dCa_SLHigh_SL;


   Ca_buffer___i_Ca_jct_tot = ICa;
   Ca_buffer___i_Ca_SL_tot = -2.0*INaCaX+ICaB+ISLCaP;

   Jpump_SR___j_pump_SR = Jpump_SR___V_max*model_parameters___Vol_cytosol/model_parameters___Vol_SR*(pow(Y[8]/Jpump_SR___Kmf, Jpump_SR___H)-pow(Y[6]/Jpump_SR___Kmr, Jpump_SR___H))/(1.0+pow(Y[8]/Jpump_SR___Kmf, Jpump_SR___H)+pow(Y[6]/Jpump_SR___Kmr, Jpump_SR___H));
   Jleak_SR___j_leak_SR = Jleak_SR___KSRleak*(Y[6]-Y[7]);
   Jrel_SR___j_rel_SR = Jrel_SR___ks*Y[24]*(Y[6]-Y[7]);

   dY[6] = Jpump_SR___j_pump_SR-(Jleak_SR___j_leak_SR*model_parameters___Vol_cytosol/model_parameters___Vol_SR+Jrel_SR___j_rel_SR)-Ca_buffer___dCalsequestrin;
   ion_diffusion___J_Ca_jct_SL = (Y[7]-Y[1])*8.2413e-13;

   dY[7] = -0.65*Ca_buffer___i_Ca_jct_tot*model_parameters___Cm/(model_parameters___Vol_jct*2.0*model_parameters___F)-ion_diffusion___J_Ca_jct_SL/model_parameters___Vol_jct+Jrel_SR___j_rel_SR*model_parameters___Vol_SR/model_parameters___Vol_jct+Jleak_SR___j_leak_SR*model_parameters___Vol_cytosol/model_parameters___Vol_jct-1.0*Ca_buffer___dCa_jct_tot_bound;

   ion_diffusion___J_Ca_SL_cytosol = (Y[1]-Y[8])*3.7243e-12;
   dY[1] = -0.65*Ca_buffer___i_Ca_SL_tot*model_parameters___Cm/(model_parameters___Vol_SL*2.0*model_parameters___F)+(ion_diffusion___J_Ca_jct_SL-ion_diffusion___J_Ca_SL_cytosol)/model_parameters___Vol_SL-1.0*Ca_buffer___dCa_SL_tot_bound;

   cytosolic_Ca_buffer___dCa_TroponinC = cytosolic_Ca_buffer___kon_TroponinC*Y[8]*(cytosolic_Ca_buffer___Bmax_TroponinC-Y[35])-cytosolic_Ca_buffer___koff_TroponinC*Y[35];
   cytosolic_Ca_buffer___dCa_TroponinC_Ca_Mg = cytosolic_Ca_buffer___kon_TroponinC_Ca_Mg_Ca*Y[8]*(cytosolic_Ca_buffer___Bmax_TroponinC_Ca_Mg_Ca-(Y[36]+Y[38]))-cytosolic_Ca_buffer___koff_TroponinC_Ca_Mg_Ca*Y[36];
   cytosolic_Ca_buffer___dMg_TroponinC_Ca_Mg = cytosolic_Ca_buffer___kon_TroponinC_Ca_Mg_Mg*model_parameters___Mgi*(cytosolic_Ca_buffer___Bmax_TroponinC_Ca_Mg_Mg-(Y[36]+Y[38]))-cytosolic_Ca_buffer___koff_TroponinC_Ca_Mg_Mg*Y[38];
   cytosolic_Ca_buffer___dCa_Calmodulin = cytosolic_Ca_buffer___kon_Calmodulin*Y[8]*(cytosolic_Ca_buffer___Bmax_Calmodulin-Y[32])-cytosolic_Ca_buffer___koff_Calmodulin*Y[32];
   cytosolic_Ca_buffer___dCa_Myosin = cytosolic_Ca_buffer___kon_Myosin_Ca*Y[8]*(cytosolic_Ca_buffer___Bmax_Myosin_Ca-(Y[33]+Y[37]))-cytosolic_Ca_buffer___koff_Myosin_Ca*Y[33];
   cytosolic_Ca_buffer___dMg_Myosin = cytosolic_Ca_buffer___kon_Myosin_Mg*model_parameters___Mgi*(cytosolic_Ca_buffer___Bmax_Myosin_Mg-(Y[33]+Y[37]))-cytosolic_Ca_buffer___koff_Myosin_Mg*Y[37];
   cytosolic_Ca_buffer___dCa_SRB = cytosolic_Ca_buffer___kon_SRB*Y[8]*(cytosolic_Ca_buffer___Bmax_SRB-Y[34])-cytosolic_Ca_buffer___koff_SRB*Y[34];
   cytosolic_Ca_buffer___dCa_cytosol_tot_bound = cytosolic_Ca_buffer___dCa_TroponinC+cytosolic_Ca_buffer___dCa_TroponinC_Ca_Mg+cytosolic_Ca_buffer___dMg_TroponinC_Ca_Mg+cytosolic_Ca_buffer___dCa_Calmodulin+cytosolic_Ca_buffer___dCa_Myosin+cytosolic_Ca_buffer___dMg_Myosin+cytosolic_Ca_buffer___dCa_SRB;

   dY[8] = -Jpump_SR___j_pump_SR*model_parameters___Vol_SR/model_parameters___Vol_cytosol+ion_diffusion___J_Ca_SL_cytosol/model_parameters___Vol_cytosol-1.0*cytosolic_Ca_buffer___dCa_cytosol_tot_bound;
	dY[11] = 0.275*Y[7]*(1.0-Y[11])-2.9e-3*Y[11];

   Jrel_SR___kCaSR = Jrel_SR___Max_SR-(Jrel_SR___Max_SR-Jrel_SR___Min_SR)/(1.0+pow(Jrel_SR___EC50_SR/Y[6], Jrel_SR___HSR));

   Jrel_SR___koSRCa = Jrel_SR___koCa/Jrel_SR___kCaSR;
   Jrel_SR___kiSRCa = Jrel_SR___kiCa*Jrel_SR___kCaSR;
   Jrel_SR___RI = 1.0-Y[25]-Y[24]-Y[23];

   dY[25] = Jrel_SR___kim*Jrel_SR___RI-Jrel_SR___kiSRCa*Y[7]*Y[25]-(Jrel_SR___koSRCa*Y[7]*Y[7]*Y[25]-Jrel_SR___kom*Y[24]);
   dY[24] = Jrel_SR___koSRCa*Y[7]*Y[7]*Y[25]-Jrel_SR___kom*Y[24]-(Jrel_SR___kiSRCa*Y[7]*Y[24]-Jrel_SR___kim*Y[23]);
   dY[23] = Jrel_SR___kiSRCa*Y[7]*Y[24]-Jrel_SR___kim*Y[23]-(Jrel_SR___kom*Y[23]-Jrel_SR___koSRCa*Y[7]*Y[7]*Jrel_SR___RI);


   dY[35] = cytosolic_Ca_buffer___dCa_TroponinC;
   dY[36] = cytosolic_Ca_buffer___dCa_TroponinC_Ca_Mg;
   dY[38] = cytosolic_Ca_buffer___dMg_TroponinC_Ca_Mg;
   dY[32] = cytosolic_Ca_buffer___dCa_Calmodulin;
   dY[33] = cytosolic_Ca_buffer___dCa_Myosin;
   dY[37] = cytosolic_Ca_buffer___dMg_Myosin;
   dY[34] = cytosolic_Ca_buffer___dCa_SRB;

   for (ii = 0; ii <= 38; ii++)
	 Y[ii] += HT*dY[ii];
}

int main()

{
double cl;
char option1,option2;
double max_bcl_restitution;

ofstream out("ap.txt");
ofstream out1("restitution.txt");
ofstream out2("initials_gates.txt");
ofstream out3("initials_cai.txt");
ofstream out4("ap_char.txt");
ofstream out10("ap_ina.txt");

int flag = 1;

    while ( flag==1){

        cout << "would you like to run restitution simulation? (Y/N) " ;
        cin >> option1;
        cout << endl << endl;
        if (option1== 'Y' || option1 == 'y'){
        max_bcl_restitution=2000;
        flag = 0;
        }

        else if (option1== 'n' || option1 == 'N'){
        max_bcl_restitution=0;
        flag =0;
        }
        else {flag=1;}
        }

cout << "insert cycle length to start (ms): ";
cin >> cl;
cout << endl << endl;
for (int loop=0;loop<=max_bcl_restitution;loop=loop+100){

cout << "Simulation @ bcl: " << loop + cl << endl << endl;
bcl=cl+loop;

//if (bcl > 600){
//loop=loop+75;
//}


double STOPTIME=1000;
double tstim=10.0;
double tend=11.0;
double stimtime=0.0;
double time=0;
double dvdt;
int step;
int i=-1;

double vmax[beats+1];
double dvdtmax[beats+1];
double apd[beats+1];
double toneapd[beats+1];
double ttwoapd[beats+1];
double trep[beats+1];
double di[beats+1];
double rmbp[beats+1];

init();


ENa =(RTF)*log(NaIono/NaIoni);
EKs = (RTF)*log((KIono + PNaK * NaIono)/(KIoni + PNaK * NaIoni));
ECa = (RTF)/2.0 * log(CaIono/(Y[8]));
EK = (RTF)*log(KIono/KIoni);
ECl = (RTF)*log(ClIoni/ClIono);

for (i = 0; i <= beats; i++) {
  dvdtmax[i] = -10000;
  vmax[i] = -10000;
} i = -1;


  STOPTIME=S2+bcl*beats;
  for(step=0;time<=STOPTIME;step++)
    {
      time+=HT;

      if (time>=tstim && time<(tstim+HT) && time <= (S2 + bcl*(beats-1))) {
       stimtime = 0;
       i = i+1;

       if (i < beats-2)
         tstim = tstim+bcl;
       else if (i == beats-2)
         tstim = tstim+S2;
       else tstim = tstim+10000;
       cout << "Stimulus " << i+1 <<" applied, Time = " << time << endl;

       rmbp[i]=V;
     }

comp_iks();
comp_ikr();
comp_ik1();
comp_ikp();
comp_ito();
comp_icl();
comp_ib();
comp_ina();
comp_ica_ipca();
comp_inak();
comp_inaca();
comp_ip();
comp_itot();
comp_cai();

if(stimtime>=0 && stimtime<1)
Istim=-40;
else Istim=0.0;

stimtime = stimtime+HT;

dvdt = Iion-Istim;
V += HT*dvdt;

if (V>vmax[i])
vmax[i] = V;
if (dvdt>dvdtmax[i])
{dvdtmax[i] = dvdt;
toneapd[i] = time;}
if (V>=(vmax[i]-0.9*(vmax[i]-rmbp[i])))
ttwoapd[i] = time;

if (step%200 == 0)
out<<time - 10<< "\t"<< V <<  "\t"<< Y[8]<<  "\t"<< INa<<  "\t"<< INa*0 << "\t"<< IKr<<  "\t"<<IKs<<  "\t"<< IK1<< "\t"<< IK1*0<< "\t"<< INaB<<  "\t"<<ICaB<<  "\t"<< ICa<< "\t"<< ICa*0<<"\t"<<INaK<<"\t"<<INaCaX<<"\t"<<ISLCaP<<"\t"<<ITO<<"\t"<<ICl<<endl;

//if (step%2 == 0)
//out10<<time << "\t"<< V <<  "\t"<< INa<<endl;

    } //end of time loop


out2<< "V = " << V << "\n"<< "m = " << m << "\n"<< "h = " << h << "\n"<< "j = " << j << "\n"<< "x = " << x << "\n"<< "d = " << d << "\n"
	<< "f = " << f << "\n"<< "xr = " << xr << "\n"<< "a = " << a<< "\n"<< "i = " << i<< "\n"<< "xto_slow = " <<xto_slow << "\n"<< "xto_fast = " << xto_fast << "\n"<< "yto_fast = " << yto_fast<< "\n"
	<< "yto_slow = " << yto_slow << "\n"<< "r_slow = " << r_slow<< endl;

out3<< "Y[0] = " << Y[0] << "\n"<< "Y[1] = " << Y[1] << "\n"<< "Y[2] = " << Y[2] << "\n"<< "Y[3] = " << Y[3] << "\n"<< "Y[4] = " << Y[4] << "\n"<< "Y[5] = " << Y[5] << "\n"
	<< "Y[6] = " << Y[6] << "\n"<< "Y[7] = " << Y[7] << "\n"<< "Y[8] = " << Y[8]<< "\n"<< "Y[23] = " << Y[23]<< "\n"<< "Y[24] = " <<Y[24] << "\n"<< "Y[25] = " << Y[25] << "\n"<< "Y[26] = " << Y[26]<< "\n"
	<< "Y[27] = " << Y[27] << "\n"<< "Y[28] = " << Y[28] << "\n"<< "Y[29] = " << Y[29] << "\n"<< "Y[30] = " << Y[30] << "\n"<< "Y[31] = " << Y[31] << "\n"<< "Y[32] = " << Y[32] << "\n"<< "Y[33] = " << Y[33] << "\n"<< "Y[34] = " << Y[34] << "\n"<< "Y[35] = " << Y[35] << "\n"<< "Y[36] = " << Y[36] << "\n"<< "Y[37] = " << Y[37] << "\n"<< "Y[38] = " << Y[38] <<"\t"<<fCa<< endl;

apd[beats-1] = ttwoapd[beats-1]-toneapd[beats-1];
apd[beats-2] = ttwoapd[beats-2]-toneapd[beats-2];
di[beats-1] = toneapd[beats-1]-ttwoapd[beats-2];

cout << "DI = " << di[beats-1] << "\t" << "APD90 = " << apd[beats -1]<< "\t" << apd[beats-2]<< endl;
out1 << bcl << "\t" << apd[beats -1]<< "\t" << apd[beats-2]<<endl;
out4 << bcl << "\t" << apd[beats -1]<< "\t" << dvdtmax[beats-1]<< "\t" << rmbp[beats -1]<< "\t" << vmax[beats-1]-rmbp[beats-1]<<endl;


}return 0;

}