
%
% Cox & Gabbiani, Mathematics for Neuroscientists
%
% Staggered Euler spike frequency
%
% usage:  stEfreq(IOlo:I0step:I0hi)
% 
% where:    IOlo:I0step:I0hi is current voltage range and step
% 
% where:    I0lo = beginning current
%           I0step = increase in I0 per time
%           I0hi = ending current
%
% example:  stEfreq(30:2:300)
%

function stEfreq(I0)

stim = struct('t1',1,'t2',2,'amp',4e-4,'loc',0.05,'Tfin',100);

for j=1:length(I0),
    
    stim.amp = I0(j);
    freq(j) = 1000/stE_isi(.01,stim);    % change (ms) to (s)

end

plot(I0,freq,'k')
xlabel('I_0   (pA)','fontsize',14)
ylabel('spikes per second','fontsize',14)
box off

%
% stE_isi.m
%
% Staggered Euler on the active sphere, returns inter spike interval
%
% usage:  isi = stE_isi(dt,stim)
% 

function isi = stE_isi(dt,stim)

Nt = ceil(stim.Tfin/dt);                         % number of time steps
A = 4*pi*(1e-6) ;                           % (cm)^2
E = struct('K', -77, 'Na', 56, 'Cl', -68);  % reversal potentials, mV
G = struct('K', 36, 'Na', 120, 'Cl', 0.3);  % channel conductances, mS/cm^2
Vr = fsolve(@(V) Iss(V,E,G),-71);           % find rest potential 

t = zeros(Nt,1); V = t;      % preallocate space               
t(1) = 0;
V(1) = Vr;
n = an(Vr)/(an(Vr)+bn(Vr));  % gating variable
m = am(Vr)/(am(Vr)+bm(Vr));  % activation voltage (up)
h = ah(Vr)/(ah(Vr)+bh(Vr));  % inactivation voltage (down)

td = 1/dt;                   % inverse distance time

up = 0;
vth = -20;                   % voltage threshold
cnt = 0;

for j = 2:Nt;

      t(j) = (j-1)*dt;

      Istim = (t(j)-dt/2>5)*stim.amp*1e-6;
      
      a = an(V(j-1));  b = bn(V(j-1)); c = (a+b)/2;
      n = ( (td-c)*n + a ) / (td + c);
      
      cK = G.K*n^4;
      
      a = am(V(j-1));  b = bm(V(j-1)); c = (a+b)/2;
      m = ( (td-c)*m + a ) / (td + c);
      
      a = ah(V(j-1));  b = bh(V(j-1)); c = (a+b)/2;
      h = ( (td-c)*h + a ) / (td + c);
      
      cNa = G.Na*m^3*h;
      
      top = 2*V(j-1)*td + cK*E.K + cNa*E.Na + G.Cl*E.Cl + Istim/A;
      
      bot = 2*td + cK + cNa + G.Cl;
      
      Vmid = top/bot;           % trapezoid scheme

      V(j) = (2*Vmid - V(j-1)); 
     
   if V(j) > vth 
       if up == 0
           up = 1;              % mark upward crossing
           cnt = cnt + 1;
           spst(cnt) = t(j);
       end
   else
       up = 0;
   end

end

isi = inf;
if cnt
   isi = diff(spst);
   if isempty(isi)
       isi = inf;
   else
       isi = isi(end);
   end
end

return

function val = Iss(V,E,G)
val = G.Na*(am(V)/(am(V)+bm(V)))^3*(ah(V)/(ah(V)+bh(V)))*(V-E.Na) + ...
      G.K*(an(V)/(an(V)+bn(V)))^4*(V-E.K) + G.Cl*(V-E.Cl);

function val = an(V)
val = .01*(10-(V+71))./(exp(1-(V+71)/10)-1);

function val = bn(V)
val = .125*exp(-(V+71)/80);

function val = am(V)
val = .1*(25-(V+71))./(exp(2.5-(V+71)/10)-1);

function val = bm(V)
val = 4*exp(-(V+71)/18);

function val = ah(V)
val = 0.07*exp(-(V+71)/20);

function val = bh(V)
val = 1./(exp(3-(V+71)/10)+1);
