%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Beeler and Reuter's model defined in Keener's Mathematical Physiology   %
%   PP. 149-151                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%  This is a Matlab 'odefile' file (BR.m). Below are the                    %
%  vectors are column vectors. The vector s = [ V m h j d f x Cai] contains %
%  V the membrane potential, m, h, j, d, f, x  gating variables and Cai     %
%  internal Calcium.                                                        %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function s_prime=BR(t,s)
 global K Na X SS nn t_old AH BH
     V    = s(1); 
     m    = s(2);
     h    = s(3);
     j    = s(4);
     d    = s(5);
     f    = s(6);
     x    = s(7);
     Cai  = s(8);
     gate = [m, h, j, d, f, x];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     pp=0;
if pp == 1

  C=[ 0,      pi,    1,-1, -10,   -47;...
     40,     -17.86, 0, 0,  pi,   -72;...
     0.126,  -4,     0, 0,  pi,   -77;...
     1.7,     pi,    0, 1, -12.2, -22.5;...
     0.055,  -4,     0, 1, -5,    -78;...
     0.3,     pi,    0, 1, -10,   -32;...
     0.095,  -100,   0, 1, -13.9,  5;...
     0.07,   -58.5,  0, 1,  20,   -44;...
     0.012,  -125,   0, 1,  6.67, -28;...
     0.0065, -50,    0, 1, -5,    -30;...
     0.0005,  12,    0, 1,  17.5, -50;...
     0.0013, -16.67, 0, 1, -25,   -20];

  a_m = (C(1,1)*exp((V-C(1,6))/C(1,2))   + C(1,3)*(V-C(1,6)))...
         /(1+C(1,4)*exp((V-C(1,6))/C(1,5)));

  b_m = (C(2,1)*exp((V-C(2,6))/C(2,2))   + C(2,3)*(V-C(2,6)))...
         /(1+C(2,4)*exp((V-C(2,6))/C(2,5)));

  a_h = (C(3,1)*exp((V-C(3,6))/C(3,2))   + C(3,3)*(V-C(3,6)))...
         /(1+C(3,4)*exp((V-C(3,6))/C(3,5)));

  b_h = (C(4,1)                          + C(4,3)*(V-C(4,6)))...
         /(1+C(4,4)*exp((V-C(4,6))/C(4,5)));

  a_j = (C(5,1)*exp((V-C(5,6))/C(5,2))   + C(5,3)*(V-C(5,6)))...
           /(1+C(5,4)*exp((V-C(5,6))/C(5,5)));
 b_j = (C(6,1)                          + C(6,3)*(V-C(6,6)))...    
        /(1+C(6,4)*exp((V-C(6,6))/C(6,5)));
 a_d = (C(7,1)*exp((V-C(7,6))/C(7,2))   + C(7,3)*(V-C(7,6)))...
         /(1+C(7,4)*exp((V-C(7,6))/C(7,5)));
 b_d = (C(8,1)*exp((V-C(8,6))/C(8,2))   + C(8,3)*(V-C(8,6)))...
         /(1+C(8,4)*exp((V-C(8,6))/C(8,5)));
 a_f = (C(9,1)*exp((V-C(9,6))/C(9,2))   + C(9,3)*(V-C(9,6)))...
		            /(1+C(9,4)*exp((V-C(9,6))/C(9,5)));
 b_f = (C(10,1)*exp((V-C(10,6))/C(10,2)) + C(10,3)*(V-C(10,6)))...
         /(1+C(10,4)*exp((V-C(10,6))/C(10,5)));
  a_x = (C(11,1)*exp((V-C(11,6))/C(11,2)) + C(11,3)*(V-C(11,6)))...  
       /(1+C(11,4)*exp((V-C(11,6))/C(11,5)));
 b_x = (C(12,1)*exp((V-C(12,6))/C(12,2)) + C(12,3)*(V-C(12,6)))...
         /(1+C(12,4)*exp((V-C(12,6))/C(12,5)));
end
  AM=-(V+47)/(exp(-.1*(V+47))-1);
  BM=40*exp(-.056*(V+72));
  AH=.5*.126*exp(-.25*(V+77));  %should be 77
  BH=.5*1.7/(exp(-.082*(V+22.5))+1);
  AJ=.055*exp(-.25*(V+78))/(exp(-.2*(V+78))+1);%should be 78, not 58
  BJ=.3./(exp(-.1*(V+32))+1);
  AD=.095*exp(-.01*(V-5))./(exp(-.072*(V-5))+1);
  BD=.07*exp(-.017*(V+44))./(exp(.05*(V+44))+1);
  AF=.012*exp(-.008*(V+28))./(exp(.15*(V+28))+1);
  BF=.0065*exp(-.02*(V+30))./(exp(-.2*(V+30))+1);
  AX=.0005*exp(.083*(V+50))./(exp(.057*(V+50))+1);
  BX=.0013*exp(-.06*(V+20))./(exp(-.04*(V+20))+1);
  % for y = m, h, j, d, f, x 
  %  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];
  A = [AM, AH, AJ, AD, AF, AX];
  B = [BM, BH, BJ, BD, BF, BX];
  Inf = A./(A+B);                        %    y_inf = a_y/(a_y + b_y);
  m = Inf(1);
  Tau = 1./(A+B);                        
  %    t_y  = 1/(a_y + b_y);
  g_Na = 4*m^3*h*j+0.003;                % Na+ conductance
  %g_katp = 0.01; % with g_katp ne 0, system can be driven at much shorter BCL
  g_katp=.0;
  Vkatp = -70;
  % The reversal potentials (mv).
  v_Na =  50;                            % sodium equilibrium potential
  % Currents
  INa = g_Na*(V-v_Na);
  IK=1.4*(exp(0.04*(V+85))-1)/(exp(0.08*(V+53))+exp(0.04*(V+53))) ...
  +0.07*(V+23)/(1-exp(-0.04*(V+23)));
  IK = IK + g_katp*(V-Vkatp);
  Ix  = 0.8*x*(exp(0.04*(V+77))-1)/exp(0.04*(V+35));
  Is  = 0.09*d*f*(V+82.3+13.0287*log(Cai));
  if t==t_old | floor(t)==nn+1 
   nn=nn+~(t==t_old);
   Na(nn) = INa;
    K(nn)  = IK;
	X(nn)  = Ix;
	SS(nn)  = Is;
	t_old=t;
	end;
	% Membrane capacitance (micro-Farads/cm^2)
	C_m = 1; % HH
	% Exterior stimulus
	%  I = 3*heaviside(t-100)*heaviside(105-t); 
	%  I = 10*(atan((t-50)*10)/pi+1/2)*(atan((55-t)*10)/pi+1/2); 
	BCL = 200;
	if (mod(t,BCL)>50&mod(t,BCL)<55)
	I = -110;
	else
		I=0.;
	end

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%
	%
	% The derivatives of V, m, h, j, d, f and x with Cai as well.  
	%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	V_prime  = -(INa + IK + Ix + Is - I); 
	%  gate_prime = (Inf - gate)./Tau;   
	gate_prime = A.*(1-gate)-B.*gate;
	 Cai_prime = 0.07*(10^(-7) - Cai) - 10^(-7)*Is;
	 s_prime = [V_prime gate_prime Cai_prime]';
	 
