%
% solve the reduced hodgkin huxley equations 
% with an autosynapse
%
% usage   v = autosyn(dt,NT,gbars,delay,VT)
%
%             dt = timestep
%             NT = number of timesteps
%             gbars = magnitude of synaptic conductance
%             delay = synaptic delay
%             VT = synaptic threshold
%
%  e.g.,  v = autosyn(0.01,5000,0.01,8,50)
%

function v = autosyn(dt,NT,gbars,delay,VT)

% parameters from table 3 on page 520 of HH

gbarK = 36;             % mS / (cm)^2
gbarNa = 120;           % mS / (cm)^2
gbarl = 0.3;            % mS / (cm)^2
VK = -12;               % mV
VNa = 115;              % mV
Vl = 10.613;            % mV
C = 1;                  % micro F / (cm)^2

%VT = 50;	% threshold
%delay = 8;	% synaptic delay
Vsyn = 60;  % synaptic reversal potential
%gbars = 1e-2;	% gbarsyn

v(1) = 0;
n = 0.3177;
gsyn(1) = 0;

for j=2:NT,

    t = (j-1)*dt;

    n = (n + dt*an(v(j-1)))/(1+dt*(an(v(j-1))+bn(v(j-1))));

    minf = am(v(j-1))/(am(v(j-1))+bm(v(j-1)));

    Istim = 0;
    if (t>=1 & t<=2)
       Istim = 20;
    end

    vback = 0;
    if t > delay,
       vback = v(j-1-delay/dt);
    end
  
    gsyn(j) = gbars*(max(VT,vback)-VT); 

    num = v(j-1) + dt*(gbarNa*minf^3*(0.8-n)*VNa +...
                       gbarK*n^4*VK +...
                       gbarl*Vl +...
                       Istim + gsyn(j)*Vsyn );

    den = 1 + dt*(gbarNa*minf^3*(0.8-n) + gbarK*n^4 + gbarl + gsyn(j));

    v(j) = num/den;

end

figure(1)
plot(dt:dt:NT*dt,v)
xlabel('t (ms)','fontsize',16)
ylabel('v (mV)','fontsize',16)

head = ['gbar_{syn} = ' num2str(gbars) '   delay = ' num2str(delay) '   V_T = ' num2str(VT)];
title(head,'fontsize',16)

figure(2)
plot(dt:dt:NT*dt,gsyn)
xlabel('t (ms)','fontsize',16)
ylabel('g_{syn} (mS)','fontsize',16)

return

% rate functions from page 519 of HH

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

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

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

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