% ===============================================================
% == Single Cell Minimal Model (4 state variables)             ==
% ==                                                           ==
% ==   Supplement for the CMACS Workshop                       ==
% ==                                                           ==
% ==                                                           ==
% ==   Author:                                                 ==
% ==                                                           ==
% ==     E. Bartocci                                           ==
% ==                                                           == 
% ==   Date:  11/05/10                                         ==
% ==                                                           ==
% ==   Free distribution with authors permission               ==
% ==                                                           ==
% ==   SUNY Stony Brook, Stony Brook, NY                       ==
% ==                                                           == 
% ===============================================================  

% The following are the parameters that you can find in the paper 
% A. Bueno-Orovio, M. Cherry, and F. Fenton, ?Minimal model for 
% human ventricular action potentials in tissue,? Journal of 
% Theoretical Biology, no. 253, pp. 544?560, 2008. 

% Questions:
%
% What is the maximum dt such that I can still observe the solution ?
% Modify the parameters in order to increase the action potential duration.
% Modify the parameters in order to obtain alternans
% Find the range of the parameters such that the cell is still excitable
% Modify the parameters in order to increase the length of the peak of the upstroke.
% Eliminate divisions from the Minimal Resistor using the Sigmoid Closure



function [out] = single_cell_mm4V ()

         global EPI_TVP; 
         global EPI_TV1M;  
         global EPI_TV2M; 

         global EPI_TWP;   

         global EPI_TW1P; %200
         global EPI_TW2P; %200

         global EPI_TW1M; %190    
         global EPI_TW2M;
         
         global EPI_TS1;
         global EPI_TS2;
         global EPI_TFI;
         global EPI_TO1;
         global EPI_TO2;
         global EPI_TSO1;
         global EPI_TSO2;

         global EPI_TSI;
         %  TSI1     1.8875
         %  TSI2     1.8875


         global EPI_TWINF;
         global EPI_THV;
         global EPI_THVM;
         global EPI_THVINF;
         global EPI_THW;
         global EPI_THWINF;
         global EPI_THSO;
         global EPI_THSI;
         global EPI_THO;

         %  KWP      5.7
         global EPI_KWM;
         global EPI_KS;
         global EPI_KSO;  
         %  KSI      97.8
         global EPI_UWM;    
         global EPI_US;  
         global EPI_UO;     
         global EPI_UU;  
         global EPI_USO;   
         global EPI_SC;
         %  WCP      0.15

         global EPI_WINFSTAR;  

         global EPI_TW2M_TW1M;
         global EPI_TSO2_TSO1;

         global EPI_TW2M_TW1M_DIVBY2;
         global EPI_TSO2_TSO1_DIVBY2;

         EPI_TVP   =     1.4506; 
         EPI_TV1M  =    60.;  
         EPI_TV2M  =  1150.; 

         EPI_TWP   =   200.0;   

         EPI_TW1P  =    60.0; %200
         EPI_TW2P  =    15.0; %200

         EPI_TW1M  =    60.0; %190    
         EPI_TW2M  =    15.;  %190   %The same with Flavio's paper
         
         EPI_TS1   =    2.7342;  %The same with Flavio's paper
         EPI_TS2   =   16.;     %The same with Flavio's paper
         EPI_TFI   =    0.11;    %The same with Flavio's paper
         EPI_TO1   =  400.;    %The same with Flavio's paper
         EPI_TO2   =    6.  ;    %The same with Flavio's paper
         EPI_TSO1  =   30.0181; %The same with Flavio's paper
         EPI_TSO2  =    0.9957;  %The same with Flavio's paper

         EPI_TSI   =    1.8875;  % We have TSI1 and TSI2 = TSI in Flavio's paper
         %  TSI1     1.8875
         %  TSI2     1.8875


         EPI_TWINF  =  0.07;    %The same with Flavio's paper
         EPI_THV    =  0.3;     %EPUM % The same of Flavio's paper
         EPI_THVM   =  0.006;   %EPUQ % The same of Flavio's paper
         EPI_THVINF =  0.006;   %EPUQ % The same of Flavio's paper
         EPI_THW    =  0.13;    %EPUP % The same of Flavio's paper
         EPI_THWINF =  0.006;    %EPURR % In Flavio's paper 0.13
         EPI_THSO   =  0.13;    %EPUP % The same of Flavio's paper
         EPI_THSI   =  0.13;    %EPUP % The same of Flavio's paper
         EPI_THO    =  0.006;    %EPURR % The same of Flavio's paper

         %  KWP      5.7
         EPI_KWM    =  65.;     %The same of Flavio's paper
         EPI_KS     =  2.0994;  %The same of Flavio's paper
         EPI_KSO    =  2.0458;  %The same of Flavio's paper
         %  KSI      97.8
         EPI_UWM    =  0.03;    %The same of Flavio's paper
         EPI_US     =  0.9087;  %The same of Flavio's paper
         EPI_UO     =  0.;     % The same of Flavio's paper
         EPI_UU     =  1.55;   % The same of Flavio's paper
         EPI_USO    =  0.65;   % The same of Flavio's paper
         EPI_SC     =  0.007;
         %  WCP      0.15

         EPI_WINFSTAR = 0.94;   % The same of Flavio's paper

         EPI_TW2M_TW1M = -45.0;
         EPI_TSO2_TSO1 = -29.0224;

         EPI_TW2M_TW1M_DIVBY2 = -22.5;
         EPI_TSO2_TSO1_DIVBY2 = -14.5112;
         

         ut1     = [];
         vt1     = [];
         wt1     = [];
         st1     = [];
         stimt1  = [];
      
         u = 0.0;
         v = 1.0;
         w = 1.0;
         s = 0.0;
         t1 = 0;
         t2 = 0;

         for i=1:20000   
             [u,v,w,s, t1, t2, stim] = nextStepNonlinear (u,v,w,s, t1, t2, 0.1);
             ut1(i) =    u;
             vt1(i) =    v;
             wt1(i) =    w;
             st1(i) =    s;
             stimt1(i)=  stim;     
         end
         
         
         %hold off;
         plot(linspace(0,2000, 20000), [ut1; vt1; wt1; st1; stimt1]);
         
         %plot(linspace(0,30, 30000), ut3);
         
