/*********************************************************************
shiffox.c

  Program integrates combination of 

  Yohannes Shiferaw's calcium model 
  (Y. Shiferaw, M. A. Watanabe, A. Garfinkel, J. N. Weiss, and A. Karma,
  Model of Intracellular Calcium Cycling in Ventricular Myocytes, 
  Biophys J 2003; 85: 3666-3686) 

  and

  Jeff Fox's canine ventricular cell ionic model 
  (J. J. Fox, J. L. McHarg, and R. F. Gilmour Jr,
  Ionic mechanism of electrical alternans,
  Am J Physiol Heart Circ Physiol 2002; 282: 516-530) 

  according to Yohannes Shiferaw 
  (personal communications and
  Y. Shiferaw, D. Sato, and A. Karma,
  Coupled dynamics of voltage and calcium in paced cardiac cells,
  Phys Rev E 2005; 71: 021903).


  Trine Krogh-Madsen, October 2004

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

#include <math.h>
#include <stdio.h>


const double Bsr = 47.0;
const double Ksr = 0.6;
const double Bcd = 24.0;
const double Kcd = 7.0;
const double k_on = 37.7*0.001;
const double k_off = 19.6*0.001;
const double KT = 0.5199; //koff/kon;
const double BT = 70.0;

const double vivs = 10.0;
const double ts = 2.0;
const double tr = 20.0;
const double ta = 50.0;
const double g = 250.0*0.3;

const double F = 96490.0;
const double F2 = 96.490;
const double R = 8.315;
const double Temp = 308.0;
const double C_sc = 1.0;

const double Ca_o = 1.8;
const double Ica_bar = 3.5;
const double Pca = 0.00000054;
const double td = 5.0;
const double tf = 30.0;
const double gama = 0.2;
const double tq = 20.0;

const double cup = 0.5;
const double vup = 250.0*0.001;
const double INACA = 8.0;
const double Km_Na = 87.5;
const double Km_Ca = 1.380;
const double ksat = 0.1;
const double xi = 0.35;

const double Na_i = 10.0;
const double Na_o = 140.0;
const double K_i = 149.4;
const double K_o = 4.0;

const double g_Na = 12.0;
const double g_K1 = 2.8;
const double Km_K1 = 13.0;
const double g_Kr = 0.0136;
const double g_Ks = 0.0245;
const double g_Kp = 0.002216;
const double g_to = 0.1;
const double cst = 1.0;

const double u = 9.0;

const double dt = 0.01;
const double T = 340.;   //BCL
const double t_end = 34000.; 
const double stim_dur = 1.0;
const double stim_amp = -50.0;
const int N_BCL = 34000;  // BCL/dt
const int N_stim = 100;    // stim_dur/dt



/******* Shiferaw functions ******/

double Q(double cjp) {
  double s, Q;
  s = 1.5*(110.0-50.0)-u*110.0;
  if (cjp < 50.0) Q = 0.0;
  else if (cjp >= 110.0) Q = u*cjp+s;
  else Q = 1.5*(cjp-50.0);
  Q*=(0.001);
  return Q;
}

double beta_ci(double ci) {
  return 1.0/(1.0+(Bsr*Ksr/((ci+Ksr)*(ci+Ksr)))+(Bcd*Kcd/((ci+Kcd)*(ci+Kcd)))+(BT*KT/((ci+KT)*(ci+KT))));
}

double beta_cs(double cs) {
  return 1.0/(1.0+(Bsr*Ksr/((cs+Ksr)*(cs+Ksr)))+(Bcd*Kcd/((cs+Kcd)*(cs+Kcd)))+(BT*KT/((cs+KT)*(cs+KT))));
}

double INaCa(double cs,double V) {
  V/=1000.0;
  cs/=1000.0;
  double a = V*F/(R*Temp);
  return INACA/(Km_Na*Km_Na*Km_Na+Na_o*Na_o*Na_o)/(Km_Ca+Ca_o)*(exp(xi*a)*Na_i*Na_i*Na_i*Ca_o-exp((xi-1.0)*a)*Na_o*Na_o*Na_o*cs)/(1.0+ksat*exp((xi-1.0)*a));
}

double Iup(double ci) {
  return vup*ci*ci/(ci*ci+cup*cup);
}

double qinf(double cs) {
  return 1.0/(1.0+pow((cs/cst),gama));
}

double Ica(double d, double f, double q, double cs, double V) {
  V/=1000.0;
  cs/=1000.0;
  return d*f*q*Ica_bar*Pca*(4.0*V*F*F/R/Temp)*(cs*exp(2.0*V*F/R/Temp)-0.341*Ca_o)/(exp(2.0*V*F/R/Temp)-1.0);
}


