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.
SPIKE TRIGGERED AVERAGE (STA)
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); end end
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); end 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); end Cspike = Cspike./nev; [vectors, values] = eig(Cspike);
Plotting the first eigenvector, we get:
