
% 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 <either excitatory or inhibitory>, so that the probability of connection is  M/N=0.1. Each inhibitory neuron is connected to M=100 excitatory neurons only.

Ne=8; Ni=6; N=Ne+Ni;           % excitatory neurons; inhibitory neurons; total number
M=4;                           % number of synapses per neuron
D=2;                           % maximal conduction delay (M/D must be integer in the model)

% 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, it specifies connection to which other neurons>. Note that 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<<


% Inhibitoryweights are not plastic, whereas excitatory weights evolve according to the STDP rule. Their derivatives are in the NxM-matrix >>sd<<, though only the NexM block of the matrix is used.


sd=zeros(N,M);


% Each synaptic connection has a fixed integer conduction delay between 1 ms and D = 20 ms, where D is a parameter (M/D (=5 here) must be integer in the model). (We do not model modifiable delays or transmission failures).

% The list of all synaptic connections from neuron i having delay j is in the cell array delays{i,j}. All inhibitory connections are assigned 1 ms delay and 1 to D ms delay to all excitatory connections.

% Although the anatomy of the model is random, reflecting the connectivity within a cortical minicolumn, one can implement an arbitrarily sophisticated anatomy by specifying the matrices post and delays.

% The details of the anatomy do not matter in the rest of the MATLAB code and do not slow the simulation.

% Once the matrices post and delays are specified, the program initializes cell arrays pre and aux. The former contains indices of all excitatory neurons presynaptic to a given neuron, and the latter is an auxiliary table of indices needed to speed up STDP implementation.

for i=1:N
   if i<=Ne
      for j=1:D                      
      delays{i,j}=(M/D)*(j-1)+(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);
end;


% D different values of delays are assigned to a total 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, 3} =                      11 12 13 14 15
.....
% delays{1,20} =
% ......
% delays{ 801, 1} =
% delays{1000,20} =
% delays{}
                                          
% delays is a >>cell array<<, and 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

% try both delay(1,1) and delay{1,1} etc.


% delays
% under Ne=8, Ni=6, M=4, D=2 => delays is 14x2 matrix before STDP is:
% 
% {1 2}     {3 4}
% {1 2}     {3 4}
% {1 2}     {3 4}
% {1 2}     {3 4}
% {1 2}     {3 4}
% {1 2}     {3 4}
% {1 2}     {3 4}
% {1 2}     {3 4}
% {1 2 3 4} {}
% {1 2 3 4} {}
% {1 2 3 4} {}
% {1 2 3 4} {}
% {1 2 3 4} {}
% {1 2 3 4} {}


%sd
% before STDP all elements of sd = 0


% 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.

%pre

% pre is 1x14 Matrix, entries are diferent from run to run
%[] [] [47] 49 17  1  5 20  4  8 21  3  2 44
%           50 18 31 34 30 36 16     7  6 48
%              19       33    32    15 46
%              22       34    35    29
%                                   43

%aux

% aux is 1x14 Matrix too, entries are diferent from run to run
% 6 [] -11 -12  2  1   3 4  7   2   1  5 3   6 
% 8     -8     -9 -6 -13   -6 -10   5  8     7
% 4                  -12       -9  -7 -7    -8
%                    -11          -13      -10


% 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

% If a spike from an excitatory presynaptic neuron arrives at a postsynaptic neuron (possibly making the postsynaptic neuron fire), then the synapse is potentiated (strengthened). In contrast, if the spike arrives right after the postsynaptic neuron fired, the synapse is depressed (weakened).

% initial values
% STDP: time evolution during 1 sec for each of all neurons
% v:    manbrane voltage of each of all neurons
% v is Nx1 Matrix representing membrane voltage of each of all neurons

STDP = zeros(N,1001+D);   
 v = -65*ones(N,1);                             
 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 = 1 msec

for sec=1:10
    for t=1:1000

       I=zeros(N,1);
       I(ceil(N*rand))=20;

% I is random thalamic input => only one out of N neurons is 20 while all else zero

       fired = find(v>=30);                    

% (again) v is Nx1 Matrix each show voltage
% fired is Nx1 Matrix each represents neuron-number of fired neurons

       v(fired)=-65;

% This resets membrane voltage after firing

       u(fired)=u(fired)+d(fired);

% Recall d was specified previously as d=[8*ones(Ne,1); 2*ones(Ni,1)];
% So, elements of u corespond to fired increased by 8 if it's excitatory, and 2 if it's inhibitory for some purposes

       STDP(fired,t+D)=0.1;

%STDP is Nx1003 matrix whose all elements are zero at the beginning. Now if a neuron fired after D msec it turns to 0.1 (Each time a neuron fires, the variable STDP is reset to 0.1.)

% The magnitude of potentiation or depression depends on the time interval between the spikes. 
% Each time a neuron fires, the variable STDP is reset to 0.1. 

%Every millisecond, STDP decreases by 0.95*STDP, so that it decays to zero as 0.1*exp(-t/20) (ms).This function determines the magnitude of potentiation or depression.


          for k=1:length(fired)
             sd(pre{fired(k)})=sd(pre{fired(k)})+STDP(N*t+aux{fired(k)});
          end;