int main() {


  /* files */  
  FILE *output;
  output = fopen("shiffox.dat","w");


  /* variables */
  double t = 0.0;
  double m, h, j, Xr, Xs, Xto, Yto, d, f, q, cs[2], ci[2], cj[2], cjp[2], Irel[2], V;
  double E_Na, E_K, K1_inf, RV, E_Ks, K_Kp;
  double I_Na, I_K1, I_Kr, I_Ks, I_to, I_Kp, I_NaCa, I_Ca, I_stim, I_total;
  double V_min = -120.0;
  int ilow, tl = 1800;
  double linext;
  double alpha_m, beta_m, alpha_h, beta_h, alpha_j, beta_j, Xr_inf, Xr_tau, Xs_inf, Xs_tau, alpha_Xto, beta_Xto, alpha_Yto, beta_Yto, d_inf, f_inf;
  double tau_m, tau_h, tau_j, tau_Xto, tau_Yto;
  double v;
  double vt[tl], alpha_mt[tl], beta_mt[tl], alpha_ht[tl], beta_ht[tl], alpha_jt[tl], beta_jt[tl], Xr_inft[tl], Xr_taut[tl], Xs_inft[tl], Xs_taut[tl], alpha_Xtot[tl], beta_Xtot[tl], alpha_Ytot[tl], beta_Ytot[tl], d_inft[tl], f_inft[tl];


  int n = 0;
  int i = 0;
  int n_BCL = 1;
  int n_stim = 1;
  int stim_on = 0;

  
  /* make look-up tables */
  for (i=0; i<tl; i++) {
    v = V_min+0.1*i;
    vt[i] = V_min+0.1*i;
    alpha_mt[i] = 0.32*(v+47.13)/(1.0-exp(-0.1*(v+47.13)));
    beta_mt[i] = 0.08*exp(-v/11.0);
    alpha_ht[i] = 0.135*exp(-(v+80.0)/6.8);
    beta_ht[i] = 7.5/(1.0+exp(-0.1*(v+11.0)));
    alpha_jt[i] = (0.175*exp(-(v+100.0)/23.0))/(1.0+exp(0.15*(v+79.0)));
    beta_jt[i] = 0.3/(1.0+exp(-0.1*(v+32.0)));
    Xr_inft[i] = 1.0/(1.0+exp(-2.182-0.1819*v));
    Xr_taut[i] = 43.0+(1.0/(exp(-5.495+0.1691*v)+exp(-7.677-0.0128*v)));
    Xs_inft[i] = 1.0/(1.0+exp(-(v-16.0)/13.6));
    Xs_taut[i] = 1.0/((0.0000719*(v-10.0)/(1.0-exp(-0.148*(v-10.0))))+(0.000131*(v-10.0)/(exp(0.0687*(v-10.0))-1.0)));
    alpha_Xtot[i] = 0.04516*exp(0.03577*v);
    beta_Xtot[i] = 0.0989*exp(-0.06237*v);
    alpha_Ytot[i] = (0.005415*exp(-(v+33.5)/5.0))/(1.0+0.051335*exp(-(v+33.5)/5.0));
    beta_Ytot[i] = (0.005415*exp((v+33.5)/5.0))/(1.0+0.051335*exp((v+33.5)/5.0));
    d_inft[i] = 1.0/(1.0+exp(-(v-5.0)/6.24));
    f_inft[i] = 1.0/(1.0+exp((v+35.0)/8.6));
  }


  /* initial conditions */
  V = -96.08;
  cs[0] = 0.1;
  ci[0] = 0.1;
  cj[0] = 150.0;
  cjp[0] = 150.0;
  Irel[0] = 0.00001;
  d = 0.00000009255;
  f = 0.9992;
  q = 0.8775;
  m = 0.00024;
  h = 0.9989;
  j = 0.9964;
  Xr = 0.0000002286;
  Xs = 0.000264;
  Xto = 0.00003675;
  Yto = 1.0;

  i = 0;     

  E_Na = (R*Temp/F2)*log(Na_o/Na_i);
  E_K = (R*Temp/F2)*log(K_o/K_i);
  E_Ks = (R*Temp/F2)*log((K_o+0.01833*Na_o)/(K_i+0.01833*Na_i));


  while (t < t_end) {


    K1_inf = 1.0/(2.0+exp((1.62/(R*Temp/F2))*(V-E_K)));
    RV = 1.0/(1.0+2.5*exp(0.1*(V+28.0)));
    K_Kp = 1.0/(1.0+exp((7.488-V)/5.98));

    I_Na = g_Na*m*m*m*h*j*(V-E_Na);
    I_K1 = g_K1*K1_inf*(K_o/(K_o+Km_K1))*(V-E_K);
    I_Kr = g_Kr*RV*Xr*pow((K_o/4.0),0.5)*(V-E_K);
    I_Ks = g_Ks*Xs*Xs*(V-E_Ks);
    I_to = g_to*Xto*Yto*(V-E_K);
    I_Kp = g_Kp*K_Kp*(V-E_K);
    I_NaCa = 3.0/4.0*INaCa(cs[0],V);
    I_Ca = 2.7/3.5*Ica(d,f,q,cs[0],V);

    /* compute stimulus current */
    I_stim = 0.0;
    if (n_BCL == N_BCL) {
      stim_on = 1;
      n_BCL = 1;
    }
    if (n_stim <= N_stim && stim_on == 1) {
      I_stim = stim_amp;
      n_stim ++;
    }
    if (n_stim == N_stim) {
      stim_on = 0;
      n_stim = 1;
    }


    I_total = I_Na+I_K1+I_Kr+I_Ks+I_to+I_Kp+0.02*1000.0*(I_NaCa+2.0*I_Ca)+I_stim;


    cs[1] = cs[0]+dt*(beta_cs(cs[0])*vivs*(Irel[0]-((cs[0]-ci[0])/ts)-I_Ca+I_NaCa));
    ci[1] = ci[0]+dt*(beta_ci(ci[0])*(((cs[0]-ci[0])/ts)-Iup(ci[0])));
    cj[1] = cj[0]+dt*(-Irel[0]+Iup(ci[0]));
    cjp[1] = cjp[0]+dt*((cj[0]-cjp[0])/ta);
    Irel[1] = Irel[0]+dt*((-g*I_Ca*Q(cjp[0])/Ica_bar)-(Irel[0]/tr));

    
    /* extrapolate rate constants from tables */
    ilow = abs((V-V_min)/0.1);
    linext = (-V+vt[ilow])/0.1;
    d_inf = (d_inft[ilow+1]-d_inft[ilow])*linext+d_inft[ilow+1];
    f_inf = (f_inft[ilow+1]-f_inft[ilow])*linext+f_inft[ilow+1];
    alpha_m = (alpha_mt[ilow+1]-alpha_mt[ilow])*linext+alpha_mt[ilow+1];
    beta_m = (beta_mt[ilow+1]-beta_mt[ilow])*linext+beta_mt[ilow+1];
    tau_m = 1./(alpha_m+beta_m);
    alpha_h = (alpha_ht[ilow+1]-alpha_ht[ilow])*linext+alpha_ht[ilow+1];
    beta_h = (beta_ht[ilow+1]-beta_ht[ilow])*linext+beta_ht[ilow+1];
    tau_h = 1./(alpha_h+beta_h);
    alpha_j = (alpha_jt[ilow+1]-alpha_jt[ilow])*linext+alpha_jt[ilow+1];
    beta_j = (beta_jt[ilow+1]-beta_jt[ilow])*linext+beta_jt[ilow+1];
    tau_j = 1./(alpha_j+beta_j);
    Xr_inf = (Xr_inft[ilow+1]-Xr_inft[ilow])*linext+Xr_inft[ilow+1];
    Xr_tau = (Xr_taut[ilow+1]-Xr_taut[ilow])*linext+Xr_taut[ilow+1];
    Xs_inf = (Xs_inft[ilow+1]-Xs_inft[ilow])*linext+Xs_inft[ilow+1];
    Xs_tau = (Xs_taut[ilow+1]-Xs_taut[ilow])*linext+Xs_taut[ilow+1];
    alpha_Xto = (alpha_Xtot[ilow+1]-alpha_Xtot[ilow])*linext+alpha_Xtot[ilow+1];
    beta_Xto = (beta_Xtot[ilow+1]-beta_Xtot[ilow])*linext+beta_Xtot[ilow+1];
    tau_Xto = 1./(alpha_Xto+beta_Xto);
    alpha_Yto = (alpha_Ytot[ilow+1]-alpha_Ytot[ilow])*linext+alpha_Ytot[ilow+1];
    beta_Yto = (beta_Ytot[ilow+1]-beta_Ytot[ilow])*linext+beta_Ytot[ilow+1];
    tau_Yto = 1./(alpha_Yto+beta_Yto);


    d = d_inf+(d-d_inf)*exp(-dt/td);
    f = f_inf+(f-f_inf)*exp(-dt/tf);
    q = q+dt*(qinf(cs[0])-q)/tq;
    m = alpha_m*tau_m+(m-alpha_m*tau_m)*exp(-dt/tau_m);
    h = alpha_h*tau_h+(h-alpha_h*tau_h)*exp(-dt/tau_h);
    j = alpha_j*tau_j+(j-alpha_j*tau_j)*exp(-dt/tau_j/5.);  //j slowed down
    Xr = Xr_inf+(Xr-Xr_inf)*exp(-dt/Xr_tau);
    Xs = Xs_inf+(Xs-Xs_inf)*exp(-dt/Xs_tau);
    Xto = alpha_Xto*tau_Xto+(Xto-alpha_Xto*tau_Xto)*exp(-dt/tau_Xto);
    Yto = alpha_Yto*tau_Yto+(Yto-alpha_Yto*tau_Yto)*exp(-dt/tau_Yto);

    V = V-dt*I_total;

    cs[0] = cs[1];
    ci[0] = ci[1];
    cj[0] = cj[1];
    cjp[0] = cjp[1];
    Irel[0] = Irel[1];


    /* write to file */
    if (i==100 && t > t_end-4.*T)  {
     fprintf(output,"%.5f\t%.5f\t%.5f\n",t,ci[0],V);
      i=0;
    }


    t = t+dt;
    n_BCL++;
    if (t > t_end-4.*T) i++;

	   
  }

  fclose(output);
  return 0;


}