

%
% iaf2net.m 
%
% a simple 3 layer 7 cell net composed of iaf2 cells
%
%  example   iaf2net(10,-50,.04,200)
%

function iaf2net(Ie,Vth,dt,T)

EL = -65;
Vreset = -65;
taum = 10;
rm = 1;		% M Ohm/ mm^2
A = .1;         % mm^2;
Rm = rm/A;	%M Ohm
Ie = 10;		%n A
EK = -70;
tausra = 10;
ginc = 50*0.06/rm;

s = taum/dt;
st = tausra/dt;
N = 7;
v = EL*ones(N,1);
gsra = zeros(N,1);
t = 0;
j = 1;

N = 7;
Esyn = zeros(N);
Eex = 0;
Einh = -70;
W = zeros(N);		% synaptic weights
W(3:5,1:2) = 1;  Esyn(3:5,1:2) = Eex;
W(6:7,3:5) = 1;  Esyn(6:7,3:4) = Eex; 
W(6:7,3:4) = 2;
W(6:7,5) = 0.1;  Esyn(6:7,5) = Einh;
W = 10*W;  

S = zeros(N,1);		% index set of spikes
I = S;
tim = -100*ones(N,1);

while t < T

    j = j + 1;
    t = (j-1)*dt;
    if t > 1,
       I(1:2) = Ie;   % drive layer 1
    end

   gsra = st*gsra/(1+st);
 
   gsyn = rm*W*S;

   v = (s*v+EL+Rm*I+rm*gsra*EK+rm*Esyn.*W*S)./(s+1+rm*gsra+gsyn);

   sind = find(v>Vth);	% indices of spiking cells

   tim(sind) = t;	% spike times of spiking cells
   S = (t-tim).*exp(3*(tim-t));    % alpha synapse
 
   v(sind) = Vreset;

   gsra(sind) = (1+st)*gsra(sind)/st + ginc;

   plot(t*ones(size(sind)),sind,'*')	% raster plot of spike times
   hold on

end

axis([0 T 0 N+1])
xlabel('t (ms)','fontsize',16)
ylabel('cell','fontsize',16)
head = ['iaf2net(' num2str(Ie) ',' num2str(Vth) ','];
head = [head num2str(dt) ',' num2str(T) ')'];
title(head,'fontsize',16)

hold off
