%% Neural Box
% Odelia Schwartz
% Partly based on Cold Spring Harbor white noise tutorial by Chichilnisky.
% This is a tutorial for figuring out the response properties
% of a neuron by presenting many random stimuli.
% The neuron here is a model simulation, but the
% same approach is used for understanding real neurons.
% We probe the model neuron as if it is a black box and
% we do not know its properties.
% This tutorial specifically focuses on spike-triggered average,
% estimating the linear filter, and estimating the nonlinearity.
%
% At the end of the tutorial, there are some questions.
% ----------------------------------------------------------------------
%% We want to first choose experimental stimuli that are random
% At each frame, the intensity of the uniform screen changes:
% it is drawn randomly, here from a Gaussian distribution.
numsamples = 25000
stimulus=(1/3*randn(numsamples,1))
%% Plot the stimulus
figure(1);
thelen = min(1000, numsamples);
plot(stimulus(1:thelen))
set(gca, 'fontsize', 14);
xlabel('Time (msec)')
ylabel('Stimulus strength')
title('First 1 second of stimulus');
%% Let's check that it is Gaussian
figure(2);
[H, x] = hist(stimulus, 20);
bar(x, H./sum(H));
set(gca, 'fontsize', 14);
xlabel('Stimulus strength');
ylabel('Histogram');
%% We're now going to simulate a model neuron
% For purposes of this demo, we are constructing the model
% neurons and so know their filters and nonlinearity
% (in an experiment with real neurons, we would be handed
% the spike trains and would not know this!)
% We've made several versions of model neurons.
% We have 3 possible linear filters.
% Toggle between
kernelSize = 60;
[linearResp, filter] = getLinear1(stimulus, kernelSize);
%[linearResp, filter] = getLinear2(stimulus, kernelSize);
%[linearResp, filter] = getLinear3(stimulus, kernelSize);
%% Let's look at the filter
figure(3);
h=plot(filter, 'o-');
set(h, 'linewidth', 2);
set(gca, 'fontsize', 14);
title('actual filter')
xlabel('time');
%% We also have two versions of nonlinearities.
nonlinearResp = getNonlinear1(linearResp);
% Toggle between
%nonlinearResp = getNonlinear2(linearResp);
%% We can use this non-linear response to simulate a
% Poisson-ish spike train...
xr=rand(size(nonlinearResp));
neuralResponse = nonlinearResp>.05*xr;
spikeCounts = neuralResponse;
%% So far, we constructed a model neuron and its response to experimental
% stimuli. Here's the first second of each step:
thelen = min(numsamples, 1000);
h = figure(4); clf;
h = subplot(3,1,1);
plot(linearResp(1:thelen),'b-')
title('Linear response')
set(h,'XTick',[],'XTickLabel',[])
h = subplot(3,1,2);
plot(nonlinearResp(1:thelen),'r-')
title('Nonlinear function of linear response')
set(h,'XTick',[],'XTickLabel',[])
h = subplot(3,1,3);
stem(neuralResponse(1:thelen))
set(gca,'Ylim',[0 2]);
title('# of Spikes (1 ms bins)');
set(h,'XTick',[],'XTickLabel',[])
%% Now we compute the spike-triggered average stimulus. This is accomplished
% by taking the 60 milliseconds of stimulus immediately preceding each spike
% and adding them together. This sum is then divided by the total number
% of spikes fired over the course of the entire experiment to determine the
% average stimulus preceding a spike.
% This spike-triggered average is, in a sense, a template for what the neuron
% is "looking for" in the stimulus.
totalCount=sum(spikeCounts)
a=zeros(kernelSize, 1);
for i=kernelSize+1:length(spikeCounts)
a = a + spikeCounts(i)*stimulus((i-kernelSize+1):i);
end
a=a/totalCount;
%% We'll first look at the answer; then unpack what we did
figure(5);
clf;
h=plot(a, 'o-');
set(h, 'linewidth', 2);
set(gca, 'fontsize', 14);
xlabel('time'); ylabel('estimated linear filter')
%% Because this is a tutorial, we *know* exactly what filtering
% was done on the stimulus to get the linear response ("linearResp").
% Below, we compare the spike-triggered average to the filter we used.
% They are similar shape up to a constant multiplication factor
figure(6);
clf;
h=plot([a/abs(sum(a)),filter/abs(sum(filter))], 'o-');
set(h, 'linewidth', 2);
set(gca, 'fontsize', 14);
legend('estimated filter', 'actual filter', 'Location', 'NorthWest')
%% Remember we summed together stimuli that led to a spike
% We can look at individual such stimuli and the average
% as we have more samples. This code is slower, so we'll
% do steps of 100
atmp = 0; totalTmp = 0;
figure(7)
thelen = min(numsamples, 10000);
for i=kernelSize+1:thelen
atmp = atmp + spikeCounts(i)*stimulus((i-kernelSize+1):i);
totalTmp = totalTmp+spikeCounts(i);
if spikeCounts(i)==1 & mod(totalTmp, 100)==0
subplot(2,1,1);
h=plot(stimulus((i-kernelSize+1):i));
set(gca, 'fontsize', 14); set(h, 'linewidth', 2);
title('stimulus that resulted in spike')
atmp = atmp + stimulus((i-kernelSize+1):i);
subplot(2,1,2);
h=plot(atmp); set(h, 'linewidth', 2);
set(gca, 'fontsize', 14);
title('average stimulus that resulted in spike')
pause(.5);
end
end
%% Extra intuition
% To get intuition about why averaging the spiked stimuli
% works, we can look at how the linear response relates to
% spikes versus no spikes. This allows us to differentiate
% between stimuli that lead to spikes or no spikes.
figure(8)
linearEst=zeros(size(spikeCounts));
for i=kernelSize+1:length(spikeCounts)
linearEst(i) = a' * stimulus((i-kernelSize+1):i);
end
plot(linearEst,spikeCounts,'o');
set(gca, 'fontsize', 12);
xlabel('Linear response');
ylabel('Spikes');
title('Relation between estimated linear responses and spikes')
%% Extra
%% We can also actually estimate the nonlinearity.
% We can plot the *average* number of spikes fired in response
% to similar linear responses. First we decide on linear
% response bins...
figure(9);
x=[-.2:.05:.3];
% ...and then calculate the averages.
N=zeros(size(x));
se=zeros(size(x));
for i=1:length(x)
L = logical((linearEst>x(i)) & (linearEst<(x(i)+.05)));
N(i)=mean(spikeCounts(L));
se(i)=sqrt(var(spikeCounts((linearEst>x(i)) & (linearEst<(x(i)+.05))))/sum(L));
end
x=x+.025;
errorbar(x,N,se,'ro')
set(gca,'xLim',[min(x), max(x)]);
xlabel('Linear response component');
ylabel('Mean spike count');
% We can compare this to the nonlinearity that we know because we constructed
% the simulation - but usually would not know.
% That is, we can superimpose the non-linear function that we actually used to
% determine the spike firing probabilities. The plot of linear response versus
% mean spike count should have the same shape as this function, but, there
% is an arbitrary scale factor relating these two quantities. Below, we estimate
% this scale factor using least-squares.
xx = getNonlinear1(x); % this we know
gamma = inv(xx*xx') * xx * N'; % find the scale factor
vals=[-.4:.05:.4];
Nth = getNonlinear1(vals);
hold on
plot(vals,Nth*gamma)
hold off
%% Questions
% Each time make only one change to the code:
%
% 1. Try lowering the numsamples in the code above
% (example, to 200, 700, and 2000)? Plot the
% filter estimates for these. When does the filter
% estimate break down? Why does it break down?
%
% 2. Try changing the linear function (choosing between getLinear1,
% getLinear2, getLinear3; see toggle comment). Can we recover each of
% the linear filters properly?
%
% 3. Try changing the nonlinear function (choose between
% getNonlinear1, getNonlinear2.m) and keeping getLinear1.
% Can we recover the linear filter of the neuron for each
% of these nonlinearities? If not, then why?
%
% Extra:
% 4. What is special about the stimuli that lead to a spike?
% What would happen if we tried to average random stimuli
% without knowing which ones lead to a spike? To do so, you can
% try to create a false spike train each time (for instance do help on
% randperm, and try to permute spikeCounts) and compute the spike-triggered
% average.