

Ne=16; Ni=4; N=Ne+Ni;                     % excitatory neurons; inhibitory neurons; total number
M= 9;                                     % number of synapses per neuron
D= 3;                                     % maximal conduction delay

a=[0.02*ones(Ne,1);0.1*ones(Ni,1)];      % used later as u=u+a.*(0.2*v-u);
d=[8*ones(Ne,1);2*ones(Ni,1)];

post=ceil([N*rand(Ne,M);Ne*rand(Ni,M)]);  % this specifies to which M neurons each neuron connects
sm=10;                                    % maximal synaptic strength                                          
s=[6*ones(Ne,M);-5*ones(Ni,M)];           % Synaptic weights in this NxM matrix >>s<<
sd=zeros(N,M);                            % used later

% Each synaptic connection has a fixed integer conduction delay (not a modifiable delays here)
%    between 1 ms and D = 3 ms.
% These delays are assigned not at random but equally distributed, i.e.,
%    M/D out of M connections are assigned the delay value of 1, then 2 to the next M/Â connections, 
%       and so on ..., automatically (Hence M/D must be integer)
%    Namely, connections from neuron 1, or any other neuron, to the neuron 1, 2, 3, ..., N
%       is assigned the delay value 1, 2, 3, ..., N, respectively, if connection exsists

% All inhibitory connections are assigned 1 ms delay  
%    Hence delays{n,1}=1,2,3,...,M  for any n>Ne
%    E.g. >>delays{13,1} = 1 2 3 4 5 6 7 8 9 that means
%       from the 13th neuron delay value 1 is assigned to the 1st, 2nd, ... 9th connection.

% All excitatory connections are assigned one of a delay value from 1 to D ms.
%    each excitatory neuron has M connections to other neurons and 
%       some of the connections has a delay of same value. 
%
% The list of all synaptic connections from neuron i having delay j is in the >>cell array<< delays{i,j}. 
%    E.g. >>delays{3,1} = 1 2 3 means that
%       from the 3rd neuron delay value 1 is assigned to the 1st, 2nd, and 3rd connection.
%    and  >>delays{3,2} = 4 5 6 means that
%       from the 3rd neuron delay value 1 is assigned to the 4th, 5th, and 6th connection.
%    and  >>delays{3,3} = 7 8 9 means that
%       from the 3rd neuron delay value 1 is assigned to the 7th, 8th, and 9th connection.

for i=1:N     
   if i<=Ne
      for j=1:D                      
      delays{i,j}=(1:M/D) + (j-1)*(M/D);
      end;
   else
      delays{i,1}=1:M;
   end;

   pre{i}=find(post==i & s>0)                                        
   aux{i}=N*(D-1-ceil(ceil(pre{i}/N)/(M/D)))+1+mod(pre{i}-1,N);    % for later usage
end;

% pre and aux are also a >>cell array<<
% pre{i} contains indices of all excitatory neurons (s>0) connected to i-th neuron
% aux{i} is an auxiliary table of indices needed to speed up STDP implementation.


STDP = zeros(N,1001+D);
v = -65*ones(N,1);                            % initial values
u = 0.2.*v;                                   % initial values
firings=[-D 0];                               % spike timings

%%for sec=1:60                                % simulation of 1 day
%%   for t=1:10                               % simulation of 1 sec
for sec=1:4
   for t=1:2
      I=zeros(N,1);
      I(ceil(N*rand))=20;                      % random thalamic input
      fired = find(v>=30);                     % indices of fired neurons
      v(fired)=-65;
      u(fired)=u(fired)+d(fired);

      STDP(fired,t+D)=0.1;
         for k=1:length(fired)
            sd(pre{fired(k)})=sd(pre{fired(k)})+STDP(N*t+aux{fired(k)});
         end;
         
% Inhibitory weights are not plastic, while excitatory ones evolve according to the STDP
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

% The model here is
   % dv/dt = 0.04 v^2 + 5 v +140 - u + I
   % du/dt = a(bv - u)
   % auxiliary after-spike resetting
   % if v >= 30 mV then v <- c and u <- u + d
   % (a,b,c,d)=(0.02, 0.2, -65, 8) for excitatory neurons, and 
   % (a,d)=(0.1, 2) for inhibitory neurons



      firings=[firings;t*ones(length(fired),1),fired];
      k=size(firings,1);

      while firings(k,1)>t-D
         del=delays{firings(k,2),t-firings(k,1)+1};
         ind = post(firings(k,2),del);
         I(ind)=I(ind)+s(firings(k,2), del)';
         sd(firings(k,2),del)=sd(firings(k,2),del)-1.2*STDP(ind,t+D)';
         k=k-1;
      end;

      v=v+0.5*((0.04*v+5).*v+140-u+I);                  % for numerical
      v=v+0.5*((0.04*v+5).*v+140-u+I);                  % stability time
      u=u+a.*(0.2*v-u);                                 % step is 0.5 ms
      STDP(:,t+D+1)=0.95*STDP(:,t+D);                   % tau = 20 ms
   end;

%%   firings                                  % for debug
%%   firings(:,1)                             % for debug
%%   firings(:,2)                             % for debug

   plot(firings(:,1),firings(:,2),'.');
   axis([0 1000 0 N]); drawnow;
   STDP(:,1:D+1)=STDP(:,1001:1001+D);
   ind = find(firings(:,1) > 1001-D);
   firings=[-D 0;firings(ind,1)-1000,firings(ind,2)];
   s(1:Ne,:)=max(0,min(sm,0.01+s(1:Ne,:)+sd(1:Ne,:)));
   sd=0.9*sd;
end;

