% function spnet
%
% NSCI 850 - Computational Approaches to Neuroscience
%
% Spiking network with axonal conduction delays and STDP
%
% Created by Eugene M.Izhikevich.                February 3, 2004
% see also: Izhikevich, Neural Computation 2006
%
% adapted by G. Blohm (2010)
%-------------------------------------------------------------------------

clear

%% initialization
M = 100;                  % number of synapses per neuron 
                          % since N = 1000, the probability of connection
                          % is M/N = 0.1 (similar to neocortex)
D = 20;                   % maximal conduction delay

Ne = 800;                 % excitatory neurons (regular spiking)
Ni = 200;                 % inhibitory neurons (fast spiking)
N = Ne+Ni;                % total number of neurons

% neuron model parameters
a = [0.02*ones(Ne,1);    0.1*ones(Ni,1)]; % excitatory neuron parameter: 0.02
                                          % inhibitory neuron parameter: 0.1
b = 0.2;                  % parameter b
c = -65;                  % resting potential parameter
d = [   8*ones(Ne,1);      2*ones(Ni,1)]; % value 8 and 2 for excitatory and inhibitory neurons
sm = 10;                  % maximal synaptic strength

post = ceil([N*rand(Ne,M);Ne*rand(Ni,M)]);% indices of post-synaptic targets
                                          % inhibitory neurons ONLY connect
                                          % to excitatory neurons!
                                          % (Ne*rand)
s = [6*ones(Ne,M);-5*ones(Ni,M)];         % corrsponsing synaptic weights
sd = zeros(N,M);                          % their derivatives

% set up synaptic conduction delays
for i = 1:N
    if i <= Ne
        for j = 1:D
            delays{i,j} = M/D*(j-1)+(1:M/D); % list of synaptic connections
                                             % from neuron i having delay j
                                             % delays vary between 1-20ms
        end;
    else
        delays{i,1} = 1:M; % all inhibitory connection delays are 1ms
    end;
    pre{i} = find(post==i & s>0);           % list of all excitatory neurons that are pre-synaptic to neuron i
    aux{i} = N*(D-1-ceil(ceil(pre{i}/N)/(M/D)))+1+mod(pre{i}-1,N); % table of indices used to speed up STDP implementation
end;
STDP = zeros(N,1001+D);
v = c*ones(N,1);                            % initialize membrane potential of neuron
u = b.*v;                                   % initialize recovery variable of neuron
firings = [-D 0];                           % spike timings

%% start simulation
for sec = 1:60*60*24                        % simulation of 1 day
    for t = 1:1000                          % simulation of 1 sec (1ms time steps)
        I = zeros(N,1);
        I(ceil(N*rand)) = 20;               % random thalamic input to one neuron
        fired = find(v>=30);                % indices of fired neurons: 30mV = spike apex
        v(fired) = c;
        u(fired) = u(fired)+d(fired);
        STDP(fired,t+D) = 0.1;              % each time a neuron fired, STDP is set to 0.1
        for k = 1:length(fired)
            sd(pre{fired(k)}) = sd(pre{fired(k)})+STDP(N*t+aux{fired(k)}); % if pre-synaptic neurons have spiked BEFORE post-synaptic neuron,
                                                                           % synaptic weights are potentiated
        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)'; % input from the other neurons is added to  the random thalamic input
            sd(firings(k,2),del) = sd(firings(k,2),del)-1.2*STDP(ind,t+D)'; % synaptic depression for pre-synaptic spiking AFER post spiking
                                                                            % depression is stonger (by factor 1.2) than potentiation
            k = k-1;
        end;
        v = v+0.5*((0.04*v+5).*v+140-u+I);    % update of v
        v = v+0.5*((0.04*v+5).*v+140-u+I);    % for numerical stability time step is 0.5 ms
        u = u+a.*(b*v-u);                     % update of u
        STDP(:,t+D+1) = 0.95*STDP(:,t+D);     % decay of STDP strength over time: tau = 20 ms
                                              % determines the magnitude of
                                              % potentiation or depression
    end;
    
    % plot results
    plot(firings(:,1),firings(:,2),'.');
    axis([0 1000 0 N]); 
    title(['Result for time = ' num2str(sec) ' sec'])
    xlabel('time (ms)')
    ylabel('neuron #')
    drawnow;
    
    % re-assign variables for next 1000ms of simulation
    STDP(:,1:D+1) = STDP(:,1001:1001+D);
    ind = find(firings(:,1) > 1001-D);
    firings = [-D 0;firings(ind,1)-1000,firings(ind,2)];
    
    % update of synaptic weights using STDP
    % inhibitory weights are not plastic, whereas excitatory weights evolve
    s(1:Ne,:) = max(0,min(sm,0.01+s(1:Ne,:)+sd(1:Ne,:))); % only excitatory weights are updated! weights cannot exceed max value sm
                                                          % value 0.01 added prevents from synaptic weight decay to zero
    sd = 0.9*sd; % decay of synaptic weight change magnitude over time
end;