%Point Neuron 10 pgm from MacGregor


clf
clear


%membrane
tMemb = 5;  %mSec 5


%threshold
ampThreshold = .5;  % range 0 to 1 -- .5
tThreshold = 25;    % mSec 25
initThreshold = 0; % mV above resting 10 -- zero=>oscillator


%potassium cond
ampGK = 5;  % 5
tGK = 5;     %mSec 5
EK = -10 ;    %mV from resting


%input current
amp = 20;    %start at 20
tstart = 1;
tend = 101;
weight = 1

%init state variables
Vmemb = 0;
Threshold = initThreshold;
Spike = 0;
GK = 0;


%Time Vector
step = 1;
time = 1:step:tend+50 ;
Vout = 0*time;  %make vectors for plotting
GKout = 0*time;
Sout = 0*time;
THout = 0*time;
Tout = 0*time;


%calculate some constants for the integrator
constTH = exp(-step/tThreshold);
constGK = exp(-step/tGK);


%step the time
tindex = 1; %pointer for output arrays
spikedur = step ; spikeStart=0; %initialize the spike state machine


for t=time


   %input current rect pulse
   if (Spike == 1)
      current = amp * weight;
   else
      current = 0;
   end


   %update state variables
   GK = GK*constGK + ampGK*Spike*(1-constGK);


   constMemb = exp(-(1+GK)*step/tMemb);
   Vmemb = Vmemb*constMemb + (current + GK*EK)*(1-constMemb)/(1+GK);


   Threshold = Threshold*constTH + (ampThreshold*Vmemb + initThreshold)*(1-constTH);
   %Threshold = initThreshold + (Threshold-initThreshold)*constTH +ampThreshold*Vmemb*(1-constTH);


   %If the V crosses throeshold then hold
   %the spike high for a constant time of
   %about one mSec.
   if (Vmemb>=Threshold & spikeStart==0)
      Spike = 1;
      spikedur = step;
      spikeStart = 1;
   elseif (spikedur<.9)
      spikedur = spikedur + step;
   else
      Spike = 0;
      spikeStart = 0;
   end


   %form plotting variables
   Vout(tindex) = Vmemb;
   GKout(tindex) = GK;
   Sout(tindex) = Spike;
   THout(tindex) = Threshold;
   Tout(tindex) = t;


   tindex = tindex+1;
end %for loop


plot(Tout,Vout,'r')
hold on
plot(Tout,5*GKout,'g')
plot(Tout,Sout,'b')
plot(Tout,THout,'r-.')


xlabel('time (mS)')
legend('Vm', 'gK*5', 'spike', 'threshold')

