%% This is a tutorial for generating Poisson spike trains.
% Original tutorial: David Heeger (for Cold Spring Harbor Lab summer school).
% Modified by: Odelia Schwartz.
% Read the comments and copy each line of code into matlab.
% In Matlab's editor, to run one block of code (starting with %%)
% you can use the Run and Advance button.
% Type 'help ' in Matlab window for any
% function you would like more information on.
% Type 'which ' in Matlab window for the
% location of a file.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generating Samples of a Poisson Process.
% Spike trains in neurons are often modeled as a Poisson random
% process.
%
% One of the easiest ways to generate a Poisson spike
% train is to rely on the approximation that in any small
% time interval the probability of a spike is proportional
% to the mean spike count (firing rate) times the time step:
%
% Prob{1 spike during (t,t+deltaT)} = r(t) deltaT
%
% where r(t) is the firing rate and deltaT is
% the time step. This approximation is only valid when deltaT is
% very short so that there is essentially no chance of there
% being more than one spike per time interval.
%% (1)
% Let's begin by choosing a time step and by choosing an average
% firing rate.
% After running this, try and also choose different firing rates.
deltaT=1e-3 % Time step: 1e-3 secs = 1 msec
rate=10 % Firing rate: spikes/sec
duration=1 % 1 sec simulation
%% Choose the time intervals
times=[0:deltaT:duration]
size(times)
%% We can generate one random number between 0 and 1
% Try this many times
xr_one = rand
%% Let's check if the mean firing rate times the time interval
% is greater than the random number
% Again, we can repeat choosing a random number and checking
% many times
intervalRate = (rate*deltaT)
xr_one
intervalRate > xr_one
%% We want to do this process not just for one interval,
% but for all the time intervals over the duration
% We'll choose a bunch of random numbers, one for each time step,
% uniformly distributed between 0 and 1.
xr=rand(size(times))
%% Check that values are reasonably random by plotting a histogram
figure(1);
hist(xr,20);
set(gca, 'fontsize', 18);
xlabel('random draws')
ylabel('Probability')
%% Now, insert a spike whenever the probability of firing
% (rate*deltaT) is greater than the corresponding random number:
% Matlab allows you to do this operation in one line, without a
% for loop.
neuralResponse = intervalRate > xr
figure(2);
% We'll plot the neural responses (spikes) using the stem command
% which is a lollipop plot...
stem(times,neuralResponse);
set(gca, 'fontsize', 18);
xlabel('Time')
ylabel('Neural Response (Spikes)')
%% To get a better understanding of what this does, we
% can look at the values for a vector of length 5
% Try changing the vector to get more spikes
tmpxr = [.001 .1 .2 .5 .06]
intervalRate
tmpResponse = intervalRate > tmpxr
%% Calculate the spike count, which are the number of spikes during the time duration
spikeCount=sum(neuralResponse)
%% TODO
%% (2) Repeat 5000 times: generating random numbers xr, computing neuralResponse and calculating
% spikeCount. Save the resulting number of spikes into spikeCounts each of the 5000 times.
len = 5000;
spikeCounts = zeros(len,1);
size(spikeCounts)
%% What is the mean of the spikeCounts? How does it
% compare to the firing rate?
mean(spikeCounts)
rate
%% Plot a spike count histogram
figure(3);
x=[0:30];
hist(spikeCounts,x)
xlabel('Spike count')
ylabel('Probability')
%% Plot a normalized histogram (that sums to 1)
% for comparison with the Poisson
spikeCountHist=hist(spikeCounts,x)/len;
y=mypoisspdf(x,mean(spikeCounts));
bar(x,spikeCountHist)
axis([min(x) max(x) 0 0.25])
xlabel('Spike count')
ylabel('Probability')
hold on
plot(x,y)
hold off
%% Here are some things for you to do in class for the code above:
%
% - Try repeating the above code by choosing a lower firing rate
% (3) or a higher rate (15) by changing the variable rate at the top.
% Plot the resulting distributions. How does the mean change?
% Also plot examples of the spike trains in the repeated runs.
%
% - Try changing the number of repeats from 1000, to 50, 10, and 10000.
% Plot the resulting distributions. How good is the fit to Poisson?
%