Goals:

• Appreciate the conceptual importance of relating LFP patterns to spiking activity
• Implement several measures of spike-field relationships
• Understand the relative strengths and weaknesses of different measures
• Cultivate awareness of some common pitfalls

Resources:

• (for reference) FieldTrip spikefield tutorial, on which this module is based
• (for reference) Vinck et al. 2012 "Improved measures of phase-coupling between spikes and the local fields potential"

Assignment:

• No assignment this week, but remember to upload your sandbox file, ideally with your replies to questions and action items (☛), so I can see your progress!

Local field potentials and related quantities such as the EEG ultimately result from the spatiotemporal summation of electrical currents across the cell membrane of neurons. Thus, in general, it is expected that spiking activity should be related to LFPs somehow. However, the exact nature of this relationship has proven to be far from trivial, and depends on many factors such as the three-dimensional arrangement of neurons, ion channel distributions, and more rapid dynamics arising from the interaction of multiple inputs. Accordingly, measures that characterize the relationship between spikes and LFPs have painted a surprisingly rich picture of how individual neurons participate in population- and systems-level activity (see e.g. Womelsdorf et al. 2005, Benchenane et al. 2010 and many others for striking examples).

A different, practical issue is that LFPs are susceptible to volume conduction, that is, they can be recorded some distance away from their source. The exact amount distance depends on numerous factors, but even intracranially in the rat the hippocampally-generated theta rhythm can be recorded several millimeters away in the cortex (Sirota et al. 2008). Thus, is is not a priori clear that a LFP recorded from a particular brain structure is generated there (is locally relevant) raising the possibility of attributing properties of the recording to the wrong site! One of the major ways to determine if a LFP is locally relevant is to establish if it is related to spiking activity, which as a fast and relatively local signal does not volume-conduct nearly as far.

In this module, we will examine basic measures of spike-field relationships useful for both approaches.

Consider the following snippet of a LFP and spike trains from three neurons, simultaneously recorded from dorsal CA1 in the rat hippocampus:

From visual inspection alone, it seems clear that the spike times have some relationship to the LFP. The spikes in the top row, for instance, appear to occur preferentially at the troughs of the prominent theta oscillation (count the peaks in one second, there are approximately 8). The neuron on the middle row seems to fire at a wider ranges of spikes, but appears to be avoiding the throughs, as is also apparent in the antiphase (alternating) relationship with the top neuron. (If you want to recreate this plot: the LFP is R016-2012-10-03-CSC02d.ncs and the three neurons are the .t files for the same session.)

How can we characterize and quantify such relationships?

Perhaps the simplest way to visualize possible relationships between spikes and a LFP is to simply average, over all spikes, the surrounding LFP. This is known as the spike-triggered average (STA) and is easy to obtain. Before loading the data, rename any *._t files to *.t; this is because FieldTrip will reject the files otherwise. Then, load the data:

% cd to R016-2012-10-03 folder
 
sd.fc = FindFiles('*.t');
sd.S = LoadSpikes(sd.fc);
 
csc = myLoadCSC('R016-2012-10-03-CSC02d.ncs');
 
cscR = Range(csc); cscD = Data(csc);
 
Fs = 2000; dt = 1./Fs;
cscD = locdetrend(cscD,Fs,[1 0.5]); % remove slow drifts in signal (this can mess up the STA)

Then, we simply grab a LFP snippet for each spike and average:

w = [-1 1]; % time window to compute STA over
tvec = w(1):dt:w(2); % time axis for STA
 
iC = 3; % only do the third neuron for now
clear sta;
 
spk_t = Data(sd.S{iC});
 
h = waitbar(0,sprintf('Cell %d/%d...',iC,length(sd.S)));
 
for iSpk = length(spk_t):-1:1 % for each spike...
 
   sta_t = spk_t(iSpk)+w(1);
   sta_idx = nearest(cscR,sta_t); % find index of leading window edge
 
   toAdd = cscD(sta_idx:sta_idx+length(tvec)-1); % grab LFP snippet for this window
   % note this way can be dangerous if there are gaps in the data
 
   sta{iC}(iSpk,:) = toAdd'; % build up matrix of [spikes x samples] to average later
 
   waitbar(iSpk/length(spk_t));
