Event / spike triggered average

Download
evTrigAvg.m

What caused an event?
Generally, event triggered averages are used to try and determine what caused an event to happen. They look at what preceded a number of events and take the average. In neuroscience events might be a spike recorded from a neuron in response to a certain stimulus. If the stimulus was white noise, for example, taking the time window from before the spike and averaging it would give the average stimulus that caused the neuron to fire. In finance the event could be a price change, and one might like to know if a consistent pattern occurs before similar price changes in order to react to them in the future.

evTrigAvg is designed to be general and agnostic to exactly what events and stimuli its dealing with are. All it needs is the "stimulus" and a logical vector of events (0s indicating no events, 1s indicating an event has occurred), and a window size to average over (in points).

Example
Let's generate some fake neural data as an example. Imagine an experiment where white noise is delivered to a subject and activity is recorded from a single neuron. The raw activity from the neuron will be noisy and will normally require preprocessing first, which will include filtering and event (spike) detection. Event detection is often done using a threshold of some value *the root mean square of the trace, where anything above this is marked as a spike. Or, if multiple units have been recorded, by much more complicated "spike sorting". I mention this mainly because it presents a fascinatingly complex problem, but it's not something we need to worry about here. The goal of preprocessing is to obtain a clean vector of events represented as 0s or 1s, we'll imagine the preprocessing has already been done, and we have this clean vector of events (spikes). We want to know what property of the stimulus caused each event.

nPoints = 10000;
stim = randn(1,nPoints);
resp = stim>2;
resp = [zeros(1,100), resp(1:end-100)];

nPoints determines the number of points in the stimulus (stim) and response (resp). stim is random noise drawn from a normal distribution. resp will be our vector of events (spikes), and for the sake of demonstration spikes are generated whenever the stim is greater than the entirely arbitrary value 2. Obviously, this means we know what stimulus property that is eliciting the spikes already, but pretend this is real data that we don't yet know anything about it. The line stim>2 returns a logical vector of 0s and 1s, which we'll also offset slightly relative to stim to emulate the latency that would exist in a real system (spikes aren't instant) by dumping 100 zeros at the start.

Let's plot out the first 1000 points of data we've just magically generated without having to do any pesky experiments.

subplot(2,1,1), plot(stim(1:1000))
line([1, 1000], [2, 2], 'LineStyle', '--', ...
    'color', 'k')
ylabel('Mag, arb'), title('Stimulus')
subplot(2,1,2), plot(resp(1:1000), 'r')
ylabel('Spikes'), title('Response')
xlabel('Time, pts')

Now that we have the data and response we can take the spike triggered average. For each spike we want to look in stim at the preceding time window and take an average of this window. To do this in Matlab, we need a function that will:

• Discard the spikes in the first time window - if we chose a time window of 200 points, and there's a spike at resp(100), we can't see the whole stim that elicited it because we'd need to try and access stim(-100:100), which is impossible.
• Count the total number of events
• Find the index of each event so we can work out the indexes of the time window preceding each spike
• Then, for each spike
• Work out the index range of its preceding time window
• Collect the data in stim from this time window
• Finally, average the data collected from all the time windows

Calculating the average
Let's use a time window of 300 points (windowSize=300) and the stimulus (stim) and response (resp) from above.

windowSize = 300;
resp(1:windowSize)=0;

First windonwSize is set to 300 and the first 300 points (1:windowSize) of resp are then set to 0 to remove the events from the first time window.

nEvs = sum(resp);
evIdx = find(resp);

Next, because response is a vector of 0s and 1s, the number of events is simply the sum of resp. The locations (indexes) of each event can be found using the find function. find simply returns the indexes of all the 1s it finds in resp and these are stored in evIdx, which is a vector listing each index. Note that find(resp,1) is not the correct command to use here, as it would just return the first index of the first 1 found in resp.

We're going to save the average in the variable avg, so let's preallocate it before the for loop.

avg = zeros(1,windowSize);

There are two ways of taking the average at this point - collect all the time windows and take the mean after the for loop, or keep a sum of the data in the windows that is added to in each iteration of the for loop and then divided by the number of iterations after the for loop. The second way is much more efficient, as avg needs to only be 1 x windowSize rather than nEvs x windowSize. It doesn't matter so much here, but this is an important consideration because most real neural data is HUGE.

% For each event
for w = 1:nEvs
    % Find the indexes of the time window preceding the event
    wIdx = evIdx(w)-windowSize : evIdx(w)-1;
    % Add the stim from this window to the average
    avg = avg + stim(wIdx);
end

Now we know the number of events and their indexes, we can start working through them and collecting the time windows from stim. This loop runs from 1:nEvs. On the first iteration (w=1) it creates a vector of indexes for the time window based on the first event index in evIdx. We want the window preceding the event, so say evIdx(w) = 400, the preceding 300 point window would be 100:399, which is calculated with the line wIdx = evIdx(w)-windowSize : evIdx(w)-1 - substituting in each value gives 400-300 : 400-1, so wIdx = 100:399. On the next iteration, if evIdx(w) is 500, the indexes calcuated would be wIdx =500-300 : 500-1, or 200:499, and so on.

At this point wIdx contains a vector of indices that we can use to collect data from stim. The line avg = avg + stim(wIdx); takes this window from stim and sums it to our preallocated variable avg.

avg isn't really an average yet, it's a sum. So the final step is finish the average (mean) by diving avg by the number of observations that contributed to it:

avg=avg./sum(resp);

Note the ./ operator, this denotes element-wise division where each value in a vector (or matrix) is divided by a corresponding element. In this case Matlab knows the division should be element-wise anyway, because we're dividing a vector by a scalar (it can't be matrix division if the sizes are different). Even if just a / were used, Matlab would still do element-wise division, but it's good practice to still write the intended operation.

Let's look at the result...

figure
plot(avg)
title(['Average of ', num2str(nEvs), ' windows'])
xlabel('Time, pts')
ylabel('Mag, arb.')

Clearly not the most exciting spike triggered average, but at least it's what we expect to see given that the criteria we specified for generating the spikes in the first place was "any point the magnitude of the stimulus goes above 2 units, 100 points before the event".

evTrigAvg.m
function [avg, nEvs] = evTrigAvg(events, stim, windowSize)

evTrigAvg.m implements the above process. It expects 3 inputs, a logical (0s and 1s) vector of events (events), the evoking stimulus (stim) and the size of the window to average over in points (windowSize). It returns two outputs, the averaged window (avg) and the number of events/windows that contributed to the average (nEvs).

For example, the command to run the example above would be:

[avg, nEvs] = evTrigAvg(resp, stim, 300);