wiki

~~DISCUSSION~~ ===== Spike train analysis II: tuning curves, encoding, decoding ===== Goals: * Understand the format of Neuralynx video tracking data * Learn to plot place cell tuning curves (rate maps), raw and smoothed * Implement a basic Bayesian decoding algorithm * Compare decoded and actual position by exporting to a movie file Resources: * (if you have not encountered Bayes' rule and conditional probability before) [[http://www.cs.ubc.ca/~murphyk/Bayes/bayesrule.html | A brief introduction ]] by Ken Murphy at UBC * (for reference) [[http://jn.physiology.org/content/79/2/1017.long | Zhang et al. 1998]], first application of decoding to place cell data, with nice explanations and derivations * (for reference) [[http://www.jneurosci.org/content/18/18/7411.long | Brown et al. 1998]], an example of a more sophisticated decoding method Assignment: * Design and implement a measure of decoding accuracy **Note**: this module uses an externally compiled function, ''ndhist''. By default, this will not work on non-Windows machines, but the source code is provided. If you are a non-Windows user and able to compile this, I would be grateful if you could share the binary here. ==== Introduction ==== To support adaptive behavior, activity in the brain must correspond in some way to relevant sensory events and planned movements, combine many sources of information into multimodal percepts, and recall traces of past events to inform predictions about the future. In other words, neural activity must somehow //encode// relevant quantities. For instance, it can be [[http://www.ncbi.nlm.nih.gov/pubmed/9817202 | demonstrated]] [[http://www.ncbi.nlm.nih.gov/pubmed/20141292 | behaviorally]] that many animals use estimates of their location and head direction to navigate towards a goal. Where, and how, are these quantities represented in the brain? What are the neural circuits that can compute and update these signals? How do place and direction estimates contribute to which way to go? This //information processing// view of the brain has been extremely influential, as highlighted by the enduring appeal of [[http://www.ncbi.nlm.nih.gov/pubmed/4966457 | Hubel and Wiesel's demonstrations]] that single cells in macaque V1 respond to bars of light not only within a particular region of visual space, but also with a specific orientation. Such cells are said to be //tuned// for orientation [of the bar] and a typical //tuning curve// would therefore look like this: {{ :analysis:course:biol377_2-3.png?300 |}} This tuning curve describes how the cell responds, on average, to different orientations of the stimulus. If the cell were to respond with the same firing rate across the range of stimulus orientations, then the cell is indifferent to this particular stimulus dimension: it does not encode it. However, because this cell clearly modulates its firing rate with stimulus orientation, it encodes, or represents (I use these terms interchangeably, but some disagree) this quantity in its activity. We can turn this idea around and note that if orientation is encoded, this implies we can also //decode// the original stimulus from the cell's activity. For instance, if we noted that this cell was firing at a high rate, we would infer that the stimulus orientation is likely close to the cell's preferred direction. Note that this requires knowledge of the cell's tuning curve, and that based on one cell only, we are unlikely to be able to decode (or reconstruct, which means the same thing) the stimulus perfectly. The more general view is to say that the cell's activity provides a certain amount of information about the stimulus, or equivalently, that our (decoded) estimate of the stimulus is improved by taking the activity of this cell into account. This module first explores some practical issues in estimating tuning curves of "place cells" recorded from the rat hippocampus. An introduction to a particular decoding method (Bayesian decoding) is followed by application to many simultaneously recorded place cells as a rat performs a T-maze task. ==== Estimating place cell tuning curves (place fields) ==== First, load the "place cell" data set also used in the previous module, which contains a number of spike trains recorded simultaneously from the dorsal CA1 area of the hippocampus: <code matlab> cd 'D:\My_Documents\My Dropbox\teaching\BIOL6xx\R042-2013-08-18'; sd.fc = FindFiles('*.t'); sd.fc = cat(1,sd.fc,FindFiles('*._t')); % also load poorly isolated cells (if you don't have these, get from database) sd.S = LoadSpikes(sd.fc); </code> Note that we are loading ''*._t'' files in addition to the familiar ''*.t'' spike time files: the underscore indicates a cell with questionable isolation quality, likely contaminated with noise, spikes from other neurons, and/or missing spikes. In general, you do not want to use such neurons for analysis, but in this case we are not concerned with properties of individual neurons. We are instead interested in the information present in a population of neurons, and for this we will take everything we can get. The spike data is loaded into a ''sd'' ("SessionData") structure to make it easier to add video tracking data to it while keeping everything organized. === Interlude: Neuralynx video data === Because we will be examining representations of place (location) in the hippocampus, we need to know the location of the rat. This is given by the ''.nvt'' video tracker (VT) file: <code> % remember to unzip VT1.zip if not yet done [Timestamps, X, Y, Angles, Targets, Points, Header] = Nlx2MatVT('VT1.nvt', [1 1 1 1 1 1], 1, 1, []); Timestamps = Timestamps*10^-6; % remove zero samples, which happen when no LED is detected toRemove = (X == 0 & Y == 0); X = X(~toRemove); Y = Y(~toRemove); Timestamps = Timestamps(~toRemove); sd.x = tsd(Timestamps,X'); sd.y = tsd(Timestamps,Y'); </code> Video tracker data relies on the detection of LEDs, mounted on the headstage which is plugged in to the rat's electrode array ("drive") during behavior. The Neuralynx software can be configured to detect multiple targets of different colors, but by default a threshold based on luminance is applied, and the center of mass of all suprathreshold pixels found. This is returned in the ''X'' and ''Y'' variables above. As with raw CSC data, the timestamps are in microseconds, so must be converted to get seconds. ☛ Using the ''Timestamps'' variable, determine the sampling rate of the video data. Beware of gaps in the data! If no pixels pass the detection threshold, Neuralynx returns a 0 position. These are removed from the data before further processing to prevent confusion. === Visual inspection === Before looking at the data we will first exclude the pre- and post-run segments: <code matlab> t = [3250 5650]; for iC = 1:length(sd.S) sd.S{iC} = Restrict(sd.S{iC},t(1),t(2)); end sd.x = Restrict(sd.x,t(1),t(2)); sd.y = Restrict(sd.y,t(1),t(2)); </code> Now we can plot the position data: <code matlab> plot(Data(sd.x),Data(sd.y),'.','Color',[0.5 0.5 0.5],'MarkerSize',1); axis off; hold on; </code> Next, we plot the spikes of a single cell (cell 7, which we met in the previous module) where the rat was when each spike was emitted: <code matlab> % get x and y coordinate for times of spike iC = 7; spk_x = interp1(Range(sd.x),Data(sd.x),Data(sd.S{iC}),'linear'); spk_y = interp1(Range(sd.y),Data(sd.y),Data(sd.S{iC}),'linear'); h = plot(spk_x,spk_y,'.r'); </code> Note the use of ''interp1()'' here: it finds the corresponding ''x'' and ''y'' values for each spike time using linear interpolation. You should see: {{ :analysis:course:week10_scatterfield.png?600 |}} This cell seems to have a place field just to the left of the choice point on the track. There are also a few spikes on the pedestals, where the rat rests in between runs on the track. This figure is a useful visualization of the raw data, but it is not a tuning curve. === Estimating tuning curves === A simple way to obtain a tuning curve is to bin the spikes (not in time, as in the previous module, but now in space): <code matlab> SET_xmin = 10; SET_ymin = 10; SET_xmax = 640; SET_ymax = 480; SET_nxBins = 63; SET_nyBins = 47; spk_binned = ndhist(cat(1,spk_x',spk_y'),[SET_nxBins; SET_nyBins],[SET_xmin; SET_ymin],[SET_xmax; SET_ymax]); imagesc(spk_binned'); axis xy; colorbar; </code> ''ndhist'', for "N-dimensional histogram" can bin data in more than one dimension, in this case a 2-D grid with ''x'' and ''y'' position as the dimensions. It needs to know the minimum and maximum values for each dimension, as well as the number of bins in each; the data should be formatted with one dimension per row. This gives: {{ :analysis:course:week10_tc0.png?600 |}} This spike count can then be divided by how much time the rat spent in each of the bins (the "occupancy", a value in seconds) to get a firing rate: <code matlab> occ_binned = ndhist(cat(1,Data(sd.x)',Data(sd.y)'),[SET_nxBins; SET_nyBins],[SET_xmin; SET_ymin],[SET_xmax; SET_ymax]); % this is a sample count, so need to convert to seconds (1/30s per sample) to get time VT_Fs = 30; tc = spk_binned./(occ_binned .* (1./VT_Fs)); % firing rate is spike count divided by time pcolor(tc'); shading flat; axis xy; colorbar; axis off; </code> Note that the occupancy is obtained by multiplying the number of samples in each bin with the sampling interval. You should get: {{ :analysis:course:week10_tc1.png?600 |}} ''pcolor()'' ensures that bins which are NaN (such as would occur from division by zero occupancy) show up as white; ''imagesc()'' shows NaNs the same way as zeros. We now have a tuning curve for our cell, sometimes also called a "rate map" because we are dealing with (x,y) position. However, it suffers from the same issues arising from binning identified in the previous module. An improvement would be to bin at a much finer scale, and then use smoothing: <code matlab> SET_nxBins = 630; SET_nyBins = 470; % more bins kernel = gausskernel([30 30],8); % 2-D gaussian for smoothing: 30 points in each direction, SD of 8 bins % spikes spk_binned = ndhist(cat(1,spk_x',spk_y'),[SET_nxBins; SET_nyBins],[SET_xmin; SET_ymin],[SET_xmax; SET_ymax]); spk_binned = conv2(spk_binned,kernel,'same'); % occupancy occ_binned = ndhist(cat(1,Data(sd.x)',Data(sd.y)'),[SET_nxBins; SET_nyBins],[SET_xmin; SET_ymin],[SET_xmax; SET_ymax]); occ_binned = conv2(occ_binned,kernel,'same'); occ_binned(occ_binned < 0.01) = 0; tc = spk_binned./(occ_binned .* (1 / VT_Fs)); tc(isinf(tc)) = NaN; figure; pcolor(tc'); shading flat; axis off axis xy; colorbar; </code> This gives: {{ :analysis:course:week10_tc2.png?600 |}} Now we have a nice smooth estimate of this cell's place field. ☛ What does the ''occ_binned(occ_binned < 0.01) = 0;'' line accomplish? Uncomment it and re-run the code. ☛ What happens if you don't do ''tc(isinf(tc)) = NaN;''? ==== Bayesian decoding ==== The procedure of Bayesian decoding is illustrated in this figure (from [[http://www.cell.com/neuron/abstract/S0896-6273(10)00507-6 | van der Meer et al. 2010]]): {{ :analysis:course:mvdm_yorkfieldtrip13.png?900 |}} For this particular experiment, the goal of decoding is to recover the location of the rat, given neural activity in some time window. More formally, we wish to know $P(\mathbf{x}|\mathbf{n})$, the probability of the rat being at each possible location $x_i$ ($\mathbf{x}$ in vector notation, to indicate that there are many possible locations) given a vector of spike counts $\mathbf{n}$. If $P(\mathbf{x}|\mathbf{n})$ (the "posterior") is the same for every location bin $x_i$ (i.e. is uniform), that means all locations are equally likely and we don't have a good guess; in contrast, if most of the $x_i$ are zero and a small number have a high probability, that means we are confident predicting the most likely location. Of course, there is no guarantee that our decoded estimate will agree with the actual location; we will test this later on. So how can we obtain $P(\mathbf{x}|\mathbf{n})$? We can start with Bayes' rule: \[P(\mathbf{x}|\mathbf{n})P(\mathbf{n}) = P(\mathbf{n}|\mathbf{x})P(\mathbf{x})\] The key quantity to estimate is $P(\mathbf{n}|\mathbf{x})$, the probability of observing $n$ spikes in a given time window when the rat is at location $x$. At the basis of estimating this probability (the "likelihood" or evidence) lies the tuning curve: this tells us the //average// firing rate at each location. We need a way to convert a given number of spikes -- whatever we observe in the current time window for which we are trying to decode activity, 3 spikes for cell 1 in the figure above -- to a probability. In other words, what is the probability of observing 3 spikes in a 250ms time window, given that for this location the cell fires, say at 5Hz on average? A convenient answer is to assume that the spike counts follow a Poisson distribution. Assuming this enables us to assign a probability to each possible spike count for a mean given by the tuning curve. Specifically, from the [[http://en.wikipedia.org/wiki/Poisson_distribution | definition of the Poisson distribution]], it follows that \[P(n_i|\mathbf{x}) = \frac{(\tau f_i(\mathbf{x}))^{n_i}}{n_i!} e^{-\tau f_i (x)}\] $f_i(\mathbf{x})$ is the average firing rate of neuron $i$ over $x$ (i.e. the tuning curve for position), $n_i$ is the number of spikes emitted by neuron $i$ in the current time window, and $\tau$ is the size of the time window used. Thus, $\tau f_i(\mathbf{x})$ is the mean number of spikes we expect from neuron $i$ in a window of size $\tau$; the Poisson distribution describes how likely it is that we observe the actual number of spikes $n_i$ given this expectation. In reality, place cell spike counts are typically not Poisson-distributed ([[http://www.ncbi.nlm.nih.gov/pubmed/9501237 | Fenton et al. 1998]]) so this is clearly a simplifying assumption. There are many other, more sophisticated approaches for the estimation of $P(n_i|\mathbf{x})$ (see for instance [[http://www.ncbi.nlm.nih.gov/pubmed/17925266 | Paninski et al. 2007]]) but this basic method works well for many applications. The above equation gives the probability of observing $n$ spikes for a given average firing rate for a single neuron. How can we combine information across neurons? Again we take the simplest possible approach and assume that the spike count probabilities for different neurons are independent. This allows us to simply multiply the probabilities together to give: \[P(\mathbf{n}|\mathbf{x}) = \prod_{i = 1}^{N} \frac{(\tau f_i(\mathbf{x}))^{n_i}}{n_i!} e^{-\tau f_i (x)}\] An analogy here is simply to ask: if the probability of a coin coming up heads is $0.5$, what is the probability of two coints, flipped simultaneously, coming up heads? If the coins are independent then this is simply $0.5*0.5$. Combining the above with Bayes' rule, and rearranging a bit, gives \[P(\mathbf{x}|\mathbf{n}) = C(\tau,\mathbf{n}) P(\mathbf{x}) (\prod_{i = 1}^{N} f_i(\mathbf{x})^{n_i}) e^{-\tau f_i (\mathbf{x})} \] This is more easily evaluated in vectorized MATLAB code. $C(\tau,\mathbf{n})$ is a normalization factor which we simply set to guarantee $\sum_x P(\mathbf{x}|\mathbf{n}) = 1$ (Zhang et al. 1998). For now, we assume that $P(\mathbf{x})$ (the "prior") is uniform, that is, we have no prior information about the location of the rat and let our estimate be completely determined by the likelihood. === Preparing tuning curves for decoding === With the math taken care of, we can now start preparing the data for the decoding procedure. First we need to make sure we have tuning curves for all neurons. Now we need to do the same for //all// cells. For now, we will revert to using the "low-resolution" version (with 63x47 bins) with a small amount of smoothing. Even though this is not as good of an estimate as the high-resolution version, our decoding will be super slow if we try to run it on a high-resolution smoothed estimate. So first, let's inspect our updated tuning curve example: <code matlab> kernel = gausskernel([4 4],2); % 2-D gaussian, width 4 bins, SD 2 SET_xmin = 10; SET_ymin = 10; SET_xmax = 640; SET_ymax = 480; SET_nxBins = 63; SET_nyBins = 47; spk_binned = ndhist(cat(1,spk_x',spk_y'),[SET_nxBins; SET_nyBins],[SET_xmin; SET_ymin],[SET_xmax; SET_ymax]); spk_binned = conv2(spk_binned,kernel,'same'); % smoothing occ_binned = ndhist(cat(1,Data(sd.x)',Data(sd.y)'),[SET_nxBins; SET_nyBins],[SET_xmin; SET_ymin],[SET_xmax; SET_ymax]); occ_mask = (occ_binned < 5); occ_binned = conv2(occ_binned,kernel,'same'); % smoothing occ_binned(occ_mask) = 0; % don't include bins with less than 5 samples VT_Fs = 30; tc = spk_binned./(occ_binned .* (1 / VT_Fs)); tc(isinf(tc)) = NaN; pcolor(tc'); shading flat; axis xy; colorbar; axis off; </code> Then, we can do the same for all cells in our data set: <code matlab> clear tc nCells = length(sd.S); for iC = 1:nCells spk_x = interp1(Range(sd.x),Data(sd.x),Data(sd.S{iC}),'linear'); spk_y = interp1(Range(sd.y),Data(sd.y),Data(sd.S{iC}),'linear'); spk_binned = ndhist(cat(1,spk_x',spk_y'),[SET_nxBins; SET_nyBins],[SET_xmin; SET_ymin],[SET_xmax; SET_ymax]); spk_binned = conv2(spk_binned,kernel,'same'); tc = spk_binned./(occ_binned .* (1 / VT_Fs)); tc(isinf(tc)) = NaN; sd.tc{iC} = tc; end </code> Note that we don't need to recompute the occupancy because it is the same for all cells. Let's inspect the resulting tuning curves: <code matlab> ppf = 25; % plots per figure for iC = 1:length(sd.S) nFigure = ceil(iC/ppf); figure(nFigure); subplot(5,5,iC-(nFigure-1)*ppf); pcolor(sd.tc{iC}); shading flat; axis off; end </code> You will see a some textbook "place cells" with a clearly defined single place field. There are also cells with other firing patterns. ☛ One cell has a completely green place map. What does this indicate, and under what conditions can this happen? Since we see cells with fields in some different locations, it seems unlikely that a single sensory cue or nonspatial source can account for this activity. Of course, numerous experiments have demonstrated that many place cells do not depend on any specific sensory cue to maintain a stable firing field. === Preparing firing rates for decoding === The tuning curves take care of the $f_i(x)$ term in the decoding equations. Now we need to get $\mathbf{n}$, which are simply spike counts: <code matlab> clear Q; binsize = 0.25; tvec = t(1):binsize:t(2); for iC = length(sd.S):-1:1 spk_t = Data(sd.S{iC}); Q(iC,:) = histc(spk_t,tvec); end nActiveNeurons = sum(Q > 0); </code> This "Q-matrix" of size ''[nCells x nTimeBins]'' is the start of a number of analyses, such as the nice ensemble reactivation procedure introduced in Peyrache et al. 2009. Let's inspect it briefly: <code matlab> % look at it imagesc(tvec,1:nCells,Q) set(gca,'FontSize',16); xlabel('time(s)'); ylabel('cell #'); </code> You should see: {{ :analysis:course:week10_qmat.png?600 |}} If you zoom in to a smaller slice of time, you will notice that there are gaps in the data, i.e. segments without any activity whatsoever. This is a quirk of this particular data set: the epochs when the rat is in transit between the pedestals and the track have been removed to facilitate spike sorting. The final step before the actual decoding procedure is to reformat the tuning curves a bit to make the decoding easier to run. Instead of keeping them as a 2-D matrix, we just unwrap this into 1-D: <code matlab> clear tc nBins = numel(occ_binned); nCells = length(sd.S); for iC = nCells:-1:1 tc(:,:,iC) = sd.tc{iC}; end tc = reshape(tc,[size(tc,1)*size(tc,2) size(tc,3)]); occUniform = repmat(1/nBins,[nBins 1]); </code> === Running the decoding algorithm === Aaandd... action! <code matlab> len = length(tvec); p = nan(length(tvec),nBins); for iB = 1:nBins tempProd = nansum(log(repmat(tc(iB,:)',1,len).^Q)); tempSum = exp(-binsize*nansum(tc(iB,:),2)); p(:,iB) = exp(tempProd)*tempSum*occUniform(iB); end </code> Note, on my 2-year old laptop this takes about 2 minutes to run. The computation can be speeded up substantially by only looping over bins that have non-zero occupancy, but this would complicate the code a bit. ☛ Compare these steps with the equations above. There is no log in there; why does it appear here? Finally we renormalize and set the posterior time windows without any neural activity to zero (done for display purposes): <code matlab> p = p./repmat(sum(p,2),1,nBins); p(nActiveNeurons < 1,:) = 0; </code> === Visualizing the results === The hard work is done. Now we just need to display the results. Before we do so, we should convert the rat's actual position into our binned form, so that we can compare it to the decoded estimate: <code matlab> xBinEdges = linspace(SET_xmin,SET_xmax,SET_nxBins+1); yBinEdges = linspace(SET_ymin,SET_ymax,SET_nyBins+1); xTempD = Data(sd.x); xTempR = Range(sd.x); yTempD = Data(sd.y); gS = find(~isnan(xTempD) & ~isnan(yTempD)); xi = interp1(xTempR(gS),xTempD(gS),tvec,'linear'); yi = interp1(xTempR(gS),yTempD(gS),tvec,'linear'); xBinned = (xi-xBinEdges(1))./median(diff(xBinEdges)); yBinned = (yi-yBinEdges(1))./median(diff(yBinEdges)); </code> Now we can visualize the decoding (press Ctrl-C to break out of the loop): <code matlab> goodOccInd = find(occ_binned > 0); for iT = 1:size(p,1) cla; temp = reshape(p(iT,:),[SET_nxBins SET_nyBins]); toPlot = nan(SET_nxBins,SET_nyBins); toPlot(goodOccInd) = temp(goodOccInd); pcolor(toPlot); axis xy; hold on; caxis([0 0.5]); shading flat; axis off; hold on; plot(yBinned(iT),xBinned(iT),'ow','MarkerSize',15); h = title(sprintf('t %.2f, nCells %d',tvec(iT),nActiveNeurons(iT))); if nActiveNeurons(iT) == 0 set(h,'Color',[1 0 0]); else set(h,'Color',[0 0 0]); end drawnow; pause(0.1); end </code> Initially, when the rat is moving around at the base of the maze (its actual position is indicated by the white ''o''), no decoding results will be shown because there are no cells active. However, after a few seconds the rat starts running up the stem of the maze, and some pixels in the plot will change color, indicating $P(\mathbf{x}|\mathbf{n})$, the "posterior probability" that is the output of the decoding procedure. As you can see, the decoding seems to track the rat's actual location as it moves. ☛ No decoding is available for those bins where no neurons are active, because we manually set the posterior to zero. However, there also seem to be some frames in the animation where some neurons are active (as indicated in the title), yet no decoded estimate is visible. What is the explanation for this? === Exporting the results to a movie file === By making the MATLAB animation into a movie file, it is often easier to explore the results. To do this, we can run the animation code above, with a few small modifications. First, before entering the main plotting loop, set the figure to be used to a specific size: <code matlab> h = figure; set(h,'Position',[100 100 320 240]); </code> This is important first, to keep the size of the resulting movie file manageable (the above sets a 320x240 pixel figure size), and second, because many movie encoders (such as the excellent [[http://www.xvid.org/ | XVid]]) will only work with certain sizes. Next, we need to store each frame into a variable that we can later write to file. Modify the last two lines inside the loop to: <code matlab> f(iT) = getframe(gcf); % store current frame drawnow; </code> If you now run the code again, each frame gets stored in the ''f'' variable as the loop runs. Break out of the loop after a few seconds to test the writing-to-file part: <code matlab> fname = 'test.avi'; movie2avi(f,fname,'COMPRESSION','XVid','FPS',10); </code> The above will only work if you have the XVid codec installed: I highly recommend this because it creates movie files that are an order of magnitude smaller than uncompressed files. If you have trouble with XVid, you can of course still save an uncompressed file for now. For longer movies, it is often required to save a file, say, every 500 frames, to prevent the ''f'' variable getting too large. These segments can then be merged with a video editing program such as [[http://www.virtualdub.org/ | VirtualDub]] (Windows only AFAIK; please suggest OSX/Linux alternatives if you know any that work well!). ==== Assignment ==== Visual inspection of the animation or movie suggests that the decoding does a decent job of tracking the rat's true location. However, especially because of the number of parameters involved in the analysis (bin size, how firing rates are computed, the Poisson and independence assumptions, etc.) it is important to quantify how well we are doing. ☛ Modify the visualization code above to also compute a //decoding error// for each frame. This should be the distance between the rat's actual location and the location with the highest posterior probability (the "maximum a posteriori" or MAP estimate). Plot this error over time, excluding those bins where no cells were active. How does this error change over the course of the session? How does it change if you reduce the bin size for decoding to 100ms?