%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Luo-Rudy 1 model (Luo & Rudy, Circ Res 68:1501-1526, 1991)
%  Mathematical model of guinea pig cardiac action potential
%
%  T. Krogh-Madsen, May 2010
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function p=LR1_alternans(w,C)

% PARAMETERS
% maximal conductances
g_Na = 16.0;
g_Ca = 0.06;
g_K = 0.432;
g_K1 = 0.6047;
g_Kp = 0.0183;

% other constants
E_b = -59.87;
RT_on_F = 26.7; 
Ki = 145;
Ko = 5.4; 
Nai = 18; 
Nao = 140;  
Cm = 1.0;

% simulation parameters
stim_amp = -100;
dt = 0.01;   
t_end = 770;
counter = 0;
T = 195;

% initial conditions
t = 0;
Ca = 1.6806e-04;
V = -83.4855;
m = 0.0020;
h = 0.9788;
j = 0.9733;
d = 0.0033;
f = 0.9638;
X = 0.1779;
APDstart = 0;
APDcounter = 1;

% compute reversal potentials
E_Na = RT_on_F*log(Nao/Nai);
E_K = RT_on_F*log((Ko+0.01833*Nao)/(Ki+0.01833*Nai));
E_K1 = RT_on_F*log(Ko/Ki);
E_Kp = E_K1;

% define vectors for storing results
v = zeros(t_end/dt/10,1);
INa = zeros(t_end/dt/10,1);
ICa = zeros(t_end/dt/10,1);
IK = zeros(t_end/dt/10,1);
IK1 = zeros(t_end/dt/10,1);
IKp = zeros(t_end/dt/10,1);
Ib = zeros(t_end/dt/10,1);
Cai = zeros(t_end/dt/10,1);

