Simulating Poisson Neuron With Refractory Period

In this post, I will again simulate homogeneous Poisson neuron, but with refractory period.

T = 60000;
rate = 40;
refractoryPeriod = 3;
dt = 1;
spikeTimes = [];
isRefractory = false;
restRate = 0;
r = rate;

for t = 1 : dt : 60000,
    if r * dt / 1000 < rand(1),
        spikeTimes(i) = t;
        isRefractory = true;
        r = restRate;
    end
    
    if spikeTimes(end) - t = refractoryPeriod,
        r = rate;
        isRefractory = false;
    end
end

now lets make TIH, Survivor and Hazard functions:

isi = diff(spikeTimes);
maxIsi = ceil(max(isi) / 100) * 100;
bins = -0.5 : maxIsi - 0.5;
isiHist = hist(isi, bins);
isiHist = isiHist / sum(isiHist);

%we need the cumultative sum for the survivor function
cumTih = cumsum(isiHist);
survivorFunc = 1 - cumTih;

%hazard functions is the ISI divided by the survivor function
hazard = isiHist ./ survivorFunc;

Finding Spike Rate Using Convolution

In the previous post, we found spike rate using non overlapping window. In this post,
I will show how to obtain spike rate using convolution with rectangular and Gaussian window.
I will also use the same variables.

%size of the window 100ms
windowSize = 100;

%the area must be 1
kernel = ones(1, windowSize ) ./ windowSize ;

output = conv(spikeTrain, kernel);

you should get something like this:

now if we want the kernel to be Gaussian window:

%width 600ms and variance 2
windowSize = 600;
windowVar = 2;

kernel = gausswin(windowSize, windowVar);

output = conv(spikeTrain, kernel);

you should get something like this:

Spikes Firing Rate

Continuing the last post, say, we would like to plot the firing rate of the spikes using non overlapping window.
I'll use the same variables as in previous post.

%window size 300ms
dt = 300;

%edges for bins
edges = 0 : dt : max(spikeTimes);
r = zeros(1, length(edges));

%we basically count the number of spikes in every bin
for i = 1 : length(spikeTimes),
   for j = 1 : length(edges) - 1,
       if edges(j) <= spikeTimes(i) && 
           spikeTimes(i) <= edges(j + 1),
          r(j) = r(j) + 1;
       end
   end
   
   if spikeTimes(i) == edges(end),
      r(end) = r(end) + 1;
   end
end

%normalizing the firing rate
r = r ./ dt;

%plot the firing rate
stairs(edges, r);

the result should be something like this:

We can now calculate the fano factor:

ff = var(r) / mean(r);

Poisson Spike Train Simulation

Recently I had to simulate spike trains of neuron which fires as a Poisson processes. Since code in neuroscience is mostly written in Matlab, I will stick to this.

%rate of firing
rate = 50;

%dt in milliseconds
dt = 1;

%total time of simulation 60 seconds
totalTime = 60000;

%declaration of variables
spikeTimes = [];
spikeTrain = [];

%loop throough whole time
for t = 1 : dt : totalTime,
    %r*dt is the poisson distribution
    %of firing in dt. if the distribution is 
    %greater the uniform sampling,
    %a spike is fired and its time
    %is kept.
    if rate * dt / 1000 >= rand(1),
        spikeTimes(end + 1) = t;
    end
end

%convert the vector of spike times
%to vector of spike trains.
spikeTrain(spikeTimes) = 1;

now lets plot the interspike interval histogram and see
if the neuron is indeed poissonian.

isi = diff(sp_times);
maxi = ceil(max(isi) / 100) * 100;
bins = 0.5 : maxi - 0.5;
isiHist = hist(isi, bins);
bar(bins, isiHist / max(isiHist), 1);

we've got:

as you can see, the TIH follows exponential distribution,
hence the neuron is indeed Poisonnian.

We can also calculate the coefficient of variation:

%for calculating the Cv we need the interspike intervals
cv = std(isi) / mean(isi);