end
 
close(h);

If you run this, you will see it is quite slow, as indicated by the waitbar. When it finishes however, you can do

plot(tvec,nanmean(sta{3}),'k','LineWidth',2); 
set(gca,'FontSize',14,'XLim',[-0.5 0.5]); xlabel('time (s)'); grid on;

You should get:

Some features of interest of this STA include:

• There is a sharp spike visible at time 0. On the one hand, this is as expected; after all, time 0 is defined as the time of a spike. However, recall that the LFP is filtered between 1 and 425Hz, as can be verified by inspecting the header (getHeader(csc)). Clearly, this filter does not eliminate all spike components.
• There is a clear oscillation in the theta range visible. This indicates that the spikes tend to occur firstly, during theta oscillations, and secondly, at a non-random phase of this oscillation (otherwise the oscillations would average out, assuming there are sufficient spikes that make up the average)

☛ There seems to be a certain asymmetry in the STA theta oscillation (clearer on the left than on the right). How might we interpret this?

In situations like this, where your code is running slowly, it is often useful to use the profiler, a built-in tool that monitors the CPU time spent in each command. We may suspect that in the above code, the loop over spikes and the many calls to nearest() that result are the cause of slowness, but we can confirm this with a test. Enable the profiler as follows:

profile on

Now run the STA code again, but now for cell 1 rather than cell 3. Wait until it completes. You can now view the profiler report with the profview command. As suspected, we see that nearest() takes up the bulk of computation time. Turn the profiler off again with profile off.

A faster way to compute the STA, that avoids the for loop and the use of nearest() is the following:

bin_edges = cscR+(dt/2);
len = length(tvec);
 
clear sta;
iC = 1;
 
spk_ts = Restrict(sd.S{iC},cscR(1)-w(1),cscR(end)-w(2));
spk_t = Data(spk_ts)+w(1); % times corresponding to start of window
 
[~,spk_bins] = histc(spk_t,bin_edges); % index into data for start of window
 
spk_bins2 = spk_bins(:, ones(len,1));
toadd = repmat(0:len-1,[length(spk_bins) 1]);
 
spk_bins3 = spk_bins2+toadd;
 
sta = cscD(spk_bins3);
 

This solution is about a hundred times faster, a speed-up accomplished by constructing a large matrix of appropriately chosen indices into the LFP: each row contains sequential indices corresponding to a LFP snippet centered on a single spike.

☛ Plot the STA by averaging the sta matrix over spikes. What differences between the STA for neurons 1 and 3 do you notice?

For a more sophisticated analysis of STAs, we move to FieldTrip. Since we now need to load not just LFPs, but also spikes, there are a few preliminaries to take care of first.

☛ Grab an updated version of read_mclust_t() to replace the version you currently have.

Now we can load some LFPs and the spikes from this neuron. Ignore the usual warnings about timestamp mismatches arising from gaps in the LFP data.

spike = ft_read_spike('R016-2012-10-03-TT02_2.t'); % needs fixed read_mclust_t.m
fc = {'R016-2012-10-03-CSC02b.ncs','R016-2012-10-03-CSC02d.ncs','R016-2012-10-03-CSC03d.ncs'};
data = ft_read_neuralynx_interp(fc);
data_all = ft_appendspike([],data, spike);

ft_appendspike() creates a binned version of the spike train, with bins set to match the sampling frequency of the LFP. You can verify this as follows:

plot(data_all.time{1},data_all.trial{1}(1,:)) % a LFP
hold on;
plot(data_all.time{1},data_all.trial{1}(4,:)*500,'r') % binarized spike train (0 means no spike, 1 means spike)

data_all.labels keeps track of what is what, with the indices matching those of data_all.trial{1}. Recall that currently there is only one cell because by default, FieldTrip treats the entire recording as one “trial”.

