Spike Triggered Analysis (STA)

This post’s notes comes from this paper, code is done in MATLAB:

Sharpee TO (2013) Computational identification of receptive fields. Ann Rev. Neurosci.

Receptive field (RF)–> optimal stimulus for the neuron

How similar is a given stimulus to the optimal stimulus? This can be answered by using a linear model to predict the firing rate (FR) of a neuron. This means taking the dot product of the stimulus and the RF, or the projection of the stimulus onto the RF.


If we assume that only the stimulus that is projected along the RF affects the FR we can look at where spikes occurred and take some time window (say like 5-10ms) before the spike occurred in the stimulus, doing this for each spike and averaging the stimulus.  The idea is that parts that are in directions that aren’t along the RF will average out. Finally, if the mean is nonzero, subtract the mean from the spike triggered average to get an estimate of the RF.


If the stimulus is circularly symmetric, STA will converge asymptotically to the RF, and ‘STA will not give a correct estimation of relevant stimulus features when stimulus values are correlated across different dimensions’. ‘In a linear model the effect of correlations is limited to pairwise correlations’. Multiplying STA by the inverse of the covariance matrix will give you relevant stimulus features. This step is known as ‘deconvolution’.

Let’s look at an example.

Design Stimulus: Generate a 1D Gaussian white noise stimulus with a time bin of 2 ms and t ranging from 0 to a T of your choosing.

T = 10000;
t = 0:2:T; 
x = randn(1, length(t) + 25);

Simulate Data: Generate synthetic data from a model neuron with a linear temporal filter and a spiking nonlinearity. Define the neuron’s linear temporal filter by:

w = 0.3 % rad/ms
Tau = 10 % ms
f = @(t)exp(-t/Tau)*sin(w*t); % where t is defined from 0 to 50 ms


To get the excitatory drive, u(t), or the projection of the stimulus on to the filter, we take the dot product of the stimulus and the filter. This step is known as convolution.

u = zeros(size(t));
for i = 1:length(u)
    for j = 0:50/2
        uf(i) = uf(i) + s(i-j+25)*f(j+1);

Last, firing rate is a function of excitatory drive, u(t).

rmax = 3.5; % Tuned so that neuron spikes around 20 Hz
theta = 5; delta = 1;
r = rmax. / (1 + exp((theta-u)/delta));
y = poissrnd(r); % This is your spike data -- your generated response

Compute the STA – aka RF – aka can we recover the linear filter?

% Spike triggered average on simulated spike train

ws = 25; % 50 ms -> 25 pts each 2ms, window size
data(1:ws)=0; % set the first window = 0
nev = nnz(data); % number of events
ind = find(data); % indices of events

% Calculating the STA: Spike Triggered Average
sta = zeros(1, ws);
for w = 1:nev
windx = ind(w)-ws: ind(w)-1;
sta = sta + data(ind(w))*x(windx + 25);

sta = sta/nev;


Now consider the case where u(t) is a quadratic non-linear term with respect to firing rate, r, such that:

theta = 15; delta = 2;
r = rmax. / (1 + exp((theta-(u)^2)/delta));

When we simulate a spike train and perform STA, we are unable to recover the linear filter f(t).

Compute the spike triggered covariance and perform and eigen-decomposition of the covariance matrix. This should recover model parameters.

% Calculate Covariance of the spike data

Cspike = zeros(ws);
for w = 1:nev
    windx = ind(w) - ws: ind(w) - 1; 
    Cspike = Cspike + x(windx + 25).'*x(windx + 25);
Cspike = Cspike./nev;
[vectors, values] = eig(Cspike);

Plotting the first eigenvector, we get:


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s