Spike tuning curve

image

Representing the neuron's preference.

Background

While neuronal activity is inherently noisy, neurons still prefer particular values of certain features. For example, a neuron from the primary auditory cortex (A1) will respond vigorously to pure tones of a particular frequency, while the activity drops for tones of increasingly different frequencies. A "visual" neuron might respond preferentially for a particular spatial orientation of a visual bar, or the velocity of a moving visual grating.

Data

Here, you will plot the tuning curve of an A1 neuron in the auditory cortex. The data is from the same neuron as the raster plot tutorial: neuronal activity (spike timing) was recorded while 13 tones of different frequencies (200 Hz to 16 kHz, logarithmically spaced) at 4 different sound levels (40, 50, 60 and 70 dB) were presented 10 times (making a total of 520 trials). The spike timings and the stimulus parameters are stored in a structure. Spike timings are stored in the field spiketime, and the stimulus values in stimvalues.

Analysis

To obtain the stimulus values of the first trial, type:

 
load('tonespike'); % load data: structure called tonespike
stim = tonespike(1).stimvalues; % stimulus values

stim is a 4x1 vector. To which stimulus features do these 4 stimulus values relate? This is conveniently stored in the field stimparams:

 
params = tonespike(1).stimparams; % stimulus feature
 

This is a cell-array of 4 strings (pieces of texts): the first stimulus feature is frequency ('Freq [Hz]'), the second is sound level ('SPL [dB]'), the third sound duration ('Duration [ms]'), and the fourth speaker channel ('Channel'). For the current experiment, these features are the same for every trial. Here, we will create a frequency tuning curve, and so we are only interested in the first feature.

 
feature = stim(1); % frequency (Hz)

This first trial contains a tone of 250 Hz. How many spikes does the neuron produce during the presentation of this tone? For that we need to select the spikes occurring during the stimulus presentation:

 
spiketime   = tonespike(1).spiketime;
% stimulus was presented between 300 and 450 ms
onset     = 300; % tone onset (ms)
offset   = 450; % tone offset (ms)
sel    = spiketime>=onset &amp; spiketime<=offset; % selection vector!
Nspikes   = sum(sel); % number of spikes during tone presentation
 

The neuron spikes only twice (Nspikes) in this trial! As we can know from the raster plot of the data, the neuron's firing pattern varies idiosyncratically across trials. We are therefore usually not interested in single-trial spike numbers, but want to know how many spikes the neuron produces on average during all trials in which a 250-Hz tone is produced. To obtain this average, we first need to determine the tone frequencies that were presented in every trial. We need to do this using a for-loop:

 
ntrials = numel(tonespike); % number of trials
Feature = NaN(ntrials,1); % always initialize vector/matrix!
for ii = 1:ntrials
  % params = tonespike(ii).stimparams; % the stimulus parameter names of trial ii
  values = tonespike(ii).stimvalues; % the stimulus parameter values of trial ii
  % From params we know that the first entry of the values-vector
  % contains the sound's frequency 
  Feature(ii) = values(1); % (Hz)
end

Note that we initialized a vector Feature before the for-loop by creating a vector containing only NaNs (Not-a-Numbers). In the loop this vector will be gradually (trial-by-trial) filled with the actual tone frequencies. Including this step is essential to reduce the time it takes for Matlab. It is also good practice: thinking about how large a vector or matrix should be increases your understanding of the data, and reduces programming errors.

Now we can select the trials in the tonespike structure with the 250 Hz tone presentation.

 
sel      = Feature == 250; % select all trials with tone frequency uFeature(ii)
spiketime    = [tonespike(sel).spiketime];
nrep      = sum(sel); % number of stimulus repeats
 
onset      = 300; % tone onset (ms)
offset      = 450; % tone offset(ms)
sel      = spiketime>=onset & spiketime<=offset; % stimulus was presented between 300 and 450 ms
Nspikes    = sum(sel); % number of spikes
 

Now we're getting somewhere: 170 spikes! Ofcourse, this number is pretty useless, as it depends on the duration of the tone, the number of stimulus repeats and the sampling rate. A much more interesting number is the average firing rate of the neuron:

 
duration    = offset-onset; % ms
Fs        = 1000; % sampling rate (Hz)
firingRate    = Nspikes/duration*Fs/nrep; % in spikes/s or Hz
 

