%Point Neuron 10 pgm three neurons feedback from spikes


%amp=200, ampGK=15 == slower overlaping bursts
%amp=200, ampGK=5 == faster bursts
%amp=40, ampGK=1.2 == faster bursts
%amp=40, ampGK=1.2, tGK=50 == slower bursts weird order



clf
clear


%membrane
tMemb = 5;  %mSec 5


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


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


%input current
amp = 40;    % 200 will cause bursting
tstart = 1;
tend = 500;
weight = -[1, 1, 1];


%init state variables
Vmemb = [0 0 0];
Threshold = initThreshold * [1 1 1];
Spike = [1 0 0];
GK = [0 0 0];


%Time Vector
step = 1;
time = 1:step:tend+50 ;


%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*[1 1 1] ; 
spikeStart=[0 0 0]; %initialize the spike state machine


for t=time
   
   for n=1:3
   
      %input current from last cell
      if ( (n==1 & Spike(3) == 1) | ...
            (n==2 & Spike(1) == 1) | ...
            (n==3 & Spike(2) == 1) )
         current = amp * weight(n);
      else
         current = 0;
      end
      
      %update state variables
      GK(n) = GK(n)*constGK + ampGK*Spike(n)*(1-constGK);
      
      constMemb(n) = exp(-(1+GK(n))*step/tMemb);
      Vmemb(n) = Vmemb(n)*constMemb(n) + ...
         (current + GK(n)*EK)*(1-constMemb(n))/(1+GK(n));
      
      Threshold(n) = Threshold(n)*constTH + ...
         (ampThreshold*Vmemb(n) + ...
         initThreshold)*(1-constTH);
      
      %If the V crosses throeshold then hold 
      %the spike high for a constant time of
      %about one mSec.
      if (Vmemb(n)>=Threshold(n) & spikeStart(n)==0)
         Spike(n) = 1;
         spikedur(n) = step;
         spikeStart(n) = 1;
      elseif (spikedur(n)<.9)
         spikedur(n) = spikedur(n) + step;
      else
         Spike(n) = 0;
         spikeStart(n) = 0;
      end
      
      %form plotting variables
      Vout(tindex,n) = Vmemb(n);  
      GKout(tindex,n) = GK(n);
      Sout(tindex,n) = Spike(n);
      THout(tindex,n) = Threshold(n);
            
   end %of cell setp
   
   Tout(tindex) = t;
   tindex = tindex+1;
   
end %main time step for-loop


hold on
plot(Tout,Sout(:,1),'r')
plot(Tout,Vout(:,1)/40,'r--')
plot(Tout,THout(:,1)/40,'r-.')
%plot(Tout,GKout(:,1)/5,'r.')


plot(Tout,Sout(:,2)-1.5,'g')
plot(Tout,Sout(:,3)-3,'b')


xlabel('time (mS)')
legend('cell 1','Vm','threshold', 'cell 2', 'cell 3')

