% Nap_current.m

% This is the Hodgkin-Huxley model with the addition of a persistent Na current.% The HH code is from page 107 of the text. The Nap current is from the
% Ross et al. JNS 2019 paper.

function [t,y] = call_RHS();
   clear;
   tmax = 1000;
   tspan = [0,tmax];

%  initial conditions
   V0=-60;
   m0=0;
   h0=1;
   n0=0;
   e0=1;
   
%  initial condition vector
   y0 = [V0 m0 h0 n0 e0] ;

%  Current pulse timing
   ton=100;
   toff=900;

%  Applied current parameter (muA/cm^2)
   Is=6; % default is 6
   gNapbar=10; % default is 10
   lambda_e=1; % default is 1
    
%  Call ODE solver, pass some parameters
   [t,y] = ode45(@(t,y)RHS(t,y,ton,toff,Is,gNapbar,lambda_e),tspan,y0);

%  change defaults for plotting
   set(0,'DefaultAxesFontSize',16)

%  set up variable vectors
   V=y(:,1);
   m=y(:,2);
   h=y(:,3);
   n=y(:,4);
   e=y(:,5);
   
%  applied current pulse for plotting
   Nsizevec=size(n);
   Nsize=Nsizevec(1);
   tme=linspace(0,tmax,Nsize);
   for j=1:Nsize
       if tme(j) > ton & tme(j) < toff
           Iap(j)=Is;  % pulse
%           Iap(j)=((tme(j)-ton)/(toff-ton))*Is;  % ramp
       else
           Iap(j)=0;
       end
   end

%   Persistent Na current for plotting
    VNa = 50;    
    gNap=gNapbar.*(m.^3).*e;

   figure(1);
%  plot gNap time course
   ha=plot(t,gNap,'black');
   set(ha,'LineWidth',2);
   axis([0 tmax 0 10]);
   box off;
   xlabel('t (ms)');
   ylabel('g_Nap');
   
   figure(2);
%  plot V time course
   ha=plot(t,V,'black');
   set(ha,'LineWidth',2);
   axis([0 tmax -80 50]);
   box off;
   xlabel('t (ms)');
   ylabel('V (mV)');

   figure(3);
   % plot AHP activation variable
   ha=plot(t,e,'black','Linewidth',2);
   axis([0 tmax 0 1]);
   box off;
   xlabel('t (ms)');
   ylabel('e');
   
%   figure(4);
   % plot Ca concentration 
%   ha=plot(t,b,'black','Linewidth',2);
%   axis([0 tmax 0 1]);
%   box off;
%   xlabel('t (ms)');
%   ylabel('b');
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Define the RHS function for the ODEs
   function dydt = RHS(t,y,ton,toff,Is,gNapbar,lambda_e);
 
%     set up the variable vectors
      V=y(1);
      m=y(2);
      h=y(3);
      n=y(4);
      e=y(5);

%     Rate coefficients
      alpha_m = -3.8*(0.1*(V+34.7)/(exp(-(V+34.7)/10)-1));
      beta_m = 3.8*4*exp(-(V+59.7)/18);
      alpha_h = 3.8*0.07*exp(-(V+53)/20);
      beta_h = 3.8/(exp(-(V+23)/10)+1);
      alpha_n = -(3.8/2)*(0.01*(V+50.7)/(exp(-(V+50.7)/10)-1));
      beta_n = (3.8/2)*0.125*exp(-(V+60.7)/80);

%  infinity functions and time constants for A current
      einf = 1/(1+exp((V+48)/6));
      taue = lambda_e*1000/cosh((V+48)/12);

%  parameters
      Cm = 1;
      gNabar = 120; gLbar = 0.3; gDRbar = 20; 
      
      VNa = 50; VK = -77; 
      VL=-70;   

%  calculate conductances
      gNa = gNabar*(m^3)*h;
      gDR = gDRbar*n^4;
      gNap = gNapbar*(m^3)*e;

%  calculate currents
      IL=gLbar*(V-VL);
      INa=gNa*(V-VNa);
      IDR=gDR*(V-VK);
      INap=gNap*(V-VNa);

%     Current applied during ramp
      if t > ton & t < toff
         Iin = Is;   % pulse
%         Iin = ((t-ton)/(toff-ton))*Is;   % ramp
      else
         Iin = 0;
      end

%     Differential equations
      Vdot = (-IL-INa-IDR-INap+Iin)/Cm;  
      mdot = alpha_m*(1-m)-beta_m*m;
      hdot = alpha_h*(1-h)-beta_h*h;
      ndot = alpha_n*(1-n)-beta_n*n;
      edot = (einf-e)/taue;
      dydt=[Vdot;mdot;hdot;ndot;edot];

   end
end 
