
%
% Gabbiani & Cox, Mathematics for Neuroscientists
%
% moli spike frequency
%
% usage:    molifreq(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:  molifreq(80:2:200)
%

function molifreq(I0)

for j=1:length(I0),
   
    freq(j) = 1000/moli_isi(.01,1000,I0(j));

end

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

%
% moli_isi.m
%
% Molineux on the active sphere
%
% usage:  isi = moli_isi(dt,Tfin,I0)
% 
% where:    dt = time step
%           Tfin = final time
%           I0 = current stimulus
%
%           isi = interspike interval
%
% example:  isi = moli_isi(.01,100,80)
%

function isi = moli_isi(dt,tfin,I0)

Nt = tfin/dt + 1;                       % number of time steps
V = zeros(1,Nt);t = V;Istim = V;        % preallocate space
A = 4*pi*(1e-6);                        % cm^2
Cm = 1.5;
g = struct('K',7,'Na',30,'L',1,'A',16);   
E = struct('K',-90,'Na',45,'L', -70);  
Vr = fsolve(@(v) irest(v,g,E),-120);     % -70.837

j = 1;
V(j) = Vr;
h1 = hinf(Vr);
n1 = ninf(Vr);
b1 = binf(Vr);

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

for j=2:Nt

    tmp = 2*tauh(V(j-1)); 
    h2 = ( h1*(tmp-dt) + 2*dt*hinf(V(j-1)) ) / ( dt + tmp );

    tmp = 2*taun(V(j-1)); 
    n2 = ( n1*(tmp-dt) + 2*dt*ninf(V(j-1)) ) / ( dt + tmp );

    tmp = 2*taub(V(j-1)); 
    b2 = ( b1*(tmp-dt) + 2*dt*binf(V(j-1)) ) / ( dt + tmp );

    a = ainf(V(j-1));
    m = minf(V(j-1));
    
    t(j) = dt*j;
    
    Istim(j) = (t(j)>2)*I0*1e-6;
    
    top = (2*Cm/dt - g.Na*m*h1 - g.K*n1 - g.A*a*b1 - g.L)*V(j-1) + ...
           g.Na*m*(h1+h2)*E.Na + g.K*(n1+n2)*E.K + g.A*a*(b1+b2)*E.K + ...
           2*g.L*E.L + (Istim(j)+Istim(j-1))/A;
    bot = 2*Cm/dt + g.Na*m*h2 + g.K*n2 + g.A*a*b2 + g.L;

    V(j) = top/bot;         % trapezoid scheme
    
    h1 = h2; n1 = n2; b1 = b2;
      
   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);
   isi = isi(end);
end

return

function val = irest(v,g,E)
val = g.Na*minf(v)*hinf(v)*(v-E.Na) + g.K*ninf(v)*(v-E.K) + ...
      g.A*ainf(v)*binf(v)*(v-E.K) + g.L*(v-E.L);

% Na functionals

function val = minf(v)
val = 1./(1+exp(-(v+35)/4));

function val = hinf(v)
val = 1./(1+exp((v+35)/4));

function val = tauh(v)
val = 2*232*28./(4*pi*(v+74).^2+28^2) - 0.15;

% K functionals

function val = ninf(v)
val = 1./(1+exp(-(v+35)/4));

function val = taun(v)
val = 0.5;

% A type K functionals

function val = ainf(v)
val = 1./(1+exp(-(v+27)/8.8));

function val = binf(v)
val = 1./(1+exp((v+68)/6.6));

function val = taub(v)
val = 15;

