Differences

This shows you the differences between two versions of the page.

--- analysis:nsb2014:week8 [2014/07/23 11:35]
mvdm [Spike autocorrelation function (acf)]
+++ analysis:nsb2014:week8 [2022/06/29 19:01] (current)
mvdm
@@ Line 1: / Line 1: @@
-~~DISCUSSION~~
+===== Time-frequency analysis: spectrograms =====
-===== Spike train analysis: firing rate, interspike interval distributions, auto- and cross-correlation =====
 Goals:
-  * Understand the need for statistical measures to describe patterns of spikes
+  * Understand the idea of the spectrogram and the inputs and outputs of typical implementations
-  * Gain familiarity with basic descriptors of spike trains: firing rate, interspike interval distribution, autocorrelation function, cross-correlation function
+  * Apply your knowledge of spectral leakage to the selection of spectrogram parameters
-  * Generate artificial Poisson spike trains and compare these to real data
+  * Cultivate an awareness of common pitfalls in spectrogram-based analysis
-  * Develop intuitions about how Poisson spike trains can be used to test interesting hypotheses about neural coding
+  * Get started with using the fieldtrip toolbox
-Resources:
+Exercise:
-  * (optional background reading, a classic book) [[http://mitpress.mit.edu/books/spikes | Spikes: reading the neural code]]
+  * Construct a event-triggered LFP plot
-  * (optional visual clarification) [[http://www.jhu.edu/signals/discreteconv2/ | convolution example]]
-  * (more mathematical treatment, for reference) [[http://icwww.epfl.ch/~gerstner/SPNM/node34.html | Gerstner & Kistler chapter]] on spike train statistics
-Exercise:
+Resources:
+  * (for reference) [[http://www.chronux.org/tutorials | Chronux toolbox slides]]
+  * (for reference) [[http://fieldtrip.fcdonders.nl/ | Fieldtrip documentation]] and [[http://fieldtrip.fcdonders.nl/tutorial/timefrequencyanalysis | spectrogram tutorial]]
-  * Explore neural coding in the rat hippocampus by comparing cross-correlations to synthetic data
-  *
 ==== Introduction ====
-Action potentials, or //spikes//, are the main unit of information processing in the brain. Spikes enable fast yet reliable electrical communication over more than the microscopic distances to which passive conduction is limited. Although there are exceptions, spikes are generally the only aspect of neural activity that is transmitted rapidly to postsynaptic targets. At the most general level, sensory input to the brain can be thought of as changing the patterns of spikes in associated sensory nuclei and brain regions; in turn, motor output is ultimately governed by patterns of spikes in motor areas. Thus, any adaptive behavior must ultimately result from the appropriate manipulation of spike patterns. As a result, the recording and interpretation of spike data is a central tool in neuroscience.
+The spectral estimation methods introduced in the previous module average spectral content over a number of different windows. This means that we get an improved estimate of the //overall// spectral content of the signal, but such an average cannot establish if the spectral content is changing over time. Neural signals rarely have stationary spectral properties, as illustrated by this snippet of a LFP from the rat ventral striatum (source: [[http://www.frontiersin.org/integrative_neuroscience/10.3389/neuro.07.009.2009/abstract | van der Meer and Redish, 2009]])
-This module introduces the most commonly used basic tools for describing spike trains. We focus first on describing spike trains from a single neuron, with a first step towards pairwise analysis at the end.
+{{ http://c431376.r76.cf2.rackcdn.com/544/fnint-03-009/image_m/fnint-03-009-g004.jpg?600 |}}
-==== Estimating firing rate ====
+The top trace illustrates some striking properties of the ventral striatal LFP: under certain conditions, gamma oscillations appear in "events" of a relatively stereotyped length, with frequencies alternating between ~80 Hz ("gamma-80" or "high gamma") and ~50 Hz ("gamma-50" or "low gamma").
-For most purposes, a spike can be considered as a unitary, punctate event that occurs at a specific time. Of course, in reality, a spike actually has a certain brief lifetime, from its generation and propagation to ultimate decay, but this is not essential to most experimental settings. Thus, a spike train can be described simply as a list of times, each indicating the time at which a given neuron emitted a spike. In statistics, this type of data is known as a //point process//: a spike can be considered a point in time.
+How can we characterize the spectral properties of such a rapidly changing signal? One approach is illustrated in the colorful panel directly below the raw LFP trace -- this is a **spectrogram**. A spectrogram estimates the spectral content of a signal over time. It has time and frequency axes, with power shown in the color dimension (note the colorbar: in this plot, red is high power, blue is low power).
-Because a point is -- at least in principle -- infinitesimally small, no two spikes technically occur at **exactly** the same time. Yet, as the figure below shows, multiple presentations of the same stimulus (panel A) can evoke spike trains that are clearly similar across trials (panel B):
+The spectrogram is constructed by sliding a window along the signal and computing the power spectrum at each position of the window. As with the PSD estimators in the previous module such as Welch's method, the window is characterized by its size and shape (rectangular, Hamming, etc.) which control the tradeoffs resulting from spectral leakage. Unlike these PSD estimators however, to obtain a spectrogram we do not average across windows. Instead, we keep the spectrum for each position of the window, filling in each "time" column of the spectrogram as it moves along.
-{{ :analysis:course:week9_fig2.png?300 |}}
+(As an aside, spectrograms are a widely used tool for time-frequency analysis; for instance, they are often the first step in voice recognition software. Plotting a spectrogram of your own voice can be done with just a few lines of MATLAB code.)
-(from [[http://robby.caltech.edu/~zoran/Reading/bialek97.pdf|de Ruyter van Steveninck et al. 1997]])
+==== Constructing a basic spectrogram ====
-We would like a way to describe the observation that this neuron's response is similar across trials, even though the exact spike times differ between trials.
+=== Understanding the spectrogram() function ===
-=== Spike counts (binned firing rate) ===
+MATLAB's signal processing toolbox includes basic functionality for constructing spectrograms.
-A simple idea is to count how many spikes occur in a given time interval (time bin). Then, the precise spike time is no longer relevant. Let's first plot a short snippet of an example spike train:
+☛ Look at the documentation for the ''spectrogram()'' function. Notice that the input arguments it takes are familiar to you from the last module: window, overlap, nfft, Fs. The usage we are looking for is the last one:
 <code matlab>
-%% loading
+[S,F,T,P] = spectrogram(...)
-cd('C:\data\R042-2013-08-18_recording');
+   P is a matrix representing the Power Spectral Density (PSD) of each segment.
+</code>
-cfg = [];
+Let's give it a spin:
-S = LoadSpikes(cfg);
-%% plot
+<code matlab>
-iC = 47;
+%% load the data
-t = [5801 5801.7];
+% remember to cd to the correct folder here, may need to get this file from the lab database
+cfg = []; cfg.fc = {'R016-2012-10-03-CSC04a.Ncs'};
-Sr = restrict(S,t(1),t(2)); % get spike times
+csc = LoadCSC(cfg);
-plot([Sr.t{iC} Sr.t{iC}],[-1 -0.5],'Color',[0 0 0]); % note, plots all spikes in one command
+Fs = csc.cfg.hdr{1}.SamplingFrequency;
+%% restrict data
+cscR = restrict(csc,3282,3286); % if you don't have this, implement it (it's one line of code!)
+plot(cscR.tvec,cscR.data); % note there are various gamma oscillations present, as well as a large negative-going transient
+%% construct and plot the spectrogram
+[S,F,T,P] = spectrogram(cscR.data,hanning(512),256,1:200,Fs);
+imagesc(T,F,10*log10(P)); % converting to dB as usual
+set(gca,'FontSize',20);
+axis xy; xlabel('time (s)'); ylabel('Frequency (Hz)');
 </code>
-To bin these spikes, we can use MATLAB's ''histc()'' function. This returns a histogram of spike counts, given a vector of bin edges:
+You should see:
+{{ :analysis:course:week6_fig1.png?600 |}}
+As in the example, color indicates signal power at a particular point (bin), corresponding to a specific frequency  and a specific time. Some power is apparent at low frequencies (<10Hz), as well as in the gamma range (~50-100Hz).
+Let's deconstruct this plot a little further. Since we specified 512-sample windows with a 256-sample overlap, we expect roughly 8 estimates per second. This is because our Fs is 2000, and therefore 8 "steps" of our moving window, moving at 256 samples per step (the window size minus the overlap equals the step size) fit into one second. Indeed we can discern approximately 8 "columns" (time bins) per second in the plot. For the frequency axis, we specified ''1:200'' as the bins of interest, so our frequency axis has a resolution of 1 Hz per bin. The centers of the time and frequency bins are returned as the ''T'' and ''F'' outputs respectively.
+An important point about the time axis is this:
 <code matlab>
-binsize = 0.1;
+>> T(1)
-tbin_edges = t(1):binsize:t(2); % vector of time bin edges (for histogram)
-tbin_centers = tbin_edges(1:end-1)+binsize/2; % vector of time bin centers (for plotting)
-spk_count = histc(Sr.t{iC},tbin_edges); % get spike counts for each bin
+ans =
-spk_count = spk_count(1:end-1); % ignore spikes falling exactly on edge of last bin.
+.1280
 </code>
-When using ''histc'', it is key to remember that element //n// of the output corresponds to the count in the interval [edge(//n//) edge(//n+1//)>. In other words, a spike occurring exactly at the time of the second edge will be assigned to output bin 2, not bin 1. As a result, ''histc()'' needs a special last bin for those spikes occurring exactly at the time of the last edge. The above code ignores this last bin, so that we have a number of output spike counts that is equal to the number of bin centers, which in turn is equal to the number of bin edges minus 1.
+That is, the time corresponding to the first column (time bin) of the spectrogram is //not// 0. It is also not the time of the first sample in our data, as is readily verified:
-We can now plot the spike counts as follows:
+<code matlab>
+>> cscR.tvec(1)
+ans =
+.2820e+03
+</code>
+In fact, the time of the first bin returned by ''spectrogram()'' is the //center of the first window for which there is enough data//. The data passed to ''spectrogram()'' is assumed to start at time 0. This makes sense when we remember that the input arguments do not contain timestamps (only voltage samples) so this function has no way of knowing the "experiment time" (''csc.tvec'').
+Thus, the time of the first spectrogram bin (''0.128'') is because this is the center of the first 512-sample window for which there is enough data (''256/2000'' is ''0.128''), and the first sample of the data is assigned time 0.
+=== Aligning the spectrogram with a LFP ===
+Before we start investigating the effects of changing the spectrogram parameters, it is helpful to be able to see the raw LFP overlaid onto the spectrogram:
 <code matlab>
+%%
 hold on;
-h = bar(tbin_centers,spk_count);
-set(h,'BarWidth',1,'EdgeColor','none','FaceColor',[1 0 0]); % reformat bar appearance
+lfp_minmax = 25; lfp_cent = 125; % range and mean of LFP plotting
+tvec0 = cscR.tvec - cscR.tvec(1); % align LFP with spectrogram
-yt = get(gca,'YTick'); ytl = get(gca,'YTickLabel');
+data = rescale(cscR.data,-lfp_minmax,lfp_minmax); data = data+lfp_cent;
-set(gca,'YLim',[-1.5 10],'YTick',yt(2:end),'YTickLabel',ytl(2:end,:)); xlabel('time (s)');
+lfp_h = plot(tvec0,data,'k');
 </code>
-This gives
+Note that to make the LFP and the spectrogram line up, the time of the first LFP sample needs to be zero. The rescaling of the LFP is just to make it easier to see when plotted overlaid onto the spectrogram; of course, if you prefer, you can plot it below or above without overlap by changing ''lfp_cent''.
-{{ :analysis:course:week9_fig1.png?600 |}}
+Now that we can see the LFP and the spectrogram together in one plot, you can observe the correspondence between them. The sharp transient early in the data shows up as a dark red in the low frequency range. This is perhaps not obvious, because such a transient can hardly be described as an oscillation -- it does not even complete a full period (cycle). So beware: if you see strong power in a spectrogram (or a PSD for that matter) this does not imply that there is a true oscillation in the signal.
-☛ Verify by visual inspection that indeed the binned spike counts, shown as red bars, correspond to the number of spikes in each 100ms bin.
+In contrast to this low-frequency transient, you can see that there are various gamma-range oscillations which extend over multiple cycles, and can therefore be regarded as true oscillations rather than transients. For instance, at the start of this data segment there is some low gamma around 50-60Hz, and later some higher frequencies appear.
-However, it should be clear that this count is a poor description of the underlying spike train. For instance, at time 5801.25 the spike count is zero, yet this time does not seem much different from, say, 5801.4. This occurs because we have sharp bin edges that are placed arbitrarily.
+Another thing you will have noticed is that the spectrogram looks pretty coarse. Let's see if we can improve this by changing the input parameters.
-As an aside, note that ''histc()'' provides a raw spike count, which is not normalized by the bin size. Thus, to convert the counts into a firing rate estimate in spikes per second (Hz), we need to divide by the bin size.
+==== Spectrogram parameters ====
-Perhaps we can mitigate the problem of arbitrary edge placement by making the bin size smaller.
+=== Frequency and time resolution ===
-☛ Reduce the bin size to 10 ms and replot.
+In the example above, we used a 512-sample Hanning window with an overlap of 256; that is, the window moves in steps of 256 samples, and as a result we get a new power estimate every 0.128 seconds. In addition, we requested power estimates at integer frequencies from 1 to 200.
-We now have a finer estimate that tracks the spike bursts and gaps more accurately, but it is still dissatisfying to see that some bins have a count of 2 and others of 1, implying that the firing rate for some bins was twice as large as those nearby, without clear justification. Again, how the bin edges happen to align with the spikes is the problem.
+☛ Change the code above so that ''spectrogram()'' uses an overlap of 384 samples (instead of 256) and returns power with 0.25 Hz resolution (instead of 1).
-☛ What happens if you reduce the bin size to 1 ms?
+The spectrogram should now look somewhat smoother. Because of the larger overlap, the spectrogram now contains an estimate every 512-384 = 128 samples, rather than every 256. However, it is important to appreciate that increasing the frequency and time resolution in this way does not change the underlying estimate. This can be seen, for instance, by setting the frequency resolution to something really small (e.g. 0.01). This will take a while to compute, but if you compare the spectrograms side-by-side, it should be clear that the frequency estimates do not become more precise.
-At some point, if you make the spike count bins small enough, there will only ever be zero or one spikes per bin. This defeats the purpose of binning in the first place: now again each spike train will be unique. What we require is a measure that avoids the arbitrariness of bin edges.
+Similarly, you can decrease the overlap to ''window_length - 1'', resulting in a spectrogram that is ultra-smooth in appearance, but this cannot escape the underlying limitation that spectral content is estimated over a window that is 512 samples wide.
-=== A brief diversion: convolution ===
+An important parameter that we know can affect our estimate is the shape of the window used.
-We saw in the previous section that binned spike counts run into problems both when bins are large (location of edges affect the count) and when bins are small (only one spike per bin at most, so no averaging). A way around this is to use small time bins to minimize the edge effects, combined with a certain amount of filtering (smoothing) to enable contributions from multiple spikes to a single bin.
+☛ Change the window used to a rectangular one (''rectwin''). What happens to the spectrogram?
-As an example, imagine that each spike in the above spike train was replaced by a Gaussian, like this:
+Tracking changes over time using a spectrogram generally means that we don't want to use very large time windows. A large window would obscure rapid changes in spectral content, such as the tendency for a low-gamma to follow high-gamma events illustrated in the first figure, above. However, as you know from the previous module, the effects of spectral leakage are more prominent in smaller windows: there is fewer data to make up for the edge effects. This means that spectrograms tend to be particularly sensitive to the choice of window. A Hanning window is often a good default choice. As we will see below, however, there are methods we can use to improve further on this.
-{{ :analysis:course:week9_fig4.png?450 |}}
+=== Frequency and time tradeoffs ===
-Then, we end up with a smoothly varying estimate of firing rate. Formally, this operation of replacing each spike (considered a [[http://en.wikipedia.org/wiki/Dirac_delta_function | Dirac delta function]]) with another function (in this example, a Gaussian) is known as //convolution//. Convolving a spike train is related to the idea of //impulse response// to a filter, which we encountered in an earlier module. Each spike is an impulse (input is zero everywhere except at the time of the spike) and the output is a filtered (convoluted) signal which is a smoothed version of the input.
+The most important choice in constructing spectrograms is the size of the time window. This controls the fundamental tradeoff between our ability to resolve //which frequencies// are present in the signal, and at //what time// those frequencies occur.
-=== Spike density function (SDF, binless firing rate) ===
+This tradeoff arises because of spectral leakage. With a large window, the edge effects are relatively small, and frequency content can be estimated relatively accurately. However, we cannot know at what time in this large window a given oscillatory event occurs. In contrast, using a small window, we may be able to localize more precisely in time when an event occurs. The price we pay is that we have increased spectral leakage, and we are now limited in our ability to determine the frequency of the event.
-To convolve our example spike train to obtain a smooth estimate of firing rate, we can use MATLAB's ''conv2()'' function:
+This idea is illustrated in panel (b) of the figure below:
-<code matlab>
+{{ http://fieldtrip.fcdonders.nl/_media/tutorial/timefrequencyanalysis/tfrtiles.png?600 |}}
-binsize = 0.001; % select a small bin size for good time resolution
-tbin_edges = t(1):binsize:t(2);
-tbin_centers = tbin_edges(1:end-1)+binsize/2;
-spk_count = histc(Sr.t{iC},tbin_edges);
+Rectangular "time-frequency tiles", with different widths and heights, correspond to the area in the spectrogram that would be highlighted by given oscillatory event of a fixed frequency for different choices of time window. For a small window, we can know with more precision when the event occurred, but we are uncertain about the frequency. In contrast, a large window constrains the frequency estimate, but now we are uncertain about time.
-spk_count = spk_count(1:end-1);
-rw_sdf = conv2(spk_count,rectwin(50),'same'); % convolve with rectangular window
+☛ For the piece of data used here, compare side by side (using ''subplot()'') the spectrogram obtained with 512- and 1024-sample windows. For a fair comparison, set the overlap such that you end up with a comparable number of time bins.
-plot(tbin_centers,rw_sdf,'b');
-gau_sdf = conv2(spk_count,gausswin(50),'same'); % convolve with gaussian window
+You should get something like this:
-plot(tbin_centers,gau_sdf,'g');
-</code>
-(The ''same'' option specified in the call to ''conv2()'' ensures that the output is the same size as the input.)
+{{ :analysis:course:week6_fig2.png?600 |}}
-When plotted on top of the earlier binned spike train, you should get
+As expected, the larger window has tighter frequency estimates, but the spectrogram is now more smeared out in time.
-{{ :analysis:course:week9_fig3.png?600 |}}
+These observations raise the question of if there are ways of picking the window size that are somehow optimal. In general, there is no substitute for looking at the raw data and tuning the parameters so that features of interest, such as the gamma events here, are well resolved in the spectrogram. However, there are some sophisticated methods available for which we have to move to the FieldTrip toolbox, which we will get to later in this module.
-Note how convolving with a rectangular function (in blue) results in an estimate that is jagged, but avoids the arbitrary bin edges imposed by the original spike count. For this estimate, each spike is effectively replaced with a rectangle centered at the spike. Convolving with a Gaussian (in green) avoids the jagged edges to make a smooth estimate.
+==== Pitfalls ====
-Convolved spike trains are often referred to as //spike density functions// (SDFs) because they are effectively a continuous estimate of firing rate, i.e. describing the density of spikes over time. Of course, the SDFs we computed here are still discretized, but the resolution could be set to something arbitrarily fine by adjusting the bin size.
+=== Mind the gaps ===
-One final adjustment we need to make is to ensure we actually end up with a firing rate estimate in spikes per second (Hz). Because it is helpful to be able to state the convolution used in a succinct manner, it is best to use the ''gausskernel()'' function, which allows us to specify not only the width of the convolution window, but also the standard deviation of the Gaussian:
+As noted above, the input to ''spectrogram()'' is not a tsd containing samples with accompanying timestamps; rather, it is simply a list of samples assumed to be acquired at a constant sampling frequency ''Fs''. If the input in fact does not have a constant ''Fs'', as can occur for instance if there are gaps in the data, there will be trouble:
 <code matlab>
-binsize = 0.001; % in seconds, so everything else should be seconds too
+cscR = restrict(csc,3300,3340); % if you don't have this, implement it (it's one line of code!)
-gauss_window = 1./binsize; % 1 second window
-gauss_SD = 0.02./binsize; % 0.02 seconds (20ms) SD
+[S,F,T,P] = spectrogram(cscR.data,rectwin(256),128,1:200,Fs);
-gk = gausskernel(gauss_window,gauss_SD); gk = gk./binsize; % normalize by binsize
+imagesc(T,F,10*log10(P)); % converting to dB as usual
-gau_sdf = conv2(spk_count,gk,'same'); % convolve with gaussian window
+set(gca,'FontSize',20);
-plot(tbin_centers,gau_sdf,'g');
+axis xy; xlabel('time (s)'); ylabel('Frequency (Hz)');
+hold on;
+lfp_minmax = 25; lfp_cent = 125; % range and mean of LFP plotting
+tvec0 = cscR.tvec - cscR.tvec(1); % align LFP with spectrogram
+data = rescale(cscR.data,-lfp_minmax,lfp_minmax); data = data+lfp_cent;
+lfp_h = plot(tvec0,data,'k');
+xlim([tvec0(1) tvec0(end)]);
 </code>
-Note that we need to normalize by the binsize (1ms) because our firing rate estimate should be independent of the bin size chosen.
+{{ :analysis:course:week6_fig3.png?600 |}}
-We thus obtain a SDF firing rate estimate with the correct units:
+As you can see, ''spectrogram()'' ignores the gap present in the LFP, causing the spectrogram and the LFP to become misaligned. In general, always verify that your data does not have gaps when computing spectrograms!
-{{ :analysis:course:week9_fig5.png?600 |}}
+=== Sampling frequency mismatches ===
-Of course, there is no single correct choice of the SD of the Gaussian kernel used for the convolution. This should be chosen by comparison with the raw spike train: in this case, we don't want to choose a Gaussian that is so wide that the modulation in firing rate is no longer visible.
+Another way for the spectrogram and LFP to become misaligned is when an incorrect sampling frequency is specified.
-It is also worth noting that convolution with a Gaussian to obtain a firing rate estimate has limitations. For instance, convolving spikes evoked by stimulus presentation might result in an increase in the firing rate estimate before the stimulus is presented. This can be undesirable and should not be interpreted as evidence of "stimulus anticipation" when in fact it simply results from the convolution! To minimize these sorts of issues, a variety of adaptive firing rate estimates have been developed, such as [[http://www.stat.cmu.edu/~kass/bars/bars.html | Bayesian Adaptive Regression Splines (BARS)]].
+☛ Return to the original data segment (''cscR = Restrict(csc,3282,3286);'') and compute the spectrogram using ''Fs-100'' instead of ''Fs''. Plot the spectrogram and overlay the LFP. What do you notice?
-If you compute SDFs frequently, it makes sense to put the above into a function, which takes a ''ts'' object as input (along with binsize, windowsize and SD parameters), and returns a ''tsd'' firing rate estimate.
+==== Loading task events ====
-=== Peri-event or -stimulus time histograms ===
+Computing a spectrogram enables a number of further analyses, such as characterizing the relationship between power at different frequency bands. An example of this is shown by the correlation matrix in panel C of the first figure on this page, which reveals that ventral striatal gamma-50 and gamma-80 tend to be slightly anticorrelated in time.
-PETHs or PSTHs show the average firing rate around a time of interest, usually a task event such as the presentation of a stimulus. Once you have a firing rate estimate for the entire session, it is relatively straightforward to average this estimate over a fixed window around each event.
+Spectrograms can also be used to obtain a time series of power in a certain band over time, similar to the filtering-based approach in the previous module. Such a time series can then for instance be used to determine if power in a certain frequency band is related to a behavioral or task variable such as running speed.
-==== Interspike intervals (ISIs) ====
+Another common spectrogram-based analysis involves aligning the spectrogram to task events of interest, such as a button press, stimulus onset, or receipt of a reinforcer. To do this, we first need to extract the times of these events from the data files.
-A different aspect of single neuron spike trains is the //regularity// of their spike timing. In principle, an average firing rate of say, 10 Hz could come about in a very regular manner (one spike exactly every 0.1s) or a more irregular manner (say, between 0.05 and 0.15s, with an average of 0.1). This idea is illustrated in the following figure, from a [[http://www.jneurosci.org/content/13/1/334.long | famous paper]]:
+=== Processing Neuralynx events files (.nev) ===
-{{ :analysis:course:softkykoch_fig15.png?600 |}}
+Apart from neural data stored in ''.ncs'' and ''.ntt'' files, the Neuralynx system also stores //events// in a file called ''Events.Nev'' (see [[analysis:nsb2014:week2|Module 2]] for details). Here is an example usage for the specific task in this data:
-Panels A-C show simulated voltage traces (the membrane potential of a single model neuron), with the trace in panel A displaying very irregular firing, C very regular firing, and B something in between. How can we capture such differences?
+<code matlab>
+cfg = [];
-=== ISI histograms and the coefficient of variation ===
+cfg.eventList = {'Feeder 0 nosepoke','Feeder 1 nosepoke', ...
+    '1 pellet cue','3 pellet cue','5 pellet cue', ...
+    '1 pellet dispensed','3 pellet dispensed','5 pellet dispensed'};
+cfg.eventLabel = {'n0','n1', ...
+    'c1','c3','c5', ...
+    'd1','d3','d5'};
+evt = LoadEvents(cfg);
+</code>
-A key difference between panels A-C above is the distribution of //interspike intervals// (ISIs). For a very regular spike train, there is a relatively stereotyped time that elapses before the next spike; for an irregular spike train, the interspike intervals are much more variable. Let's plot the ISI histogram for our own example neuron:
+This result is the following:
 <code matlab>
-iC = 47;
+>> evt
-spk_t = S.t{iC}; % spike times
-isi = diff(spk_t); % interspike intervals
-dt = 0.001; % in s, because spike times are in s
+evt =
-isi_edges = 0:dt:0.25; % bin edges for ISI histogram
-isi_centers = isi_edges(1:end-1)+dt/2; % for plotting
-isih = histc(isi,isi_edges);
+        t: {1x8 cell}
+    label: {'n0'  'n1'  'c1'  'c3'  'c5'  'd1'  'd3'  'd5'}
+      cfg: [1x1 struct]
+>> evt.t
+ans =
+  Columns 1 through 4
+    [1x58 double]    [1x57 double]    [1x35 double]    [1x19 double]
+  Columns 5 through 8
+    [1x34 double]    [1x40 double]    [1x25 double]    [1x42 double]
-bar(isi_centers,isih(1:end-1)); % remember to ignore the last bin of histc() output
-set(gca,'FontSize',20,'XLim',[0 0.25]); xlabel('ISI (s)'); ylabel('count'); grid on;
 </code>
-You should get
+Thus, for each item in ''cfg.eventList'' we extract the corresponding event timestamps, placing them in an output structure field with a name specified in ''cfg.eventLabel''.
+The events returned by ''LoadEvents()'' above are already present in the Events file; this function merely presents them in a more usable format. However, for many tasks, some more processing is needed to obtain event times of interest, for instance "only those nosepokes that follow presentation of a cue predicting 5 pellets, and for which the rat does not un-poke for at least one second". Creating appropriate "trial functions" that return usable event times is unglamorous, but a key analysis step for any data set. (Remember: "Garbage in, garbage out"!)
-{{ :analysis:course:week9_fig6.png?600 |}}
+In any case, armed with some event times, we are now in a position to ask questions of the sort, "does reward delivery elicit changes in the ventral striatal LFP?" and "is the ventral striatal LFP sensitive to reward-predictive cues?" To do so, we will use the %%FieldTrip%% toolbox, which massively expands MATLAB's suite of built-in functions for spectral analysis.
-This histogram, which simply counts how many interspike intervals there are in this spike train for every interval from 0-250ms, shows that the most common interspike intervals are around 5ms. There is a hint of an increase in interspike intervals between 100-150ms, in the theta range typical for a hippocampal "place cell". This ISI distribution can also be seen in the spike train snippet we used for computing firing rates, above: the short ISIs occur in bursts, and then there is a longer interval (in the theta range) until the next burst.
+==== Event-triggered spectrograms using FieldTrip ====
-Clearly, this cell does not emit spikes at random intervals; instead, there is structure in the spike train. A further indication of this is the //ISI return plot//, a scatterplot of each interspike interval against the next.
+[[http://fieldtrip.fcdonders.nl/ | FieldTrip]] is a toolbox originally developed for the analysis of human EEG and MEG data. Because these settings have historically faced complex signal processing challenges, such as eyeblink artifact removal, source reconstruction, and mapping across many electrodes, the methods developed to address these challenges are mature and have much to offer to rodent electrophysiologists.
-☛ Using the ''isi'' variable created above, plot each ISI against its successor. Set the axes to run from 0 to 0.25 seconds.
+The FieldTrip toolbox continues to be actively developed and maintained, primarily at the [[http://www.ru.nl/donders/ | Donders Institute for Brain, Cognition and Behaviour]] in Nijmegen. It is supported by a popular [[http://fieldtrip.fcdonders.nl/discussion_list | mailing list]], extensive documentation with [[http://fieldtrip.fcdonders.nl/tutorial | tutorials]], and various workshops organized across the globe (such as [[http://fieldtrip.fcdonders.nl/workshop/toronto | this one]] in Toronto).
-You should get something like this:
+=== FieldTrip setup ===
-{{ :analysis:course:week9_fig7.png?600 |}}
+☛ If you haven't already done so, clone the lab's [[https://github.com/mvdm/fieldtrip|FieldTrip repository on GitHub]]. We currently run a custom version modified to allow the processing of Neuralynx data, so the official %%FieldTrip%% release will not work! See [[analysis:nsb2014:week1|Module 1]] for more details on setting this up.
-This plot shows that a short ISI (in the 5ms range) tends to be followed either by another short ISI, or by a long one in the theta range. It however seems rare for a theta-range ISI to be followed by another theta-range ISI, again in accordance with our visual impression of the spike train.
+☛ Also make sure you have a MATLAB shortcut to create a path that includes the nsb-2014 and fieldtrip code (but not %%MClust%%!). Run this shortcut before proceeding.
+=== Loading Neuralynx data into FieldTrip ===
-=== Interlude: Poisson point process ===
+Let's load our familiar ventral striatal LFP into FieldTrip:
-A useful tool for many spike train analysis questions is to make a comparison with a random, synthetic spike train. Such comparisons can bring into focus salient aspects of real spike trains, and can be used for tests of significance.
+<code matlab>
+%% remember to cd to your data folder
+fc = {'R016-2012-10-03-CSC04a.ncs'};
+data = ft_read_neuralynx_interp(fc);
+</code>
-A simple way of generating a random spike train is to essentially flip a coin at each time step, scoring a spike for heads, and and no spike for tails:
+You will get a number of warnings about discontinuous recording and expected sample mismatches. This is because FieldTrip expects an exactly constant ''Fs'' in the data, which we know Neuralynx typically does not provide.
+Next, let's look at the contents of the ''data'' structure:
 <code matlab>
-dt = 0.001;
+>> data
-t = [0 10]; % time interval (length) of spike train to generate
-tvec = t(1):dt:t(2);
+data =
-pspike = 0.5; % probability of generating a spike in bin
+           hdr: [1x1 struct]
-rng default; % reset random number generator to reproducible state
+         label: {'R016-2012-10-03-CSC04a'}
-spk_poiss = rand(size(tvec)); % random numbers between 0 and 1
+          time: {[1x5256982 double]}
-spk_poiss_idx = find(spk_poiss < pspike); % index of bins with spike
+         trial: {[1x5256982 double]}
-spk_poiss_t = tvec(spk_poiss_idx)'; % use idxs to get corresponding spike time
+       fsample: 2000
+    sampleinfo: [1 5256982]
+           cfg: [1x1 struct]
-line([spk_poiss_t spk_poiss_t],[-1 -0.5],'Color',[0 0 0]); % note, plots all spikes in one command
-axis([0 0.1 -1.5 5]); set(gca,'YTick','');
 </code>
-This gives:
+A quirk of FieldTrip is that it treats all data as organized into //trials//, a legacy of its event-related potential roots in human work. However, a quick ''plot(data.time{1},data.trial{1})'' should assure you that all the data is there. At this point, FieldTrip thinks of this full recording session as one single "trial", the first cell of ''data.label'', ''data.time'', and ''data.trial''.
-{{ :analysis:course:week9_fig8.png?600 |}}
+=== Creating trial data in FieldTrip ===
-This spike train looks different from our real one.
+To convert this single "trial" containing the entire recording into proper trials surrounding events of interest, %%FieldTrip%% requires a ''trl'' variable that specifies the start and end indices (into the data) which are then passed to the ''ft_redefinetrial()'' function, which cuts up the data into trials:
-☛ Plot its ISI histogram. Some differences are apparent; for instance, this synthetic spike train clearly has much shorter ISIs overall.
+<code matlab>
+cfg = [];
+cfg.t = cat(2,getd(evt,'n0'),getd(evt,'n1'));
+cfg.mode = 'nlx';
+cfg.hdr = data.hdr;
+cfg.twin = [-1 4];
-This difference occurs because the mean firing rates between the two spike trains are not equal. Our real spike train has a mean firing rate of about 0.47 spikes/s (as ''1./mean(diff(spk_t))'' will verify). In contrast, our synthetic spike train has a theoretical mean firing rate of 500 spikes/s (50% probability of a spike in each 1ms bin).
+trl = ft_maketrl(cfg);
-☛ What should the probability per 1ms bin of a spike be, in order for our synthetic spike train to have a mean firing rate of 0.47 spikes/s? Recompute the ISI histogram for a spike train generated with this probability instead of 0.5. You will need to increase the length of the spike train to 4500s in order to get a number of spikes which is similar to that in the real spike train.
+cfg = [];
+cfg.trl = trl;
+data_trl = ft_redefinetrial(cfg,data);
+</code>
-You should get something like:
+The above reorganizes the data into 4-second windows (from -1 to 4 seconds) around nosepoke times (''evt.n0'' and ''evt.n1''). ''ft_maketrl()'' converts the event times from ''LoadEvents()'' into the data indices that ''ft_redefinetrial()'' needs.
-{{ :analysis:course:week9_fig9.png?600 |}}
+☛ Inspect the ''data_trl'' structure. Notice that there are now 115 trials instead of 1, and that each trial has its own time and data vectors. Plot the 47th trial as an example to verify this.
-The number of "short" ISI counts (i.e. the ones shown in this plot) is way smaller in our synthetic spike train than in the real one. So, even though the mean firing rates are the same, our real spike train has relatively many short ISIs, which must be balanced with longer ones to end up with the same mean.
+=== Constructing the event-triggered spectrogram ===
-Generating artificial spike trains in this way, by flipping a coin at each time step, is known as a //Poisson point process// (Poisson process for short): it generates interspike intervals drawn from the [[http://en.wikipedia.org/wiki/Poisson_distribution | Poisson distribution]]. An important property of this distribution is that its values are never negative. This is useful when dealing with counts and interspike intervals, which cannot be negative. For large means, the Poisson distribution will approach a Gaussian distribution, but a Gaussian with a mean close to zero is likely to generate negative values, which are nonsensical when dealing with ISIs or spike counts.
+Now we are ready to construct the event-triggered spectrogram:
-A further property of the Poisson process that each ISI is independent of the last.
+<code matlab>
+cfg              = []; % start with empty cfg
+cfg.output       = 'pow';
+cfg.channel      = 'R016-2012-10-03-CSC04a';
+cfg.method       = 'mtmconvol';
+cfg.taper        = 'hanning';
+cfg.foi          = 10:2:100; % frequencies of interest
+cfg.t_ftimwin    = ones(size(cfg.foi)).*0.5;  % window size: fixed at 0.5s
+cfg.toi          = -0.5:0.05:3.5; % times of interest
-☛ Generate the ISI return plot for your Poisson spike train. If you only see a few spikes, increase the axes limits.
+TFR = ft_freqanalysis(cfg, data_trl);
-You should see a total absence of structure (correlations) in the plot, unlike the equivalent for the real spike train, indicating that a given ISI does not tell you anything about what the next one might be.
+figure
+cfg = []; cfg.channel = 'R016-2012-10-03-CSC04a';
+ft_singleplotTFR(cfg, TFR);
+</code>
-A final difference with our actual spike train is that pure Poisson spike trains have no refractory period. Verify that the ISIh of the real spike train does!
+That is quite a mouthful. Almost any FieldTrip function requires a ''cfg'' struct as input that specifies the parameters of the task to be completed. To avoid unexpected results, it is good practice to start with an empty cfg (''cfg = [];'') when moving on to a new type of operation, such as from trial segmentation to spectrogram construction in the above example. The function ''ft_freqanalysis()'' computes the spectrogram based on the data and the parameters included in the ''cfg'' file. ''ft_singleplotTFR()'' plots the result, again based on the information in the ''cfg'' file and the output of ''ft_freqanalysis()''.
-Our synthetic Poisson spike train had a fixed probability of emitting a spike in each time bin. This results in a //homogenous// Poisson process, meaning that theoretically, the firing rate does not vary over time. Of course, there will be random fluctuations, but the underlying probability of a spike remains constant. In contrast, an //inhomogenous// Poisson process involves changes over time of the underlying spike probability.
+You should get:
-☛ In the above code for generating a spike train, ''pspike'' was set to a fixed probability of 0.5. Change ''pspike'' to be a vector of the same size as ''tvec'' describing a 1 Hz sinusoid which oscillates between 0 and 0.5. Use this changing spike probability to generate your artificial (inhomogenous Poisson) spike train. Verify by visual inspection that indeed the firing rate of this spike train is varying over time.
+{{ :analysis:nsb2014:week6_fig4.png?600 |}}
-==== Spike autocorrelation function (acf) ====
+This plot is an //average// spectrogram over 115 trials, with time 0 corresponding to the time of nosepoking into the reward receptacle. As with the examples using ''spectrogram()'', we can control the smoothness of the spectrogram by specifying the time and frequency steps. For the frequency steps this is accomplished the same way in FieldTrip, but for the time axis we can now simply specify a vector with the window centers instead of a window size and overlap.
-A different way of characterizing the properties of a spike train is to compute its autocorrelation function, which describes the probability of finding a spike at time $t+t'$ given a spike at time $t$, for some range of lags $t'$. This is analogous to the autocorrelation of gamma-band power in a ventral striatal LFP that we computed in the previous module.
+Looking at the results, it appears there are some beta (~20Hz) and gamma events visible in the average spectrogram.
-MATLAB doesn't provide a built-in function that can take in a spike train and compute its autocorrelation, so we have to make our own. Here is a simple implementation:
+☛ Change the ''cfg.toi'' and ''cfg.foi'' parameters to verify you can change the resolution of the time and frequency axes. What happens if you try to extend the time axis from -2 to 5? What would you do if you wanted data for that extended range?
-<code matlab>
+=== Options and parameters for FieldTrip spectrograms ===
-function [ac,xbin] = acf(spk_t,binsize,max_t)
-% function [ac,xbin] = acf(spike_times,binsize,max_t)
-%
-% estimate autocorrelation of input spike train
-%
-% INPUTS:
-% spike_times: [1 x nSpikes] double of spike times (in s)
-% binsize: acf bin size in s
-% max_t: length of acf in s
-%
-% OUTPUTS:
-% ac: autocorrelation estimate (normalized so that acf for the zero bin is
-% 0)
-% xbin: bin centers (in s)
-%
-% MvdM 2013
-xbin_centers = -max_t-binsize:binsize:max_t+binsize; % first and last bins are to be deleted later
+FieldTrip provides many powerful tools for spectral analysis. This section highlights a few useful ones, but the [[http://fieldtrip.fcdonders.nl/tutorial/timefrequencyanalysis | tutorials]] and documentation of the functions involved (e.g. [[http://fieldtrip.fcdonders.nl/reference/ft_freqanalysis | ft_freqanalysis]]) give a more complete overview.
-ac = zeros(size(xbin_centers));
-for iSpk = 1:length(spk_t)
+== Baseline correction ==
-   relative_spk_t = spk_t - spk_t(iSpk);
-   ac = ac + hist(relative_spk_t,xbin_centers); % note that hist() puts all spikes outside the bin centers in the first and last bins! delete later.
-end
-xbin = xbin_centers(2:end-1); % remove unwanted bins
+Often we are not interested in absolute power levels, as provided per default by a spectrogram. Rather, we wish to compare spectral content following some event of interest to a baseline. This is easily accomplished:
-ac = ac(2:end-1);
-ac = ac./max(ac); % normalize
+<code matlab>
+figure
+cfg = [];
+cfg.baseline     = [-2 0];
+cfg.baselinetype = 'relative';
+cfg.channel      = 'R016-2012-10-03-CSC04a';
+ft_singleplotTFR(cfg, TFR);
 </code>
-☛ Using this autocorrelation function and the code for generating Poisson spike trains above, compute and plot the autocorrelation of a Poisson spike train with a mean firing rate of 0.47 Hz (exactly as done above), a 10ms time bin, from -1 to 1 seconds. Note that the ''acf()'' function requires a ''ts'' object as input, so you will need to convert the spike times to a ts object using ''ts(spk_poiss_t)''.
+Note that the beta and gamma events now stand out clearly.
-You should get something like:
+== Frequency-dependent windowing ==
-{{ :analysis:course:week9_fig10.png?600 |}}
+Up until now, we have used a //fixed// window size (number of samples, or equivalently, time) in computing spectrograms. The above example specified a 0.5s window.
-The autocorrelation of a Poisson spike train appears to be essentially flat over different time lags.
+☛ Recompute and plot the event-triggered spectrogram using a //frequency-dependent// window of 20 cycles by setting ''cfg.t_ftimwin = 20./cfg.foi;'' and a frequency range of 1:100.
-☛ Why? What differences to you notice with the acorr of the actual spike train (''[acorr,xbin] = acf(S{47},0.01,1);'')?
+You should now have a beautiful image that balances time and frequency resolution:
-==== Spike cross-correlation function (ccf) ====
+{{ :analysis:nsb2014:ft_specgram2.png?600 |}}
-The intuition behind the cross-correlation between two spike trains is the following: given that neuron 1 emits a spike at time $t$, what is the probability of seeing a spike from neuron 2 at different lags $t+t'$? For instance, if neuron 2 reliably fires 10ms after neuron 1, the cross-correlation between 1 and 2 would have a peak at +10ms. The autocorrelation is a special case of this with the two spike trains being the same.
+☛ What causes the whitespace at the bottom of the image?
-Computing the cross-correlation is a simple modification of our autocorrrelation function:
+This frequency-dependent windowing can be thought of as computing many spectrograms with different windows in parallel, and then keeping different pieces of each so that low frequencies are estimated with good precision in the frequency domain (but poor temporal precision) and high frequencies are resolved from good temporal precision (but relatively poor frequency precision). This tradeoff is illustrated by the time-frequency tiles in the "B" panel of the figure earlier in this module.
+(As a technical note, this frequency-dependent windowing is similar to what can be obtained with wavelet transforms.)
+== Statistical tests ==
+To set up a statistical comparison of the pre-nosepoke baseline and the oscillation patterns following the nosepoke, we perform two different frequency analyses, with different times of interest:
 <code matlab>
-function [cc,xbin] = ccf(spike_times1,spike_times2,binsize,max_t)
+cfg              = [];
-% function [cc,xbin] = ccf(spike_times1,spike_times2,binsize,max_t)
+cfg.output       = 'pow';
-%
+cfg.channel      = 'R016-2012-10-03-CSC04a';
-% estimate cross-correlation function of input spike train
+cfg.method       = 'mtmconvol';
-%
+cfg.taper        = 'hanning';
-% INPUTS:
+cfg.foi          = 1:1:100;
-% spike_times1: [1 x 1 ts]
+cfg.keeptrials   = 'yes'; % need this for stats later
-% spike_times2: [1 x 1 ts]
+cfg.t_ftimwin    = 20./cfg.foi;  % 20 cycles per time window
-% binsize: ccf bin size in s
-% max_t: length of acf in s
-%
-% OUTPUTS:
-% cc: cross-correlation estimate (relative to spike_times1)
-% xbin: bin centers (in s)
-%
-% MvdM 2013
-spk_t1 = Data(spike_times1);
+cfg.toi          = -1:0.05:0; % pre-nosepoke baseline
-spk_t2 = Data(spike_times2);
+TFR_pre = ft_freqanalysis(cfg, data_trl);
+cfg.toi          = 0:0.05:1; % post-nosepoke
+TFR_post = ft_freqanalysis(cfg, data_trl);
+TFR_pre.time = TFR_post.time; % time axes should be identical for comparison
+</code>
+Now we run ''ft_freqstatistics()'' to do the comparison between the two frequency analyses:
+<code matlab>
+%% t-test
+cfg = [];
+cfg.channel     = 'R016-2012-10-03-CSC04a';
+cfg.latency     = 'all';
+cfg.trials      = 'all';
+cfg.frequency   = 'all';
+cfg.avgoverchan = 'no';
+cfg.avgovertime = 'no';
+cfg.avgoverfreq = 'no';
+cfg.parameter   = 'powspctrm';
+cfg.method      = 'stats';
+cfg.statistic   = 'ttest2';
+cfg.alpha       = 0.05;
-xbin_centers = -max_t-binsize:binsize:max_t+binsize; % first and last bins are to be deleted later
+nTrials1 = size(TFR_pre.powspctrm,1); nTrials2 = size(TFR_post.powspctrm,1);
-cc = zeros(size(xbin_centers));
-for iSpk = 1:length(spk_t1)
+cfg.design = cat(2,ones(1,nTrials1),2*ones(1,nTrials2)); % two conditions
+cfg.ivar = 1; % dimension of design var which contains the independent variable (group)
-   relative_spk_t = spk_t2 - spk_t1(iSpk);
-   cc = cc + hist(relative_spk_t,xbin_centers); % note that histc puts all spikes outside the bin centers in the first and last bins! delete later.
-end
-xbin = xbin_centers(2:end-1); % remove unwanted bins
+stat = ft_freqstatistics(cfg,TFR_post,TFR_pre);
-cc = cc(2:end-1);
-cc = cc./(length(spk_t1)); % normalize by number of spikes of first input
+cfg.parameter = 'stat';
+ft_singleplotTFR(cfg,stat); % plot the t-statistic
 </code>
-The only potentially tricky part of cross-correlations is how to do the normalization (there are some different choices that determine if the output is a conditional probability, a correlation, or something else), but we will not worry about this for now. All these normalizations produce outputs that are qualitatively similar.
+The result shows the t-statistic of the bin-by-bin comparison between pre- and post-nosepoke power. Depending on what correction for multiple comparisons is to be applied, some time-frequency bins may approach significance. This is however only a single recording session.
-☛ What do you expect the cross-correlation function of two Poisson spike trains to look like? Why? Verify your intuition.
+The Statistics section of [[http://fieldtrip.fcdonders.nl/walkthrough | this walkthrough]] has some background on the different steps performed here; the [[http://fieldtrip.fcdonders.nl/tutorial/eventrelatedstatistics || event-related statistics]] tutorial has further examples of different statistics.
-Let's compute the cross-correlation between two hippocampal place cells, simultaneously recorded from a rat running a T-maze:
+== Multichannel data ==
+FieldTrip excels at processing many data channels in parallel. Processing all 16 channels from this session is as simple for the user as processing 1 channel. (Be aware however, that this will be mild torture for older computers; the following took about 5 minutes to run on my 3-year old laptop.)
 <code matlab>
-cell1_id = 5; cell2_id = 42;
+%% load the data
+fc = FindFiles('*.ncs'); % get filenames of all LFPs recorded
+data = ft_read_neuralynx_interp(fc); % load them all -- this will take a while
+data_all.hdr.Fs = data_all.fsample; % for some reason this is missing from the header
-s1 = Restrict(S{cell1_id},3200,5650); % restrict to on-track times only
+%% define layout for later plotting
-s2 = Restrict(S{cell2_id},3200,5650);
+cfg = [];
+cfg.layout = 'ordered'; cfg.channel = data.label;
+layout = ft_prepare_layout(cfg, data);
+</code>
-[xcorr,xbin] = ccf(s1,s2,0.01,1);
+Defining a channel layout is the only addition when dealing with multichannel data. Here, we are using one of the default layouts ('ordered'). However, layouts specific to particular EEG/MEG and intracranial probes are also available, so that results can be plotted topographically (i.e. matching the physical sensor layout).
-plot(xbin,xcorr);
+Remarkably, the trial definition and spectrogram computations are identical to those for a single channel:
-set(gca,'FontSize',20); xlabel('lag (s)'); ylabel('xcorr');
-title(sprintf('%d-%d',cell1_id,cell2_id));
-</code>
-You should get:
+<code matlab>
+%%
+cfg = [];
+cfg.t = cat(2,getd(evt,'n0'),getd(evt,'n1'));
+cfg.mode = 'nlx';
+cfg.hdr = data_all.hdr;
+cfg.twin = [-1 4];
-{{ :analysis:course:week9_fig11.png?600 |}}
+trl = ft_maketrl(cfg);
-This is an example of one of my favorite plots in neuroscience. It holds the key to revealing a fundamental and unique property of the rodent hippocampus, thought to reflect a specialization for the rapid encoding of episodic memories. The plot shows (something akin to) the probability of cell 42 spiking at various time lags (between -1 and 1 second) given that cell 5 spikes at time 0. This ccf is clearly asymmetric, in that cell 42 is more likely to spike after cell 5 than before, as evident from the crosscorrelation values on the right side of the plot (5-->42) being larger overall than those on the left (42-->5).
+cfg = [];
+cfg.trl = trl;
+data_trl = ft_redefinetrial(cfg,data_all);
-There is something else going on, too: on top of this slow component of the ccf, there is a faster component with sharp peaks repeating at theta frequency (about 8-9 per second). The peak closest to zero is in fact slightly offset, occurring at approximately 30-40ms. Thus, cell 42 is likely to spike just a few tens of ms after cell 5; in contrast, the reverse, where cell 5 spikes a few tens of ms after cell 42, almost never happens.
+%%
+cfg              = []; % start with empty cfg
+cfg.output       = 'pow';
+cfg.method       = 'mtmconvol';
+cfg.taper        = 'hanning';
+cfg.foi          = 1:100; % frequencies of interest
+cfg.toi          = -0.5:0.05:3.5; % times of interest
+cfg.t_ftimwin = 20./cfg.foi;
+TFR = ft_freqanalysis(cfg, data_trl);
+</code>
-As it turns out, this highly precise ccf shape results from the fact that each theta cycle in the hippocampal LFP is in fact associated with the sequential activation of a number of place cells, tracing out a trajectory on the maze, as illustrated in this figure from [[http://itb.biologie.hu-berlin.de/~kempter/HippoJC/Articles/Phase_Precession/skaggs96.pdf | Skaggs et al. 1996]]:
+Only the final plot command is slightly different:
-{{ :guides:skaggs_mcnaughton_pp.png?600 |}}
+<code matlab>
+figure
+cfg = [];
+cfg.baseline     = [-2 0];
+cfg.baselinetype = 'relative';
+cfg.layout = layout;
-A different way of showing this idea is as follows (from [[http://www.vandermeerlab.org/SM_RWC_MvdM_pp.pdf | Malhotra et al. (2012)]]:
+ft_multiplotTFR(cfg, TFR); % note this is now multiplot rather than singleplot
+</code>
-{{ :analysis:course:week9_fig12.png?800 |}}
+You should get:
-This figure illustrates the ccf (labeled 'CCG' here, for cross-correlogram, which is the same thing) that would be expected from generating inhomogenous Poisson spike trains based on a firing rate profile given by two place fields (the blue and red lines in the top panel). In this case, with the animal running from left to right, the blue spikes would tend to come after the red, because the red field is entered before the blue field. This is reflected in the asymmetric shape of the ccf. However, this lacks the sharp theta peaks apparent in the ccf we just computed from our two real place cells. The figure here shows that if we enable an effect called "phase precession" (see [[guides:ppintro|here]] for more optional info on this) then the shape of the ccf does reproduce what is found experimentally. You can also see that accordingly, the blue spikes tend to come just after the red spikes within each theta cycle, forming sequences of place cells.
+{{ :analysis:nsb2014:ft_specgram3.png?600 |}}
-An important question is to establish if such sequences could come about by just having each place cell //independently// emit spikes in a manner described by some average spike density function, or if the sequences that are actually seen require some sort of //coordination// between place cells.
+This shows the baseline-corrected, event-aligned spectrograms for all 16 channels in this session. There is also an average shown (the subplot on the lower right). Note that a subset of the channels, presumably located in the ventral striatum, show a very similar time-frequency pattern, whereas another subset do not show this at all. The channels that do not show this same pattern are likely located in the hippocampus.
+==== Exercise ====
-To test this, we can first convert our spike trains to a SDF over time, and then use that SDF to generate inhomogenous Poisson spike trains. Since these two spike trains are generated by independent coin flips, i.e. whether or not there was a spike for neuron 1 cannot directly affect the probability of a spike in neuron 2, we can determine if this is sufficient to reproduce the ccf seen in the real data (in which there may be some dependency between neuron 1 and neuron 2).
+Because the event-triggered spectrograms above show an average over trials, it is possible to lose touch with the properties of the raw data that generated it. For instance, the average may arise from a distribution of very similar looking individual trials, or it may be dominated by one atypical but extreme trial. For this reason it is important to compare the spectrogram to the raw data.
-=== Exercise ===
+An example of a study where this was done well is [[http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3463873/ | Leventhal et al]]. Neuron 2012. [[http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3463873/figure/F3/ | Figure 3]] in this paper shows a number of event-triggered LFP traces, one for each trial.
-☛ Compute the ccf of two spike trains generated from the SDF of our two hippocampal place cells (5 and 42) and compare it to their actual ccf. Use a 50 ms standard deviation for the Gaussian, and 1ms bin size. The sequence of steps for this would look something like this:
+☛ Implement a function eventLFPplot(cfg,csc) to make such a plot. The ''cfg'' variable should contain an ''eventTimes'' field which is required, and optionally accept a ''window'' field to override the default window of [-1 3] seconds relative to the event times.
+You can use a skeleton like the following:
 <code matlab>
-% for each of the two neurons, restrict the data to [3200 5650] (the time interval when the rat was running on the track)
+% for each event:
+% *extract the corresponding piece of LFP
+% *replace the original time axis with a new one based on the time window asked for
+% *optionally, rescale the LFP
-% compute the spike density function for each, making sure that your tvec runs from 3200 to 5650 also, and that you have a 50ms SD for the Gaussian convolution kernel
+% *add a y-offset to the LFP to plot one above the other
-% to use these SDFs to generate Poisson spike trains, convert the firing rates given by the SDF to a probability of emitting a spike in a given bin. (As you did above for a 0.47 Hz constant firing rate.)
+% *plot the current piece of LFP
-% generate Poisson spike trains, making sure to use the same tvec
+% add a vertical line to the plot indicating time zero
-% convert Poisson spike trains to ts objects and compute the ccf
+% further niceties: add an option to decimate, to color, to filter before plotting
 </code>
-(This comparison of independently generated spikes with the actual data is the essence of [[http://www.ncbi.nlm.nih.gov/pubmed/17663452 | this]] excellent paper.)
+Test your code with the same events used to generate the event-triggered spectrogram above. Do you think the average is a reasonable charaterization of what is happening in the LFP for these events?

User Tools

Site Tools

Differences

Page Tools