
                                          
% The network consists of N=1000 neurons with the first Ne=800 of excitatory RS type, and the remaining Ni=200 of inhibitory PS type.

% Each excitatory neuron is connected to M=100 random neurons, so that the probability of connection is  M/N=0.1. Each inhibitory neuron is connected to M=100 excitatory neurons only.

Ne=800; Ni=200; N=Ne+Ni;                  % excitatory neurons; inhibitory neurons; total number
M=100;                                    % number of synapses per neuron
D=20;                                     % maximal conduction delay

                                          % Each synaptic connection has a fixed integer
                                          % with conduction delay beeing between 1 ms and D=20 ms,
                                          % where D is a parameter (M/D must be integer)

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)]);  % The indices of postsynaptic targets are in the 
                                          % NxM-matrix >>post<<. That is to say,
                                          % <specify each neuron connect to which neurons> 
                                          % Note that7/2
                                          % ceil(N*rand) is a random number form 0 to N-1

sm=10;                                    % maximal synaptic strength                                          
s=[6*ones(Ne,M);-5*ones(Ni,M)];           % Synaptic weights are in the NxM matrix >>s<<

sd=zeros(N,M);                            % their derivatives

for i=1:N
   if i<=Ne
      for j=1:D                      
      delays{i,j}=(M/D)*(j-1)+(1:M/D);    % The list of all synaptic connections from neuron i
      end;                                % having delay j is in the cell array delays{i,j}.
   else                                   % All inhibitory connections are assigned 1 ms delay
      delays{i,1}=1:M;                    % and 1 to D ms delay to all excitatory connections.
   end;

% D different values of delays are assigned to a totla of M connections
% Therefore each of the D values is given to M/D connections, such as
% delays{1,1} = 1 2 3 4 5
% delays{1,2} =           6 7 8 9 10
% delays{1,2} =                      11 12 13 14 15
% etc
                                          
% delays is a cell array, delay(i,j) is the cell containing element in the (i,j)-location
% while delay{i,j} is the element itself. E.g., assume cell = {2+3i, 'akira'; [1:3]', hilb(5)] 
% then cell(1,1) = {2+3i} and cell{1,1} = 2+3i. So, cell(1,1)^5 is wrong but c{1,1}^5 is correct

   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);
end;

                                          % Matrix >>pre<< contains indeces of all excitatory 
                                          % presynaptic neurons to a given neuron i.

                                          % Matrix >>aux<< is an auxiliary table of indices 
                                          % needed to speed up STDP.

% 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


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*60*24                            % simulation of 1 day
   for t=1:1000                               % simulation of 1 sec
      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;
      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;
   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;