Contents

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 demonstrated 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? Decoding approaches can help address such questions, and many more besides.

To develop an intuition for how decoding works, consider Hubel and Wiesel's demonstration 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:

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, described in van der Meer et al. (2017).

Setup

To view this notebook, I used the following steps:

• install Anaconda with default Python 3 environment
• you should now be able to open this notebook (jupyter notebook launches jupyter and opens a browser)

The notebook uses a number of default Python libraries (e.g., numpy, matplotlib) which come in default anaconda installations. It also uses a library we are developing for use with neural data, nelpy. To install nelpy, you should be able to pip install nelpy from a terminal window.

Setup continued

Finally, you'll need the data we use for this notebook. The source for this data is the CRCNS.org neural data repository, where it was deposited by the laboratory of Loren Frank. A large compressed archive can be found here. You should download it into a data directory. We then need to set the data location.

Loading and inspecting the data

Next, we load the spike and position data. These have been stored the 'FrankLabPickle.npz' file. You can load it using numpy.load. The data are in a dictionary called data, organized by animal nickname, experiment day, and data type (pos or spikes).

The Frank Lab archive does not distinguish between the quality of sorted cells. In general, you would not want to use neurons which were not well-sorted to make claims about single neuron properties; however, for this analysis 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!

In these data, there are multiple runs each day in different mazes. (You can find out the the description of these runs by looking in the taskinfo dataframe, i.e., jar.taskinfo.query('Animal=="Bon" & Day==3').) We're interested in spatial coding, so we want to look only at the data taken while running in the same maze. For the example of Bon/day 3, epochs 2 and 4 are runs in "TrackB".

For time series variables, the nelpy library incorporates a concept of temporal support (i.e., periods during which data are defined). This is defined using the EpochArray type. For the position data, the different sleep, run epochs define a 7 member EpochArray, which can be accessed using pos.support.

Using temporal support in nelpy

The position data is of type AnalogSignalArray, which nelpy uses for time series data. In this case, the two "channels" correspond to the x and y position dimensions. Note that an AnalogSignalArray is not just a simple numpy array (though you can access the data that way). Rather, the two dimensions of slicing for AnalogSignalArrays are epoch and channel. This means that asa[[2,4]] corresponds to data during the 3rd and 5th epochs in the support of asa, asa[:,3:6] corresponds to data from all epochs from the 4th through 7th channels, and asa[[2,4],3:6] implies both conditons. (Remember that Python is 0 indexed rather than 1 indexed like Matlab!)

Let's get the position for the two runs and plot them.

In these sets of experiments, spike sorting was done combining all sessions in a day. We store the spike times for each neuron in a nelpy SpikeTrainArray object. SpikeTrainArray is a similar structure to AnalogSignalArray in that it has an EpochArray as support, and can be sliced. The slicing for a SpikeTrainArray object is slightly different, however. The first dimension is still epochs. The second dimension has numeric labels (UnitIds) - this allows information about neurons to be preserved while still allowing for things like reordering for presentation. To index using the UnitId labels, use sta.loc[], but if you just want to index by position use sta.iloc[]. For example, if you want to select the first 10 units, you should use sta.iloc[:,:10], but if you want UnitIDs 1 and 5, you can use sta.loc[:,[1,5]] regardless of whether the order of the units has been changed or some have been removed.

Now we can plot the position data (location of the rat as it runs the maze, in gray), and the position of the rat when an example neuron fired a spike (in red). We'll use a fun AnalogSignalArray feature, the at option for support.

How does the at option work?

The at option is a convience wrapper which creates a scipy.interpolate.interp1d() and then calls it. In this case this amounts to something like

from scipy import interpolate
tnew = S[[1,3],unit_id].time
t,x = pos.asarray()
f = interpolate.interp1d(x=t, y=x, axis=-1, kind='linear', copy=True, bounds_error=False)
xnew = f(tnew)

This cell seems to have a place field on the corner leading to the left arm of the maze. There are also a few spikes at the ends of the maze, where the rat receives his reward as he alternates on the track.

Estimating tuning curves