% integrate model
while t < t_end

  % m-gate (sodium activation)   
  alpha_m = 0.32*(V+47.13)/(1.0-exp(-0.1*(V+47.13)));
  beta_m = 0.08*exp(-V/11.0);
  tau_m = 1./(alpha_m+beta_m);
  m_inf = alpha_m*tau_m;

  % h-gate (sodium inactivation)   
  if (V >= -40.0)
    alpha_h = 0.0;
	beta_h = 1.0/(0.13*(1.0+exp((V+10.66)/(-11.1))));
  else
	alpha_h = 0.135*exp((80.0+V)/(-6.8));
	beta_h = 3.56*exp(0.079*V)+3.1e5*exp(0.35*V);
  end
  tau_h = 1./(alpha_h+beta_h);
  h_inf = alpha_h*tau_h;

  % j-gate (sodium inactivation)   
  if (V >= -40.0)
	alpha_j = 0.0;
    beta_j = 0.3*exp(-2.535e-7*V)/(1.0+exp(-0.1*(V+32.0)));
  else
	alpha_j = (-1.2714e5*exp(0.2444*V)-3.474e-5*exp(-0.04391*V))*(V+37.78)/(1.0+exp(0.311*(V+79.23)));
	beta_j = 0.1212*exp(-0.01052*V)/(1.0+exp(-0.1378*(V+40.14)));
  end
  tau_j = 1./(alpha_j+beta_j);
  j_inf = alpha_j*tau_j;
  
  % d-gate (calcium activation)   
  alpha_d = 0.095*exp(-0.01*(V-5.0))/(1.0+exp(-0.072*(V-5.0)));
  beta_d = 0.07*exp(-0.017*(V+44.0))/(1.0+exp(0.05*(V+44.0)));
  tau_d = 1./(alpha_d+beta_d);
  d_inf = alpha_d*tau_d;

  % f-gate (calcium inactivation)   
  alpha_f = 0.012*exp(-0.008*(V+28.0))/(1.0+exp(0.15*(V+28.0)));
  beta_f = 0.0065*exp(-0.02*(V+30.0))/(1.0+exp(-0.2*(V+30.0)));
  tau_f = 1./(alpha_f+beta_f);
  f_inf = alpha_f*tau_f;

  % X-gate (potassium activation)   
  alpha_X = 0.0005*exp(0.083*(V+50.0))/(1.0+exp(0.057*(V+50.0)));
  beta_X = 0.0013*exp(-0.06*(V+20.0))/(1.0+exp(-0.04*(V+20.0)));
  tau_X = 1./(alpha_X+beta_X);
  X_inf = alpha_X*tau_X;

  % compute gating variables
  m = m_inf+(m-m_inf)*exp(-dt/tau_m);
  h = h_inf+(h-h_inf)*exp(-dt/tau_h);
  j = j_inf+(j-j_inf)*exp(-dt/tau_j);
  d = d_inf+(d-d_inf)*exp(-dt/tau_d);
  f = f_inf+(f-f_inf)*exp(-dt/tau_f);
  X = X_inf+(X-X_inf)*exp(-dt/tau_X);

  % integrate intracellular calcium 
  dCa = -1.0e-4*(g_Ca*d*f*(V-(7.7-13.0287*log(Ca))))+0.07*(1.0e-4-Ca);
  Ca = Ca+dt*dCa;

  % time-independent gates
  K1_inf = (1.02/(1.0+exp(0.2385*(V-E_K1-59.215))))/((1.02/(1.0+exp(0.2385*(V-E_K1-59.215))))+((0.49124*exp(0.08032*(V-E_K1+5.476))+exp(0.06175*(V-E_K1-594.31)))/(1.0+exp(-0.5143*(V-E_K1+4.753)))));
  if (V > -100.0)
	Xi = 2.837*(exp(0.04*(V+77.0))-1.0)/((V+77.0)*exp(0.04*(V+35.0)));
  else
	Xi = 1.0;
  end
  
  % membrane currents
  I_Na = g_Na*m*m*m*h*j*(V-E_Na);
  I_Ca = g_Ca*d*f*(V-(7.7-13.0287*log(Ca)));
  I_K = g_K*sqrt(Ko/5.4)*X*Xi*(V-E_K);
  I_K1 = g_K1*sqrt(Ko/5.4)*K1_inf*(V-E_K1);
  I_Kp = g_Kp*(V-E_Kp)/(1.0+exp((7.488-V)/5.98));
  I_b = 0.03921*(V-E_b);

  % stimulus current  
  if mod(t,T)<=0.5,
	I_stim = stim_amp;
  else
    I_stim = 0.0;
  end

  % total current      
  I_total = I_Na + I_Ca + I_K + I_K1 + I_Kp + I_b + I_stim;

  % integrate voltage      
  V = V - dt*I_total/Cm;

  % update time 
  counter = counter+1;
  t = dt*counter;
  
  % compute action potential duration (APD)   
  if (V>-75. && V+dt*I_total/Cm<=-75.)
      APDstart=t;
  end
  if (V<-75. && V+dt*I_total/Cm>=-75.)
      APD(APDcounter,1)=t;
      APD(APDcounter,2)=t-APDstart;     
      APDcounter = APDcounter+1;
  end
  
  
  % save results (downsampled) in vectors
  if (mod(counter,10)==0)
    v(counter/10)=V;
    INa(counter/10)=I_Na;
    ICa(counter/10)=I_Ca;
    IK(counter/10)=I_K;
    IK1(counter/10)=I_K1;
    IKp(counter/10)=I_Kp;
    Ib(counter/10)=I_b;
    Cai(counter/10)=Ca;
  end
  
end

t = dt:dt*10:t_end;
% plot 
figure(w), hold on
subplot(221), hold on, plot(t,v,'Color',C), ylabel('V (mV)')
subplot(222), hold on, plot(APD(:,1),APD(:,2),'*-','Color',C), ylabel('APD (ms)')
subplot(223), hold on, plot(t,ICa,'Color',C), ylabel('I_{Ca} (pA/pF)'), xlabel('Time (ms)')
subplot(224), hold on, plot(t,IK,'Color',C), ylabel('I_{K} (pA/pF)'), xlabel('Time (ms)')

