function poissonTutorial %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % poissonTutorial.m % % NOTE: For compatibility with different web browsers, % this file has been renamed 'poissonTutorial.txt'. In % order to use it in MATLAB, you should rename it as % 'poissonTutorial.m' and then open it up with the % MATLAB editor. % % This tutorial introduces the Poisson model of stochastic % neuronal firing. It will demonstrate the type of variability % seen in the response of sensory neurons, and explore the % properties of that variability % % This tutorial is closely based on a tutorial created by % David Heeger for the Cold Spring Harbor course on Computational % Neuroscience: Vision, which was subsequently extended by % Greg Horwitz. % % INSTRUCTIONS: % % While other approaches are possible, it is recommended that you % put a breakpoint on the first line of executable code below, and % then use the Matlab debugger to step through the lines while % reading the text that is interspersed with code. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; % clear, but don't move the windows %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % INTRODUCTION: % % It is not important that you memorize the equations for a Poisson % process, but you should know something about its general properties % how they relate to neuronal spiking. % % A Poisson process is a member of the class of stochastic (random) % point (binary) processes. It is formally defined by the following % set of equations (here, N(t) is the number of events that have occurred % by time 't'). % % 1) N(0)=0 % This means that there are no spikes at the begining of time. % 2) The process is stationary % This means that the rate of spike firing doesn't change over time, % although neuronal responses can be modeled using Poisson processes with % different rates during different epochs (e.g. background and stimulus periods). % 3) The process has independent increments % This means that the number of spikes occurring in any interval is independent % of the number of spikes occurring in any other (non-overlapping) % interval. That is, what happened in the previous interval doesn't affect % the current interval % 3) Prob{N(h)=1} = lamba*h + o(h) % Translation: The probability of one event in a "brief" window of % time is a constant, "lambda" multiplied by the duration of the window, "h". % "o(h)" is the probability of more than one spike in a brief window, and it % is vanishingly small. That is, spike counts in brief intervals are all zeros % and ones. For neurons, this applies if we use ~1 ms intervals. % % The Poisson model occurs very naturally in situations in which a % large number of independent components all have the capacity % for generating an event, but none of them are very likely to do so % at a given time. For example, it does an excellent job of describing % radioactive decay. % % Early neurophysiologists were struck by the irregularity spike trains % recorded from neurons in vivo and thus the Poisson model of neural % firing was born. Even these early neurophysiologists were aware % that the firing rates of real neurons are not stationary and don't have strictly % independent increments. Still, these violations aside, the Poisson % model does a remarkable job of describing firing variability of cortical % neurons. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SIMULATING POISSON SPIKING: % % One of the easiest ways to generate a Poisson spike train is to rely % on the approximation: % % Prob{1 spike during (t, t + timeStepS)} = r(t) * timeStepS % % where r(t) is the instantaneous firing rate (in spikes per second) and timeStepS % is the time step. The right half of equation is essentially "lamba * h" from % equation 3 above, with the assumption that "o(h)" is zero -- timeStepS is short % enough so that there is no chance that the neuron would fire more than one spike % in any given interval. % Let's begin by choosing a time step and by choosing an average % firing rate. timeStepS = 0.001; % 1 msec spikesPerS = 50; % 50 spikes per second, on average durationS = 1.000; % 1 sec simulation times = 0:timeStepS:durationS; % a vector with each time step % Now we choose random numbers, one for each time step, unformly distributed between % 0 and 1. We will use these to decide whether a spike has occurred at each time step. vt = rand(size(times)); % Finally, create a vector of ones and zeros depending on whether the probability of firing % (spikesPerS*timeStepS) is greater than the corresponding random number. We % have a function rasterPlot() that will plot the result as raster, using lines % to show which time steps contain a spike. This should appear in Figure 1, % which you can move to a convenient position on your screen. You might want to % expand the window to see spikes that closely spaced. spikes = (spikesPerS*timeStepS) > vt; rasterPlot(spikes, timeStepS); % The main point is that spikes occur at highly irregular times. You might seen % some signs of a pattern, but there is no reliable pattern in the occurrence of % spikes. If we make a new set of spike times, and plot that, you will see that % it is different. To make things simpler, we have a function makeSpikes() that % will create spike trains in the manner described above, taking arguments for % the time step, spike rate, and sample duration. The function is at the end of this % file if you want to see the code. If you run the next two lines % of code, a new spike train will appear in Figure 1. spikes = makeSpikes(timeStepS, spikesPerS, durationS); rasterPlot(spikes, timeStepS); % Because the Poisson process is stochastic, each spike train will look % different, even though they have the same average number of spikes. The % variability is more easily seen if we display a set of spike trains together % in the plot, in format know as a raster plot (each spike train is one line in % the raster). spikes = makeSpikes(timeStepS, spikesPerS, durationS, 25); rasterPlot(spikes, timeStepS); % Your eye might pick out patterns in the raster, but there is not reliable % structure. Each spike is independent of others in its train, and each spike % train is independent of all others. When testing conditions are held constant % in neurophysiological experiments, the spiking of individual neurons often % looks like this %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % INTER-SPIKE INTERVAL DISTRIBUTIONS: % % An interesting property of the Poisson process is that the interspike % intervals are distributed exponentially. This is a straightforward % consequence of the definition of the Poisson process. If the probability of a % spike in each interval is p, the probability of seeing an interspike interval of % 1 interval is p. The probability that there will not be a spike in the first % interval is (1-p), so the probability of seeing an interspike interval of 2 is % (1-p) * p. The probability of seeing an interspike interval of 3 is (1-p)^2 * % p. Because the probability of going an extra interval without a spike % decreases the probability by a factor of (1-p), the probability distribution % for interspike interval will decrease by a fixed fraction for successive % intervals, which will be an exponentially decaying function. % Let's plot the interspike intervals for our Poisson process. First, get the % spike times. durationS = 1000; % get lots of intervals spikes = makeSpikes(timeStepS, spikesPerS, durationS); % make a spike train spikeTimes = find(spikes) * timeStepS * 1000; % get times when spikes occurred (ms) % We can get the inter-spike intervals by subtraction the times of successive % spikes. spikeIntervals = spikeTimes(2:length(spikeTimes)) - spikeTimes(1:length(spikeTimes) - 1); % Compute a histogram of the spike intervals, normalized to unit volume. figure(2); % use Figure 2 binSize = 1; % 1 ms bins x = [1:binSize:100]; intervalDist = hist(spikeIntervals(spikeIntervals < 100), x); intervalDist = intervalDist / sum(intervalDist) / binSize; % normalize by dividing by spike number bar(x, intervalDist); % As promised, the histogram of interspike intervals looks like an exponential % probability distribution. For verificiation, we'll superimpose the exponential % probability density function. y = exppdf(x, 1 / (spikesPerS * timeStepS)); % exponential function, scaled to predicted max axis([min(x) max(x) 0 max(y) * 1.1]); % make sure everything shows on the plot xlabel('Interspike interval'); ylabel('Probability'); hold on; plot(x, y, 'r'); % add exponential function hold off; % The exponential interspike interval distribution for a Poisson process % points to one way that real neurons will not be perfectly Poisson. Real neurons will % have interspike intervals that depart from a perfectly exponential % distribution. % % QUESTION 1: Why will real neurons not have a perfectly exponential % interspike interval distribution? % % Before leaving the interval distribution, it is worth mentioning that although % a Poisson process will have an exponential interval distribution, the % existence of an exponential interval distribution does NOT necessarily mean that it arose % from a Poisson process. You could get an exponential distribution in a % situation where the interval had some special dependence on each other, % which is not allowed for a Poisson process. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SPIKE COUNT VARIANCE: % % The Poisson model of spiking is important because it makes it possible % to explore the variability of neuronal spiking and how that affect the % reliability of neuronal signals. One notable property of a Poisson process % is that for a given spike rate mean of the spike counts over repeated counts % is equal to the variance of those spike counts. By "variance", we mean the % formal statistical definition of variance, which is the square of the standard % deviation of a distribution. NB: The equality of the mean and variance holds % only for spike COUNTS, and is generally NOT the case for the spike RATES or % interspike intervals of Poisson processes. % % This relationship between spike count mean and spike count variance can be % seen by plotting some distributins of spike counts. Let's see what spike % counts are for 1 second counting intervals when the average rate is 25 % spikes/s. timeStepS = 0.001; % set the spike train parameters spikesPerS = 20; % 20 spike/s rate durationS = 1.000; % 1 s simulation startS = 0.001; % count from for the full 1 second endS = 1.000; spikes = makeSpikes(timeStepS, spikesPerS, durationS, 250); % make 250 spike trains counts = countSpikes(spikes, timeStepS); % get the 250 spike counts plotSpikeCounts(counts); % plot the 250 spike counts (this function is below) % The counts are random, but you should see that they range from as low as ~10 % to as high as ~30 -- a broad range around the expected 20 spikes. The plot in % the lower half of the figure has a point showing the spike count mean for this distribution % (x-axis) and the spike variance (y-axis). The point should fall close to the % identity line. % % Let's add some more distributions with mean counts that are lower (2 and ~6). for i = 1:2 spikesPerS = spikesPerS / sqrt(10); spikes = makeSpikes(timeStepS, spikesPerS, durationS, 250); counts(i+1,:) = countSpikes(spikes, timeStepS, startS, endS); end plotSpikeCounts(counts); % The lower plot shows that the equality between variance and mean holds % regardless of what the mean count is. % % The plots of the three distributions show how the Poission distribution differs from the % normal (also known as bell-shaped or Gaussian) distribution. When a Poisson distribution has % a large mean, as with the blue distribution (mean 20 spikes), it is similar to a normal % distribution. For small means, however, the Poisson distribution differs % greatly. The red distribution (mean 2) is not clipped: Poisson distributions % do not include negative numbers (there are no negative counts). When the mean % counts are small, the bins for 0 and 1 contain most of the counts. % % Variance is a somewhat abstract quantity, but you should be able to convince % yourself that the means are approximately equal to the variances by % considering the standard deviations, for which you should have some intuition. % The standard deviation is the square root of the variance, so Poisson distributions % with means of 2, 6 and 20 should have standard deviations of about 1.5, 2.5 % and 4.5. You should see that the red, green and blue distributions have % roughly these standard deviations. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % EQUIVALENCE OF TIME, NEURONS AND RATE FOR SPIKE COUNTS: % % The last point we will treat here is how Poisson processes produce an equivalence between % time, neurons and rate. Because the occurrence of every spike is independent % of every other spike, adding two Poisson distributions will give a % distribution that is also Poisson. Thus, if a neuron is firing an average of 10 spikes a second % and you repeatedly count spikes in pairs of 1 second periods, the distribution of (spikes from period 1) + % (spikes from period 2) will be a Poisson distribution with a mean count of 20 % spikes. If you instead had two neurons that were both firing an average of 10 % spikes a second and you repeated counted spikes from both over 1 second, their combined counts % would also yield a Poisson distribution with a mean of 20. Both distributions % wuold be the same as that obtained by counting for 1 second periods from a % neuron that was firing an average of 20 spikes per second. It doesn't matter % if doubled the mean by counting spike for twice as long, counting from % twice as many (identical) neurons, or from a neuron firing at twice the rate, % the result is the same. % % Let's look at this directly. The following code will plot 4 distributions. % Each includes 250 counts. The first (blue) distribution is spikes counts for 1 s from a % neuron firing an average of 10 spikes/s. The second distribution is spike counts for 1 s % come a neuron firing an average of 20 spikes/s. The third distribution is % spike counts for 2 s from a neuron firing an average of 10 spikes/s. The % final distribution is counts from 1 s combined from 2 neurons firing 10 % spikes/s. You should be able to see that the three final distributions all % have a mean of 20 and are indistinguishable. counts = zeros(4, 250); spikesPerS = 10; durationS = 1.000; spikes = makeSpikes(timeStepS, spikesPerS, durationS, 250); counts(1, 1:250) = countSpikes(spikes, timeStepS); spikesPerS = 20; durationS = 1.000; spikes = makeSpikes(timeStepS, spikesPerS, durationS, 250); counts(2, 1:250) = countSpikes(spikes, timeStepS); spikesPerS = 10; durationS = 2.000; spikes = makeSpikes(timeStepS, spikesPerS, durationS, 250); counts(3, 1:250) = countSpikes(spikes, timeStepS); spikesPerS = 10; durationS = 1.000; spikes1 = makeSpikes(timeStepS, spikesPerS, durationS, 250); spikes2 = makeSpikes(timeStepS, spikesPerS, durationS, 250); counts(4, 1:250) = countSpikes(spikes1 + spikes2, timeStepS); plotSpikeCounts(counts); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This is the end of the Poisson Tutorial. The code below includes only % functions that were used in the tutorial. You should continue by running % neuronalPerformance and following along in the Lesson.pdf in that folder. function counts = countSpikes(spikes, timeStepS, startS, endS) if (nargin < 4) endS = length(spikes) * timeStepS; end if (nargin < 3) startS = 0; end trains = size(spikes, 1); counts = zeros(1, trains); startBin = startS / timeStepS + 1; endBin = floor(endS / timeStepS); for train = 1:trains counts(train) = sum(spikes(train, startBin:endBin)); end end function plotSpikeCounts(counts) t = counts'; figure(3); subplot(2, 1, 1); hist(t, 0:max(max(t)) * 1.1); s = axis; axis([0 s(2) s(3) s(4)]); subplot(2, 1, 2); xlabel('Spike Counts'); m = mean(t); v = var(t); plot(m, v, 'bo'); s = axis; axisLimit = max(s(2), s(4)); axis([0 axisLimit 0 axisLimit]); hold on; line([0, axisLimit], [0, axisLimit]); hold off; axis('square'); xlabel('Mean Spike Count'); ylabel('Spike Count Variance'); end function rasterPlot(spikes, timeStepS) figure(1); times = [0:timeStepS:timeStepS * (length(spikes) - 1)]; axes('position', [0.1, 0.1, 0.8, 0.8]); axis([0, length(spikes) - 1, 0, 1]); trains = size(spikes, 1); ticMargin = 0.01; % gap between spike trains (full scale is 1) ticHeight = (1.0 - (trains + 1) * ticMargin) / trains; for train = 1:trains spikeTimes = find(spikes(train, :) == 1); yOffset = ticMargin + (train - 1) * (ticMargin + ticHeight); for i = 1:length(spikeTimes) line([spikeTimes(i), spikeTimes(i)], [yOffset, yOffset + ticHeight]); end end xlabel('Time (s)') title('Raster plot of spikes'); end function spikes = makeSpikes(timeStepS, spikesPerS, durationS, numTrains) if (nargin < 4) numTrains = 1; end times = [0:timeStepS:durationS]; spikes = zeros(numTrains, length(times)); for train = 1:numTrains vt = rand(size(times)); spikes(train, :) = (spikesPerS*timeStepS) > vt; end end