%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Two coupled cells each having ionics defined by 
%   Luo and Rudy's I model defined in Cir Res. Vol 68 No.6 June 1991        
%
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%                                                                           
%
%  This is a Matlab 'odefile' file (LR1_coupODE.m). Below are the           
%
%  vectors are column vectors. The vector                                   %
%
%  s = [V1 m1 h1 j1 d1 f1 x1 Cai1 V2 m2 h2 j2 d2 f2 x2 Cai2] contains 
%  V1,2 the membrane potentials,
%  m1,2 Na+ activating gate, 
%  h1,2 Na+ fast inactivating gate, 
%  j1,2 Na+ slow inactivating gate,
%  d1,2 Ca2+ activating gate, 
%  f1,2 Ca2+ inactivating gate,
%  x1,2 K+ gate, and 
%  Cai1,2 intracellular calcium variables.                                  %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function s_prime=LR1_coupODE(t,s)
 global iK iK1 iNa iSi ib iKp KKo1 Tm Th Tj Td Tf Tx step t_old nn t_end
     V(1)    = s(1);
     m(1)    = s(2);
     h(1)    = s(3);
	 j(1)    = s(4);
     d(1)    = s(5);
     f(1)    = s(6);
     x(1)    = s(7);
     Cai(1)  = s(8);
     V(2)    = s(9);
     m(2)    = s(10);
     h(2)    = s(11);
     j(2)    = s(12);
     d(2)    = s(13);
     f(2)    = s(14);
     x(2)    = s(15);
     Cai(2)  = s(16);
     gate = [m, h, j, d, f, x];
     w    = 1/2*(h+d);
     coup = 0.1; %coupling coefficient (units of 1/mohms)
	 coup = 0.;
	
     m = 0.9;     %relative mass of two regions
	 % Concentrations
	 % Ko = 5.4; %30.4 for oscillations
	 Ko = [5.4,5.4];
	 g_katp = [0.,0.];
	 %if t > 200
	 %  Ko = 5.4 + 144.6*t/1000;
     %  Ko = 5.4 + 144.6*0.0648;
     %  Ko = 5.4 + 144.6*(t-200)/(t_end-200)
	 %end
	 Ki = 145;
     Nao = 140;
     Nai = 18;
	 PRNaK = 0.01833;  % Na - K exchanger rate
	 R = 8.315; %(J/(mol*K)
	 T = 310;   %(K)  
	 F = 96.49; %(kC/mol where kC = J/mV);
	 % The reversal potentials (mv).
	 v_Na =  54.4;                                          % Na+ eq. pot
	 v_Si =  7.77-13.0287*log(Cai);
	 % v_K's are 2-vectors
	 v_K  =  R*T/F*log((Ko + PRNaK*Nao)/(Ki + PRNaK*Nai));  % K+ eq pot w/ xchg
	 v_K1 =  R*T/F*log(Ko/Ki);                              % K+ eq pot
	 v_Kp =  v_K1;                                          % K+ eq pot
	 v_b  = -59.87;                                         % bg K+ eq pot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The steady state probability functions m_inf, h_inf, and n_inf 
% (similar to the expressions of the Hodgkin-Huxley equations), and
% the relaxation time t_n (sec).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if V(1)<-40 
     a_h(1) = 0.135*exp((80+V(1))/(-6.8));
	 a_j(1) = (-1.2714*10^5*exp(0.2444*V(1))-...
	       3.474*10^-5*exp(-0.04391*V(1)))*(V(1)+37.78)/(1+exp(0.311*(V(1)+79.23)));
	 b_h(1) = 3.56*exp(0.079*V(1))+3.1*10^5*exp(0.35*V(1));
	 b_j(1) = 0.1212*exp(-0.01052*V(1))/(1+exp(-0.1378*(V(1)+40.14)));
	 else
	 a_h(1) = 0;
	 a_j(1) = 0;
	 b_h(1) = 1/(0.13*(1+exp((V(1)+10.66)/(-11.1))));
	 b_j(1) = 0.3*exp(-2.535*10^-7*V(1))/(1+exp(-0.1*(V(1)+32)));
	 end
if V(2)<-40 
	 a_h(2) = 0.135*exp((80+V(2))/(-6.8));
	 a_j(2) = (-1.2714*10^5*exp(0.2444*V(2))-...
	       3.474*10^-5*exp(-0.04391*V(2)))*(V(2)+37.78)/(1+exp(0.311*(V(2)+79.23)));
     b_h(2) = 3.56*exp(0.079*V(2))+3.1*10^5*exp(0.35*V(2));
	 b_j(2) = 0.1212*exp(-0.01052*V(2))/(1+exp(-0.1378*(V(2)+40.14)));
	 else
	 a_h(2) = 0;
	 a_j(2) = 0;
     b_h(2) = 1/(0.13*(1+exp((V(2)+10.66)/(-11.1))));
	 b_j(2) = 0.3*exp(-2.535*10^-7*V(2))/(1+exp(-0.1*(V(2)+32)));
	 end
	 a_m = 0.32*(V+47.13)./(1 - exp(-0.1*(V+47.13)));
	 b_m = 0.08*exp(-V/11);
	 a_d = 0.095*exp(-0.01*(V-5))./(1+exp(-0.072*(V-5)));
	 b_d = 0.07*exp(-0.017*(V+44))./(1+exp(0.05*(V+44)));
	 a_f = 0.012*exp(-0.008*(V+28))./(1+exp(0.15*(V+28)));
	 b_f = 0.0065*exp(-0.02*(V+30))./(1+exp(-0.2*(V+30)));
if V(1) > - 100     
	xi(1) = 2.837*(exp(0.04*(V(1)+77))-1)./((V(1)+77)*exp(0.04*(V(1)+35)));
	else
	xi(1) = 1;
	end
if V(2) > - 100     
	xi(2) = 2.837*(exp(0.04*(V(2)+77))-1)./((V(2)+77)*exp(0.04*(V(2)+35)));
	else
    xi(2) = 1;
	end
a_x = 0.0005*exp(0.083*(V+50))./(1+exp(0.057*(V+50)));
b_x = 0.0013*exp(-0.06*(V+20))./(1+exp(-0.04*(V+20)));
a_k1 = 1.02./(1+exp(0.2385*(V-v_K1-59.215)));
b_k1 = (0.49124*exp(0.08032*(V-v_K1+5.476))+...
	              exp(0.06175*(V-v_K1-594.31)))./(1+exp(-0.5143*(V-v_K1+4.753)));
Kp = 1./(1+exp((7.488-V)/5.98));
A = [a_m, a_h, a_j, a_d, a_f, a_x];
B = [b_m, b_h, b_j, b_d, b_f, b_x];
Inf = A./(A+B);        
Tau = 1./(A+B);        
	%   g_Na = 23*Inf(1)^3*h*h;
g_Na = 23*m.^3.*h.*j;                  % Cell1 Na+ conductance
%   g_Si = 0;
%   g_Si = 0.09*d*(1-1.3*x)/0.8;              % Ca2+ conductance
%   g_Si = 0.09*(1-h)*f;                      % Ca2+ conductance
g_Si = 0.09*d.*f;                           % Ca2+ conductance
% g_K and g_K1 are 2-vectors
g_K  = 0.282*(Ko./5.4).^(1/2).*x.*xi;       % K+ conductance (t-dep)
g_K1 = 0.6047*(Ko./5.4).^(1/2).*(a_k1./(a_k1+b_k1));%K+ conductance (t-indep)
g_Kp = 0.0183*Kp;                                % K+ conductance (t-indep)
g_b  = 0.03921;                                  % K+ conductance (bg)
% Membrane capacitance (micro-Farads?)
C_m = 1; % HH
% Exterior stimulus
%    I = 0;
%    I = 10;
%    I = 10*(heaviside(t-1200)*heaviside(1300-t)+...
%            heaviside(t-1000)*heaviside(1100-t)+...
%            heaviside(t-1)*heaviside(101-t)); 
%    I = Iapp(t);
%    I = 10*(atan((t-100)*10)/pi+1/2)*(atan((120-t)*10)/pi+1/2); 
	BCL = 650;
	if (mod(t,BCL)>50&mod(t,BCL)<55)
	I(1) = 10;
	I(2) = 6;
	else
		I=0.;
	end
	
	
	%if (t > 100 & t < 105)|(t > 650 & t < 655)
   % I(1) = 0; %10 for oscillations
	%else
	% I(1) = 0;
	%end
	% I(2) = 0; %10 for oscillations
% Currents
INa = g_Na.*(V-v_Na); 
ISi = g_Si.*(V-v_Si);
IK  = g_K.*(V-v_K);
IK1 = g_K1.*(V-v_K1);
IKp = g_Kp.*(V-v_Kp);
Ib  = g_b.*(V-v_b);
I_Katp = g_katp.*(V-v_K1);
%pause;
%t
if t==t_old | floor(t)==nn+1 
    nn=nn+~(t==t_old);  
%(t==t_old) = 1 if true; =0 if not true
%if (t==t_old) = 1, then ~(t==t_old) = 0
%if (t==t_old) = 0, then ~(t==t_old) = 1
iNa(nn) = INa(1);
iK(nn)  = IK(1);
iK1(nn) = IK1(1);  
iKp(nn) = IKp(1);
ib(nn)  = Ib(1);
iSi(nn) = ISi(1);
KKo1(nn)  = Ko(1);
%    KKo2(nn)  = Ko(2);
Tm(nn) = Tau(1,1);
Th(nn) = Tau(1,2);
Tj(nn) = Tau(1,3);
Td(nn) = Tau(1,4);
Tf(nn) = Tau(1,5);
Tx(nn) = Tau(1,6);
t_old=t;
  end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The derivatives of V, m, h, j, d, f, x, and Cai.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
V_prime=-1/C_m*(INa+ISi+IK+IK1+I_Katp+IKp+Ib)+I+...
          [-1/m*coup*(V(1)-V(2)),1/(1-m)*coup*(V(1)-V(2))];
 gate_prime = (Inf - gate)./Tau;
 %  Cai_prime  = [0,0]; 
 Cai_prime  = -10^(-4)*ISi+0.07*(10^(-4)-Cai);
 s_prime = [V_prime(1,1), gate_prime(1,1:2:11), Cai_prime(1,1),...
             V_prime(1,2), gate_prime(1,2:2:12), Cai_prime(1,2)]'; 
			 
