at the end of the day

Simulating Poisson Neuron With Refractory Period

2011-04-11T15:36:00.000-07:00

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

2011-04-11T14:55:00.000-07:00

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

2011-04-11T11:15:00.000-07:00

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

2011-04-10T14:38:00.000-07:00

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

install g77 in Ubuntu 9.04 (Jaunty Jackalope)

2009-10-23T09:13:00.000-07:00

Recently i needed to compile some programs with g77, ubuntu 9.04 repositories has gcc 4.3 and gfortran. Unfortunately, gfortran is not compatible with many g77 dependent programs.
There's an easy way to install gcc 3.4 and g77:
add to /etc/sources.list the ubuntu 8.04 hardy repositories:

deb http://hu.archive.ubuntu.com/ubuntu/ hardy universe
deb-src http://hu.archive.ubuntu.com/ubuntu/ hardy universe
deb http://hu.archive.ubuntu.com/ubuntu/ hardy-updates universe
deb-src http://hu.archive.ubuntu.com/ubuntu/ hardy-updates universe

then in the terminal:

sudo aptitude update
sudo aptitude isntall g77

remove or comment the ubuntu 8.04 repositories in sources.list
that's it