Vectorized approach to compute PSTH (peristimulus time histogram) in MATLAB

Question

I have a vector of spike times (action potentials from a neuron) and a vector of stimulus event timestamps. I want to create a PSTH to see if the stimulus influences the spike rate of the neuron. I can do this by looping through each stimulus event (see simple example below), but this is very slow for long experiments where there are over 30,000 stimulus events and many neurons are being recorded.

How can this be done without for loops?

Example of the slow way:

% set variables
spikeTimes = [0.9 1.1 1.2 2.5 2.8 3.1];
stimTimes = [1 2 3 4 5];        
preStimTime = 0.2;
postStimTime = 0.3;
for iStim = 1:length(stimTimes)
    % find spikes within time window
    inds = find((spikeTimes > (stimTimes(iStim) - preStimTime)) & (spikeTimes < (stimTimes(iStim) + postStimTime)));
    % align spike times to stimulus onset
    stimONtimes = spikeTimes(inds) - stimTimes(iStim);
    % store times in array for plotting
    PSTH_array(iStim,1:length(stimONtimes)) = stimONtimes;
end

You might have to tell us what a PSTH does. In a normal histogram you only want the counts for each bin, but in your case it seems you're putting the individual values in each bin. Is this what you want? — beaker, Aug 14 '16 at 15:35
@beaker I'm not putting the values into bins in the example code, I'm merely storing the spike times that occurred in the defined time window for each stimulus presentation. This is what I want to optimize. One could then make a histogram using that array and define time bins of any size. — raz, Aug 14 '16 at 15:44
Ah, I see. That's a shame because it would be easier to do the sum or count or whatever. (Or, at least, I can see a way to do those more directly.) Still, a description would be quite welcome. — beaker, Aug 14 '16 at 15:46
@beaker How would you do it using the sum or count? I would definitely like to know how to ask how many spike times are between 0.8 and 1.3, and how many spike times are between 1.8 and 2.3, without looping. The output should be 3 and 0, and should be stored somehow. Here is a reference for PSTH if you're interested in the application https://en.wikipedia.org/wiki/Peristimulus_time_histogram — raz, Aug 14 '16 at 16:38
Well, you can try `histogram` or `discretize` and pass the bin edges, but I'm not sure what they do when the bins overlap. — beaker, Aug 14 '16 at 17:06
Sorry, that was silly. The bins aren't allowed to overlap. My original thought was to use `bsxfun` to create a `|stimTimes| x |spikeTimes|` matrix, but depending on how large these get, you might run into space problems. — beaker, Aug 14 '16 at 17:57

Uri Cohen · Answer 1 · 2016-08-14T20:37:05.873

Here is a solution with a single loop over all spikes and two assumptions:

stimuli times are in fixed intervals
stimuli intervals are larger than the PSTH interval

Assuming that stimuli times are in fixed intervals:

delta_times = mean(diff(stimTimes));
assert(max(abs(diff(stimTimes)-delta_times))<1e-3);

Now we align spike times to preStimTime before the first stimulus:

spikeTimes0 = spikeTimes - stimTimes(1) + preStimTime;

Now we wish to calculate for each spike the stimuli to which it is counted, using the second assumption:

assert((postStimTime-preStimTime)<dekta_times);
stimuli_index = floor(spikeTimes0 / delta_times);

Calculate relative to that stimuli:

spike_time_from_stimuli = spikeTimes0 - stimuli_index*delta_times;

Now lets build the PSTH, in presicion of 0.01 (in the same units as all other times):

dt = 0.01;
times_around_stimuli = preStimTime:dt:postStimTime;
n_time_bins = length(times_around_stimuli);
n_stimuli = length(stimTimes);
PSTH = zeros(n_stimuli, n_time_bins)
for i=1:length(spikeTimes)
     time_index = ceil(spike_time_from_stimuli(i) / dt);
     % Ignore time-bins far from the event
     if time_index > n_time_bins
           continue;
     end
     PSTH(stimuli_index(i),time_index) = PSTH(stimuli_index(i),time_index) + 1;
end

thanks for the response. would you please test a full version of your proposed solution and then send it to me? there are a number of errors with the code you have pasted here (e.g. stimuli_before_index, n_time_bins, and n_stimuli are all undefined). — raz, Aug 14 '16 at 20:33
I've fixed the problems I've found, but there is no fuller version, I wrote it for your answer, not copy-pasted it. You should first understand it so that you can fix any such issue yourself! Good luck. — Uri Cohen, Aug 14 '16 at 20:42
by full version, i mean that you were able to successfully run something on your computer and get my desired output. when i try to use your code it does not execute without error. i'd rather not debug your code when you don't even know that it works... — raz, Aug 15 '16 at 05:59

beaker · Accepted Answer · 2016-08-14T22:57:45.053

The best way is to probably just use the existing histogram functions. They're very fast and should give you all of the information you need. Of course, this assumes that the bins don't overlap. Given your example data:

spikeTimes = [0.9 1.1 1.2 2.5 2.8 3.1];
stimTimes = [1 2 3 4 5];        
preStimTime = 0.2;
postStimTime = 0.3;

you can construct the bins like so:

bins = sort([stimTimes - preStimTime, stimTimes + postStimTime])

or

bins = [stimTimes - preStimTime; stimTimes + postStimTime];
bins = bins(:).'

bins =
   0.80000   1.30000   1.80000   2.30000   2.80000   3.30000   3.80000   4.30000   4.80000   5.30000

Then you can use histcounts, discretize, or histc, depending on your desired results and which version of MATLAB you have. I'm going to use histc (because I don't have all that fancy stuff) but the inputs are the same for all three functions. You have one extra output for histcounts (the edges, which is useless to us) and one less for discretize (the actual counts).

[N, IDX] = histc(spikeTimes, bins)

N =    
   3   0   0   1   2   0   0   0   0   0

IDX =    
   1   1   1   4   5   5

Since the bins include the times between (T(i) + postStimTime) and (T(i+1) - preStimTime) we need to take every other bin:

N = N(1:2:end)

N =
   3   0   2   0   0

Likewise, we're only interested in the spikes that happened in the odd-numbered timeslots, and we need to adjust the indices to match the new IDX:

v = mod(IDX, 2)

v =
   1   1   1   0   1   1

IDX = ((IDX+1)/2).*v

IDX =
   1   1   1   0   3   3

The results agree with what we got originally: There are 3 spikes in bin 1 and 2 spikes in bin 3.

this changed the execution time from anywhere between 20 mins to a few hours down to **3 seconds**. bravo beaker! — raz, Aug 15 '16 at 15:26

Vectorized approach to compute PSTH (peristimulus time histogram) in MATLAB

2 Answers2