end


function [u,v,w,s, t1, t2, stim] = nextStepNonlinear (u,v,w,s, t1, t2, dt)
         global EPI_TVP; 
         global EPI_TV1M;  
         global EPI_TV2M; 

         global EPI_TWP;   

         global EPI_TW1P; %200
         global EPI_TW2P; %200

         global EPI_TW1M; %190    
         global EPI_TW2M;
         
         global EPI_TS1;
         global EPI_TS2;
         global EPI_TFI;
         global EPI_TO1;
         global EPI_TO2;
         global EPI_TSO1;
         global EPI_TSO2;

         global EPI_TSI;
         %  TSI1     1.8875
         %  TSI2     1.8875


         global EPI_TWINF;
         global EPI_THV;
         global EPI_THVM;
         global EPI_THVINF;
         global EPI_THW;
         global EPI_THWINF;
         global EPI_THSO;
         global EPI_THSI;
         global EPI_THO;

         %  KWP      5.7
         global EPI_KWM;
         global EPI_KS;
         global EPI_KSO;  
         %  KSI      97.8
         global EPI_UWM;    
         global EPI_US;  
         global EPI_UO;     
         global EPI_UU;  
         global EPI_USO;   
         global EPI_SC;
         %  WCP      0.15

         global EPI_WINFSTAR;  

         global EPI_TW2M_TW1M;
         global EPI_TSO2_TSO1;

         global EPI_TW2M_TW1M_DIVBY2;
         global EPI_TSO2_TSO1_DIVBY2;


         t1 = t1 + dt;
         t2 = t2 + dt;
         stim = HEAVISIDE(t1 - 0)*(1 - HEAVISIDE(t2 - 1)) + HEAVISIDE(t1 - 350)*(1 - HEAVISIDE(t2 - 351)) + HEAVISIDE(t1 - 700)*(1 - HEAVISIDE(t2 - 701)) + HEAVISIDE(t1 - 1100)*(1 - HEAVISIDE(t2 - 1101));
         

         
         jfi = 0.0;
         jso = 0.0;
         jsi = 0.0;
         
         if u < 0.006
              
              %w = w + (1.0 -(u/EPI_TWINF) - w) * (0.04166666666666667 + 0.025 * tanh(EPI_KWM*(u- 0.0406638027778453)))*dt;
              w = w + ((1.0 -(u/EPI_TWINF) - w)/(EPI_TW1M + (EPI_TW2M - EPI_TW1M) * 0.5 * (1.+tanh(EPI_KWM*(u-EPI_UWM)))))*dt;     
              v = v + ((1.0-v)/EPI_TV1M)*dt;
              s = s + ((((1.+tanh(EPI_KS*(u - EPI_US))) * 0.5) - s)/EPI_TS1)*dt;
              jfi = 0.0;
              jso = u/EPI_TO1;
              jsi = 0.0;  
              
         elseif u < 0.13
              %w = w + ((0.94-w) * (0.04166666666666667 + 0.025 * tanh(EPI_KWM*(u- 0.0406638027778453 ))))*dt;
              w = w + ((0.94-w)/(EPI_TW1M + (EPI_TW2M - EPI_TW1M) * 0.5 * (1.+tanh(EPI_KWM*(u-EPI_UWM)))))*dt;
              
              v = v + (-v/EPI_TV2M)*dt;
              s = s +((((1.+tanh(EPI_KS*(u-EPI_US))) * 0.5) - s)/EPI_TS1)*dt;
              jfi = 0.0;
              jso = u/EPI_TO2;
              jsi = 0.0;
         
         elseif u < 0.3
              w = w + (-w/EPI_TWP)*dt;
              v = v + (-v/EPI_TV2M)*dt;
              s = s + ((((1.+tanh(EPI_KS*(u-EPI_US))) * 0.5) - s)/EPI_TS2)*dt;
              jfi = 0.0;
              jso = 0.518815902099483 + 0.485502667750873 * tanh(EPI_KSO * (u -  1.48246402499052 ));
              jsi = -w * s/EPI_TSI;
         else
              w  = w + (-w/EPI_TWP)*dt;
              v  = v + (-v/EPI_TVP)*dt;
              s  = s +((((1.+tanh(EPI_KS*(u - EPI_US))) * 0.5) - s)/EPI_TS2)*dt;
                       
              jfi = -v * (u - EPI_THV) * (EPI_UU - u)/EPI_TFI;
              jso = 0.518815902099483 + 0.485502667750873 * tanh(EPI_KSO * (u -  1.48246402499052 ));
              jsi = -w * s/EPI_TSI;  
         end

         u = u  - (jfi+jso+jsi-stim)*dt;
                 
end


