% A_current.m

% This is the Hodgkin-Huxley model with the addition of an A current. From
% page 107 of the text.

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

%  initial conditions
   V0=-60;
   m0=0;
   h0=1;
   n0=0;
   a0=0.5;
   b0=0.2;
   
%  initial condition vector
   y0 = [V0 m0 h0 n0 a0 b0];

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

%  Applied current parameter (muA/cm^2)
%   Is=6.53; % this gives a delay before a single spike
    Is=8; % this gives a train of spikes
   gAbar = 47.7; 
%  gAbar = 0; 
   VL=-22;  % value with A-current
%  VL=-73;   % value without A-current
    
%  Call ODE solver, pass some parameters
   [t,y] = ode45(@(t,y)RHS(t,y,ton,toff,Is,gAbar,VL),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);
   a=y(:,5);
   b=y(:,6);
   
%  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

%  A current conductance for plotting
   gA = gAbar.*a.^3.*b;

   figure(1);
%  plot A current conductance time course
   ha=plot(t,gA,'black');
   set(ha,'LineWidth',2);
   axis([0 tmax 0 10]);
   box off;
   xlabel('t (ms)');
   ylabel('g_A');
   
   figure(2);
%  plot V time course
   ha=plot(t,V,'black');
   set(ha,'LineWidth',2);
   axis([0 tmax -80 30]);
   box off;
   xlabel('t (ms)');
   ylabel('V (mV)');

   figure(3);
   % plot AHP activation variable
   ha=plot(t,a,'black','Linewidth',2);
   axis([0 tmax 0 1]);
   box off;
   xlabel('t (ms)');
   ylabel('a');
   
   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,gAbar,VL);
 
%     set up the variable vectors
      V=y(1);
      m=y(2);
      h=y(3);
      n=y(4);
      a=y(5);
      b=y(6);

%     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
      lambda_a=1;  % scaling factor for "a" time constant
      lambda_b=1;  % scaling factor for "b" time constant
      ainf = ((0.0761*exp((V+99.22)/31.84))/(1+exp((V+6.17)/28.93)))^(1/3);
      taua = lambda_a*(0.3632+1.158/(1+exp((V+60.96)/20.12)));
      binf = (1/(1+exp((V+58.3)/14.54)))^4;
      taub = lambda_b*(1.24+2.678/(1+exp((V-55)/16.027)));

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

%  calculate conductances
      gNa = gNabar*(m^3)*h;
      gDR = gDRbar*n^4;
      gA = gAbar*a^3*b;

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

%     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-IA+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;
      adot = (ainf-a)/taua;
      bdot = (binf-b)/taub;
      dydt=[Vdot;mdot;hdot;ndot;adot;bdot];

   end
end 