28 spikes/second! Finally, we need to do this for every tone to be able to determine the frequency tuning curve. Let's use a for-loop, and determine the firing rate for one frequency in every step:

 
uFeature    = unique(Feature); % the unique frequencies presented during the recordings
nFeature    = numel(uFeature); % the number of freqs
firingRate  = NaN(nFeature,1); % Initialization of firingRate vector
Fs      = 1000; % sampling rate (Hz)
onset    = 300; % tone onset (ms)
offset    = 450; % tone offset(ms)
duration  = offset-onset; % ms
for ii = 1:nFeature
  sel        = Feature == uFeature(ii); % select all trials with tone frequency uFeature(ii)
  spiketime    = [tonespike(sel).spiketime];
  nrep      = sum(sel); % number of stimulus repeats
 
  sel        = spiketime>=onset & spiketime<=offset; % stimulus was presented between 300 and 450 ms
  Nspikes      = sum(sel); % number of spikes
  firingRate(ii)  = Nspikes/duration*Fs/nrep; % in spikes/s or Hz
end
 

And plot this:

 
h = semilogx(uFeature,firingRate); % Plot
set(h,'Color','k','LineWidth',2,'LineStyle','-');
set(gca,'XTick',uFeature(1:2:end),'XTickLabel',round(uFeature(1:2:end)));
box off
xlabel('Frequency (Hz)');
ylabel('Firing rate (Hz)');
axis square;
xlim([0.9*min(uFeature) 1.1*max(uFeature)]);
ylim([0 60]);
 

We have plotted the data with semilogx. This is similar to the plot function but the scale of the x-axis is logarithmic. The reason we use this is because auditory neurons respond to tones in a logarithmic fashion.

image 

We can see in the tuning curve that a frequency of 1000 Hz produces on average the highest firing rate of 50 spikes/s. Or we can determine this in Matlab:

 
[mx,indx] = max(firingRate); % highest firing rate (spikes/s)
BF = uFeature(indx); % best frequency (Hz)
 

Matlab Code

Putting everything together:

 
clc;
close all hidden;
clear all hidden;
 
load('tonespike'); % load data: structure called tonespike
 
%% Vectorize 
ntrials = numel(tonespike); % number of trials
Feature = NaN(ntrials,1); % initialize frequency vector
for ii = 1:ntrials
  params = tonespike(ii).stimparams; % the stimulus parameter names of trial ii
  values = tonespike(ii).stimvalues; % the stimulus parameter values of trial ii
  Feature(ii) = values(1); % (Hz)
end
 
uFeature    = unique(Feature); % the unique frequencies presented during the recordings
nFeature    = numel(uFeature); % the number of freqs
firingRate  = NaN(nFeature,1); % Initialization of firingRate vector
firingRateSE = firingRate;
Fs      = 1000; % sampling rate (Hz)
onset    = 300; % tone onset (ms)
offset    = 450; % tone offset(ms)
duration  = offset-onset; % ms
for ii = 1:nFeature
  sel        = Feature == uFeature(ii); % select all trials with tone frequency uFeature(ii)
  spike    = tonespike(sel);
  nrep      = sum(sel); % number of stimulus repeats
  FR = NaN(nrep,1);
  for jj = 1:nrep
    t = spike(jj).spiketime;
    sel        = t>=onset & t<=offset; % stimulus was presented between 300 and 450 ms
    Nspikes      = sum(sel); % number of spikes
    FR(jj) = Nspikes/duration*Fs;
  end
 
  firingRate(ii)    = mean(FR); % in spikes/s or Hz
  firingRateSE(ii)  = std(FR)./sqrt(nrep); % in spikes/s or Hz
 
end
 
%% plot
[hp,hl] = pa_errorpatch(uFeature,firingRate,firingRateSE,'b'); % Plot
set(hl,'Color','k','LineWidth',2,'LineStyle','-');
 
set(hp,'FaceColor',[.7 .7 .7]);
 
set(gca,'XTick',uFeature(1:2:end),'XTickLabel',round(uFeature(1:2:end)),'Xscale','log');
box off
xlabel('Frequency (Hz)');
ylabel('Firing rate (Hz)');
axis square;
xlim([0.9*min(uFeature) 1.1*max(uFeature)]);
ylim([0 60]);
 
%% Save
print('-depsc','-painter',mfilename);
 
%% Best Frequency
[mx,indx] = max(firingRate); % highest firing rate (spikes/s)
BF = uFeature(indx); % best frequency (Hz)
str = ['Best frequency = ' num2str(round(BF)) ' Hz'];
title(str);