**Reference**

HowToGuides

Manuals

LabAlumni

DataAnalysis

Advice for...

**Admin**

analysis:nsb2015:week13

~~DISCUSSION~~

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"

The code in this module works, but needs to be converted to use the most recent data types

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.

For this section to work, you need to use the version of `read_mclust_t.m`

that is in the `nsb2015`

codebase, **NOT** the one in fieldtrip!

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.

analysis/nsb2015/week13.txt · Last modified: 2018/07/07 10:19 (external edit)

Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Share Alike 4.0 International