This figure is a useful visualization of the raw data, but it is not a tuning curve, which captures the relationship between a variable of interest (e.g. position) to firing rate. A set of tuning curves can be thought of as an encoding model that specifies how the position variable is encoded by our population of place cells. As a first step to estmiating this model, we restrict our data to those times the rat is actually running on the track. Note the use of the complement of an EpochArray in the second figure. We were required to create a new variable defined just over the two run epochs ([1,3]), so that the domain of the support would be properly defined for the complement. (Try what happens if you don't do this to see how it pulls in position data from other epochs!)

The nelpy plotting routines include a very nice raster plot function as well as functions which make it easy to color in different epochs. Lets show all of our spikes, with a light blue window over run periods.

Binned spike trains

We will calculate tuning curves using binned spike trains. This is not necessary (and, in the 1D case, tuning curve estimation is easy without it). However, if we want to smooth the spike trains, binning makes this much more easy. nelpy includes a BinnedSpikeTrain class, which is a special time series for spike trains. Below, we bin spikes into narrow time windows, smooth with a Guassian kernel, and then rebin into larger timebins. The nelpy functions are careful to just bin/smooth in defined epochs, respecting epoch boundaries.

Estimating tuning curves

Now, we can compute tuning curves for each of our cells. The TuningCurves() function does this in three steps:

• Compute an occupancy histogram (amount of time spent in each 2-D bin)
• Compute a spike histogram (count number of spikes fired in each 2-D bin) - we'll do this for smoothed/binned spike trains
• Compute firing rates by dividing the spike count by the occupancy

The below code plots the output of the third step for our example cell:

Smoothing tuning curves

As can be seen from the plot above, the (arbitrary) binning results in a noisy-looking tuning curve. As shown in van der Meer et al. (2017) applying some judicious smoothing can increase decoding accuracy substantially. This smoothing is kind of a hack, and can be viewed as an approximation of a more principled estimation process that starts with a prior on the likely shape of single cell tuning curves which is then updated using the encoding data. For now, we'll just do what works:

LoadMetadata; % loads experimenter-annotated file associated with each data session

% ENCoding variables: used for estimating tuning curves (encoding model) ENC_S = restrict(S,metadata.taskvars.trial_iv); % trial_iv contains the start and end times of trials ENC_pos = restrict(pos,metadata.taskvars.trial_iv);

% check for empties and remove keep = ~cellfun(@isempty,ENC_S.t); ENC_S = SelectTS([],ENC_S,keep);

% also set up DECoding variables for use later DEC_S = SelectTS([],S,keep);

You can inspect the tuning curves of all cells as follows:

You will see a some textbook “place cells” with a clearly defined single place field. There are also cells with other firing patterns.

This notebook computes tuning curves for location, but the idea is of course more general. For continuous variables in particular, it is a natural and powerful way to describe the relationship between two quantities – spikes and location in this case, but there is no reason why you couldn't do something like pupil diameter as a function of arm reaching direction, for instance!

Note: The tc object returned by TuningCurve2D() contains a number of other fields in addition to the actual tuning curves plotted above. One very useful field is tc.occupancy. Later, we will need to parse bins by occupancy in order to speed up decoding computations. In a nutshell, we will want to only consider spatial bins with sufficient sampling (occupancy). In addition, you can find the boundaries of the x and y axis bins (tc.xbins, tc.ybins) as well as bin centers, which we've already used for plotting. Spatial information and spatial sparsity can also be calculated to look for place-iness.

Bayesian decoding

As noted in the introduction above, given that we have neurons whose activity seems to encode some stimulus variable (location in this case, as evident from the tuning curves above), we can attempt to decode that variable based on the neurons' time-varying activity.

A popular approach to doing this is “one-step Bayesian decoding”, illustrated in this figure (from van der Meer et al. 2010):

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})$

If you have not come across Bayes' rule before, or the above equation looks mysterious to you, review the gentle intro by linked to at the top of the page. In general, it provides a quantitative way to update prior beliefs in the face of new evidence.

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 firing rate given by the tuning curve. For instance, here are the probabilities of observing different numbers of spikes $k$ (on the horizontal axis) for four different means ($\lambda = $1, 4 and 10):

In general, from the 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 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 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 \sum_{i = 1}^N 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.

The tuning curves take care of the $f_i(x)$ term in the decoding equations. Next, we need to get $\mathbf{n}$, the spike counts.

Obtaining the Q-matrix

This “Q-matrix” is of size [nCells x nTimeBins] and contains the spike count for each neuron in a given time bin. We've alreayd calculated this above as bst_run!

Each column of this matrix is a vector of spike counts ($\mathbf{n}$ in the section above). You can use the interactive features in the plotly version to zoom in to specific columns in the above plot.

Note that this type of matrix is shared by many different kinds of data in neuroscience -- this one contains binned spike counts, but a similar matrix can be constructed for spike density functions (firing rates), BOLD signals for voxels measured by fMRI, calcium levels from fluorescent microscopy, and so on. Analyses such as representational similarity analysis, functional and effective connectivity, discovery of network structure (hubs, ensembles, cliques, etc.) take a matrix of this type as a starting point.

Running the decoding algorithm

For the purposes of this tutorial, this matrix completes the set of ingredients required to perform decoding. Although the codebase provides configurable wrapper functions for this, a minimal version of the raw code follows below so you can see how it works:

Now we want 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. We can again use the nifty at flag in asarray().

Visualizing decoder output

Now we can visualize the decoding, shown here as a movie (see below for the code):

You can see the rat's true location, rendered as a yellow circle, moving around in the continuous-alternation pattern characteristic of the "W"-maze. The output of the decoder is plotted using the "seismic" colormap, with red colors corresponding to maximum probability. You should see that as the rat runs along the track, the decoder output tends to follow the rat's true location. In other words, based on neural activity in a given time bin, we can infer the location of the rat.

Note that there are times, especially when the rat is at the ends of the arms, when the decoding is particularly off. Why might this be? (Note that we have not restricted data for decoding to only running periods!)

Also note that in the Karlsson, Frank 2008 and 2009 papers that use this data, Mattias Karlsson used "linearized" position rather than 2D position. By clamping decoding to a 1D track, quality can be improved!

The code that generates the video is below.

Quantifying decoding error

The above movie gives a useful visual impression of what the decoder output looks like, and that it roughly seems to track the rat's position. A more systematic approach to evaluating the performance of the decoder is to plot the error, i.e. the distance between the actual and decoded positions. The cells below will estimate the decode error and a histogram of error as a function of position.

In this data, the results are not quite as cool as in MVDM's. It looks like the decoding error is on average larger on right side of the maze.

☛ What could be some reasons for this? Can you think of ways to test your suggestions?