Nonhomogenous Poisson spiking train in a spiking neural network framework

nonhomogenousPoisson.m

%% Non-homogenous Poisson spike train in a spiking neural network framework
% We generate inhomogenous Poisson spike trains for a set of output neurons
% in a spiking neural network. We use the "thinning" approach introduced in
% "Simulation of nonhomogeneous poisson processes by thinning" by Lewis and
% Shedler (1979).

% Number of output neurons
numOut = 10;

% Define a supremum for lambda(t). I.e., this should be the smallest upper
% bound on the values that lambda(t) can take.
lambdaMax = 1;

% Runtime of simulation
tSim = 2500;

% Initialize (but note that we don't need to store these because we only
% depend on their values at the current time t)
lambda = zeros(numOut, tSim);
u = zeros(numOut, tSim);

% Inhibitory function
I = 2000 * exp(-5 * (0: tSim + 4)') + 550;

% Initialize a time counter for the inhibitory function
tInh = 5 * ones(numOut, 1);

% Initialize the time of next spike for the output neurons
tNext = tSim * ones(numOut, 1);

for t = 1: tSim
    % Generate EPSP u(t)
    % YOUR CODE HERE
    
    % If there is a spike, reset the time of the inhibitory function for
    % the neurons that did not spike
    if any(tNext == t)
        % Find the index of the first neuron that is spiking this time t
        ind = find(tNext == t, 1);
        
        % Reset the time of the inhibitory function on the other neurons
        tInh([1:ind-1, ind+1:10]) = 1;
    end
    
    % Poisson rate at time t for each of the output neurons
    lambda(:, t) = exp(u(:, t) - I(tInh));
    
    % Find the indices to update: if the next spiking time is the default
    % initialization or the current time, then we have to calculate a new
    % next spiking time.
    update = tNext == tSim | tNext == t;
    
    % Find the next time of spike assuming a homogenous Poisson rate 
    % lambdaMax. This uses inverse transform sampling of the exponential 
    % distribution to find the time interval to the next Poisson event. 
    % Note that the time between Poisson events is exponentially 
    % distributed.
    tNext(update) = t - log(rand(sum(update), 1)) / lambdaMax;
    
    % Reject the next time of spike if a Unif ~ [0, 1] random number is
    % greater than the ratio lambda(t)/lambdaMax
    tNext(update & (rand(numOut, 1) > lambda(:, t) / lambdaMax)) = tSim;
    
    % Perform postsynaptic potential weight update
    if any(tNext == t)
        % YOUR CODE HERE
    end
    
    % Update time on inhibition function
    tInh = tInh + 1;
end