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

Answer 1

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.