% Previously sd was specified as sd=zeros(N,M) where M is number of synapses per neuron
% pre, 1xN Matrix, contains indeces of all excitatory presynaptic neurons to a given neuron i.
% fired (fired = find(v>=30)) indices of fired neurons, i.e., fired is Nx1 Matrix each represents No. of fired neurons
% aux is cell array which is an auxiliary table of indices needed to speed up STDP implementation.
% sd NxM, pre 1xN, fired Nx1

% So, sd(pre{fired(k)})=sd(pre{fired(k)})+STDP(N*t+aux{fired(k)}) means ...
% sd is defiened as sd=zeros(N,M) at the beginning, and as N=14, M=4 now, sd is 14x4 matrix
% pre is 1xN cell whose i-th entry is indeces of all presynaptic neurons connected to a  neuron i.
% fired is Nx1 Matrix each represents neuron number of fired neurons
% STDP is Nx1003 matrix whose all elements are zero at the beginning. Now if a neuron fired after D msec it turns to 0.1 (Each time a neuron fires, the variable STDP is reset to 0.1.)

% The magnitude of potentiation or depression depends on the time interval between the spikes. Each time a neuron fires, the variable STDP is reset to 0.1. Every millisecond, STDP decreases by 0.95*STDP, so that it decays to zero as 0.1*exp(-t/20) (ms).This function determines the magnitude of potentiation or depression.

% aux is ... aux{i}=N*(D-1-ceil(ceil(pre{i}/N)/(M/D)))+1+mod(pre{i}-1,N);
% So ...
%koko



% length(M) returns the number of rows of matrix M. So length(fired) is the number of neurons fired at this moment


       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

%% Every millisecond, STDP decreases by 0.95*STDP, so that it decays to zero as 0.1*exp(-t/20).


    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;




% Here variable v represents the membrane potential of the neuron, and u represents a membrane recovery variable, which accounts for the activation of K+ ionic currents and inactivation of Na+ ionic currents, and it provides negative feedback to v. After the spike reaches its apex at +30 mV, which is not to be confused with the firing threshold, the membrane voltage and the recovery variable are reset.

% We use (b, c) = (0.2,.65) for all neurons in the network. For excitatory neurons, we use the values (a, d) = (0.02, 8) corresponding to cortical pyramidal neurons exhibiting regular spiking (RS) firing patterns. For inhibitory neurons, we use the values (a, d) = (0.1, 2) corresponding to cortical interneurons exhibiting fast spiking (FS) firing patterns.

% Variable I in the model combines two kinds of input to the neuron: (1) random thalamic input and (2) spiking input from the other neurons. This is implemented via N-dimensional vector I.

% ---------------------------------
% Spike-Timing-Dependent Plasticity
% ---------------------------------
% The synaptic connections in the model are modified according to the spike-timing-dependent plasticity (STDP) rule (Song et al., 2000). We use the simplest and the most effective implementation of this rule,

% If a spike from an excitatory presynaptic neuron arrives at a postsynaptic neuron (possibly making the postsynaptic neuron fire), then the synapse is potentiated (strengthened). In contrast, if the spike arrives right after the postsynaptic neuron fired, the synapse is depressed (weakened).

% If pre- and postsynaptic neurons fire uncorrelated Poissonian spike trains, there are moments when the weight of the synaptic connection is potentiated, and there are moments when it is depressed. We chose the parameters of the STDP curve so that depression is stronger than potentiation and the synaptic weight goes slowly to zero. Indeed, such a connection is not needed and should be eliminated. In contrast, if the presynaptic neuron often fires before the postsynaptic one, then the synaptic connection slowly potentiates. Indeed, such a connection causes the postsynaptic spikes and should be strengthened. In this way, STDP strengthens causal interactions in the network.

% The magnitude of potentiation or depression depends on the time interval between the spikes. Each time a neuron fires, the variable STDP is reset to 0.1. Every millisecond, STDP decreases by 0.95*STDP, so that it decays to zero as 0.1*exp(-t/20) (ms). This function determines the magnitude of potentiation or depression.

% For each fired neuron, we consider all its presynaptic neurons and determine the timings of the last excitatory spikes arrived from these neurons. Since these spikes made the neuron fire, the synaptic weights are potentiated according to the value of STDP at the presynaptic neuron adjusted for the conduction delay. This corresponds to the positive part of the STDP curve. Notice that the largest increase occurs for the spikes that arrived right before the neuron fired, that is, the spikes that actually caused postsynaptic spike.

% In addition, when an excitatory spike arrives at a postsynaptic neuron, we depress the synapse according to the value of STDP at the postsynaptic neuron. This corresponds to the negative part of the STDP curve. Indeed, such a spike arrived after the postsynaptic neuron fired, and hence the synapse between the neurons should be weakened. (The same synapse will be potentiated when the postsynaptic neuron fires.)

% Instead of changing the synapticweights directly,we change their derivatives sd, and then update the weights once a second according to the rule s <- s + 0.01 + sd, and sd <- 0.9*sd, where 0.01 describes activityindependent increase of synaptic weight needed to potentiate synapses coming to silent neurons. Thus, the synaptic change is not instantaneous but slow, taking many seconds to develop.

% We manually keep the weights in the range between 0 and sm, where sm is a parameter of the model, typically less than 10 (mV).