Clearly, some precision about spike times is lost when converting spike times to this binarized format: as a result, spike times are now only known with a precision of 1/2000 s, rather than Neuralynx's native 0.1 us. However, when examining spike-field relationships for oscillations below 100Hz or so, 1/2000 is generally sufficient. FieldTrip has ways of maintaining full-precision spike representations explained in the spikefield tutorial, but for now we will use this fast and convenient method.

Let's compute and plot the STA for this neuron:

cfg              = [];
cfg.timwin       = [-0.5 0.5]; %
cfg.spikechannel = spike.label{1}; % first unit
cfg.channel      = data.label(1:3); % first 3 LFPs
staAll           = ft_spiketriggeredaverage(cfg, data_all);
 
% plot
figure
plot(staAll.time, staAll.avg(:,:)');
legend(data.label); h = title(cfg.spikechannel); set(h,'Interpreter','none');
set(gca,'FontSize',14,'XLim',cfg.timwin,'XTick',cfg.timwin(1):0.1:cfg.timwin(2)); 
xlabel('time (s)'); grid on;

You should see:

FieldTrip makes it easy to process several LFPs simultaneously. Note that as in our manually computed STA above, the LFP channels on the same tetrode as the spikes (“CSC02a-d” are the LFP channels for tetrode “TT02”) have the clear spike artifact at time zero. Also, the red STA is from a different tetrode (“TT03”) and consequently looks very different. In fact, this tetrode was positioned not in the hippocampus (as TT02 was) but in the ventral striatum. If you look carefully, you can just make out a faint theta modulation in the vStr STA as well.

One common question is to ask if the spike-LFP relationship is related to task events. To explore this, we first need to segment our data into trials, as we have done in previous modules:

cfg = [];
cfg.trialfun = 'ft_trialfun_lineartracktone2';
cfg.trialdef.hdr = data.hdr;
cfg.trialdef.pre = 2.5;
cfg.trialdef.post = 5;
 
cfg.trialdef.eventtype = 'nosepoke'; % could be 'nosepoke', 'reward', 'cue'
cfg.trialdef.location = 'both'; % could be 'left', 'right', 'both'
cfg.trialdef.block = 'both'; % could be 'value', 'risk'
cfg.trialdef.cue = {'c1','c3','c5'}; % cell array with choice of elements {'c1','c3','c5','lo','hi'}
 
[trl, event] = ft_trialfun_lineartracktone2(cfg);
cfg.trl = trl;
 
data_trl = ft_redefinetrial(cfg,data_all);

Notice that the contents of data_trl now have multiple cells (88) corresponding to trials, rather than just one for the entire session as before.

Now, we can compute the STA for specific task segments, defined relative to task event time zero, the nosepokes into the reward receptacle:

cfg              = [];
cfg.timwin       = [-0.5 0.5];
cfg.latency      = [-2.5 0];
cfg.spikechannel = spike.label{1}; % first unit
cfg.channel      = data.label(1:3); % first 3 LFPs
staPre           = ft_spiketriggeredaverage(cfg, data_trl);
 
% plot
figure
plot(staPre.time, staPre.avg(:,:)');
legend(data.label); h = title(cfg.spikechannel); set(h,'Interpreter','none');
set(gca,'FontSize',14,'XLim',cfg.timwin,'XTick',cfg.timwin(1):0.1:cfg.timwin(2)); 
xlabel('time (s)'); grid on;

Note the use of the latency field: this restricts the STA to the “reward approach” from 2.5 seconds before to the nosepoke.

☛ How does this pre-nosepoke STA compare to the non-restricted STA? Confirm your impressions by also computing and plotting a “post-nosepoke” STA of equivalent window length.

The spike artifact at time zero is going to cause us trouble later on. FieldTrip provides functionality for removing it:

cfg              = [];
cfg.timwin       = [-0.002 0.006]; % remove 4 ms around every spike
cfg.spikechannel = spike.label{1};
cfg.channel      = data.label(2); % only remove spike in the second LFP ('02d')
cfg.method       = 'linear'; % remove the replaced segment with interpolation
data_trli        = ft_spiketriggeredinterpolation(cfg, data_trl);

If you now recompute and plot the pre-nosepoke STA as above (remembering to use data_trli, not the original data_trl) you should get:

The spike artifact is now no longer visible in the green '02d' STA.

Note that a possibly better way, avoiding potential issues arising from the removal procedure, is to compute the STA against a nearby channel that does not have the spike waveform; however, this is not always possible.

The STA gives a useful visualization of spike-field relationships, but it is not well suited for statistical tests, or visualizations of changes over time. To do this, we start with computing the Fourier transform of a small window of data centered on every spike, to obtain a periodogram with magnitude and phase information. This idea can be visualized as follows:

(source: FieldTrip spike-field tutorial; “DFT” is a specific implementation of the Fourier transform for sampled data)

The code:

cfg              = [];
cfg.method       = 'convol';
cfg.foi          = 1:1:100;
cfg.t_ftimwin    = 5./cfg.foi; % 5 cycles per frequency
cfg.taper        = 'hanning';
cfg.spikechannel = spike.label{1};
cfg.channel      = data.label{1};
stsConvol        = ft_spiketriggeredspectrum(cfg, data_trli);

To understand the output struct stsConvol better, we can plot the spike-triggered PSD by averaging over spikes:

plot(stsConvol.freq,nanmean(sq(abs(stsConvol.fourierspctrm{1}))))

You will see a clear peak in the theta range, indicating that spikes are associated with elevated theta power.

Similarly, we can plot the distribution of phases in the theta frequency range:

hist(angle(stsConvol.fourierspctrm{1}(:,:,9)),-pi:pi/18:pi)

The spike phases are not distributed uniformly across the phase range, indicating a relationship between the spike times and the LFP. Indeed, the spike phases form the basis for statistical measures of spike-field relationships, such as the pairwise phase consistency (PPC; Vinck et al. 2010):

cfg               = [];
cfg.method        = 'ppc0'; % compute the Pairwise Phase Consistency
cfg.spikechannel  = stsConvol.label;
cfg.channel       = stsConvol.lfplabel; % selected LFP channels
cfg.avgoverchan   = 'unweighted'; % weight spike-LFP phases irrespective of LFP power
cfg.timwin        = 'all'; % compute over all available spikes in the window
cfg.latency       = [-2.5 0];
statSts           = ft_spiketriggeredspectrum_stat(cfg,stsConvol);

% plot the results
figure
plot(statSts.freq,statSts.ppc0')
xlabel('frequency')
ylabel('PPC')

You should get:

A clear peak in the theta frequency is visible, suggesting a spike-field relationship for that range. However, there is a strange increasing value with increasing frequency. Could this be due to the spike artifact at time zero?

☛ Modify the cfg for ft_spiketriggeredspectrum() and subsequent steps, but now using the LFP channel from which we had previously removed the spike artifact. What does the the resulting PPC over frequency look like?

The pairwise phase consistency is a modified version of the phase-locked value, the length of the vector mean of all the phases. PLV is intuitive, but suffers from an upward bias for low numbers of spikes, because if there is only one spike, the resultant mean vector length is also 1.

We can compute the PPC using a sliding window:

param = 'ppc0'; % set the desired parameter
 
cfg                = [];
cfg.method         = param;
 
cfg.spikechannel  = stsConvol.label;
cfg.channel       = stsConvol.lfplabel; % selected LFP channels
cfg.avgoverchan    = 'unweighted';
cfg.winstepsize    = 0.01; % step size of the window that we slide over time
cfg.timwin         = 0.5; % duration of sliding window
statSts            = ft_spiketriggeredspectrum_stat(cfg,stsConvol);
 
figure
cfg            = [];
cfg.parameter  = param;
cfg.refchannel = statSts.labelcmb{1,1};
cfg.channel    = statSts.labelcmb{1,2};
cfg.xlim       = [-2 3]; cfg.ylim = [2 30];
ft_singleplotTFR(cfg, statSts)

This results in:

This plot indicates that there is some phase-locking in the high theta range, which disappears from -0.5 seconds to about 0.5 seconds after the nosepoke, after which a complex pattern appears. Some of this may be due to chewing artifacts.

An alternative method of quantifying spike-field relationships is to not extract a phase for each spike, as done above, but instead compute a spectrum over multiple spikes. This uses the same approach as what we did previously for LFP-LFP coherence, relying on the idea of a cross-spectral density:

cfg              = [];
cfg.output       = 'powandcsd';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.foi          = 1:1:100; % frequencies to use
cfg.t_ftimwin    = 5./cfg.foi;  % frequency-dependent, 5 cycles per time window
cfg.keeptrials   = 'yes';
cfg.channel      = {'R016-2012-10-03-CSC02b', 'R016-2012-10-03-TT02_2'};
cfg.channelcmb   = {'R016-2012-10-03-CSC02b', 'R016-2012-10-03-TT02_2'}; % channel pairs to compute csd for
 
cfg.toi          = -2:0.05:3;
 
TFR_pre = ft_freqanalysis(cfg, data_trl);
 
cfg            = [];
cfg.method     = 'ppc'; % compute coherence; other measures of connectivity are also available
fd             = ft_connectivityanalysis(cfg,TFR_pre);
 
iC = 1; % which signal pair to plot
lbl = [fd.labelcmb{1,:}]; % get the label of this pair
imagesc(fd.time,fd.freq,sq(fd.ppcspctrm(iC,:,:))); axis xy; colorbar
xlabel('time (s)'); ylabel('Frequency (Hz)'); title(lbl);

This gives:

As noted previously, the Fourier transform of the autocorrelation function of a LFP is related to the Fourier transform of the signal itself, and the Fourier transform of the cross-correlation function gives the cross-spectral density, the basis for the coherence measure that characterizes the degree to which two signals maintain a consistent phase relationship. Since we can compute the acf and ccf of spike trains, this means we can similarly compute a cross-spectral density between a spike train and a LFP. However, the spike spectrum does not maintain information about how many spikes were used to estimate it, and as a result the spike-field coherence is also subject to biases. The PPC method based on a separate FFT for each spike is the preferred current method.

The idea of assigning a phase to each spike also lends itself to other interesting analyses. We can for instance ask if spike phase varies systematically as a function of task:

csc = myLoadCSC('R016-2012-10-03-CSC02d.ncs');
cscR = Range(csc); cscD = Data(csc);
 
% filter in theta range
Fs = 2000;
Wp = [ 6 10] * 2 / Fs;
Ws = [ 4 12] * 2 / Fs;
[N,Wn] = cheb1ord( Wp, Ws, 3, 20); % determine filter parameters
[b_c1,a_c1] = cheby1(N,0.5,Wn); % builds filter
 
csc_filtered = filtfilt(b_c1,a_c1,cscD);
phi = angle(hilbert(csc_filtered));

Now we have a time series of theta phases, which we can use to get a theta phase for each spike:

S = LoadSpikes({'R016-2012-10-03-TT02_1.t'});
spk_t = Data(S{1});
spk_phi = interp1(cscR,phi,spk_t,'nearest');
 
hist(spk_phi,-pi:pi/18:pi)

As before, the phase histogram is non-uniform, indicating a spike-LFP relationship. Now we can plot phase as a function of position:

[Timestamps, X, Y, Angles, Targets, Points, Header] = Nlx2MatVT('VT1.nvt', [1 1 1 1 1 1], 1, 1, []);
Timestamps = Timestamps*10^-6;
toRemove = (X == 0 & Y == 0);
X = X(~toRemove); Y = Y(~toRemove); Timestamps = Timestamps(~toRemove);
 
spk_x = interp1(Timestamps,X,spk_t,'linear');
spk_y = interp1(Timestamps,Y,spk_t,'linear');
 
plot(X,Y,'.','Color',[0.5 0.5 0.5],'MarkerSize',1); axis off; hold on;
h = scatterplotC(spk_x,spk_y,spk_phi,'Scale',[-pi pi],'solid_face',1,'plotchar','.');

You should get:

Each spike is now plotted according to where the rat was when it occurred, and the color indicates the theta phase. Note that there is a systematic phase change visible as the rat runs along the bottom edge of the track: this is the famous theta phase precession phenomenon characteristic of hippocampal place cells.