Linear track demo¶

Here we load sorted unit data obtained from MatClust to identify putative place cells. The session has about 15 minutes of linear track running, followed by 15 minutes of rest.

Data was recorded by Joshua Chu, with a Spikegadgets wireless headstage, on July 8th, 2017, from the CA1 area of a male Long-Evans rat named install.

Notebook was created by Etienne Ackermann.

Preliminaries¶

We use the following Python packages; they can be installed either with conda install <pkg> or with pip install <pkg>. Most of these packages are very common, and it is more than likely that you have them installed already.

It doesn't hurt to type pip install <pkg> if the package is already installed, so don't be nervous to try it out!

pip install nelpy

For downloading example data from the web:

1. requests   # used to download data from web
2. tqdm       # used to show progress of download

the rest of this notebook:

3. numpy      # numerical powerhorse for Python
4. matplotlib # used to make plots and figures
5. scipy      # signal processing, stats, etc.; used here to smooth signals
6. sklearn    # machine learning in Python; used here to do train-test split 
              # (under the hood) when evaluating Bayesian decoding performance

and of course:

7. nelpy      # Ephys object models, and analysis routines

Now that we have all the packages we nee, let's get started! First, we need to get the sample data...

1. Obtain example data¶

We will look for data in the example-data\linear-track\ directory inside your current working directory. If the data doesn't exist, we will download it from https://github.com/nelpy/example-data, and save it to your local machine.

If you already have the data, it won't be downloaded again.

In particular, we will download two files, namely

1. trajectory.videoPositionTracking which is a binary file with (x,y) position coordinate pairs and timestamps, and
2. spikes.mat which is a Matlab file containing information about sorted units (cells) obtained by using MatClust (https://bitbucket.org/mkarlsso/matclust).

2. Import packages¶

3. Extract position data¶

The position data is stored in a binary file, and we already made a reference to it as filename_pos, which should point to ./example-data/linear-track/trajectory.videoPositionTracking, inside your current working directory.

More information about the format of the file can be obtained at https://github.com/nelpy/example-data/tree/master/linear-track .

Here we define a simple function to read the plaintext header from the file, since the header also informs us of the format of the rest of the data.

Armed with this new information, we read over the header, and then extract (1) the 32 bit unsigned timestamp, as well as (2) the 16 bit unsigned position data into lists x1 and y1 (x2 and y2 were not used in this recording, and can be ignored).

Remark: up to this point, we have not used nelpy yet. We have our position as two lists x1 and y1, along with a list of timestamps. But working with these lists of coordinates and timestamps directly can be tedious and error-prone. Consequently, we will wrap the data into nelpy objects to make our life easier.

3.a) Get session boundaries¶

First, we will estimate epochs during which the animal was supposed to be running on the track, and those for which the animal was in its sleep box. Note that this task can be a little tricky to estimate directly from the data, since the animal is not always running while on the track, and may in fact just be stationary for significant periods.

We can (and probably should) of course get the experimantal epochs (that of running on the track, and that of resting in the sleep box) from our lab notes, but it's always good to make sure that your data agrees with your notes, and moreover, it demonstrates how we can use nelpy to estimate these epochs in case our notes were lost, or inaccurate.

Here we've used nelpy to get event boundaries, in sample numbers, where the animal spent at least 10 seconds with zero change in its x1 coordinate. We could do this, because the tracker does not track the animal when it is in the sleep box (it is outside of the camera's view).

Next, we create our first nelpy object, namely an EpochArray containing all the epochs defined by the bounds returned above.

In general it is not necessary to sort epochs like we do below, but we have to do that here because the way we build our session_epochs is a little convoluted: we build it by first estimating sleep box epochs, and then combine the sleep box epochs with the complement of the sleep box epochs. When we take a union of epochs, they are no longer guaranteed to be in sorted order, and so we sort them for good measure.

We can always inspect nelpy objects by printing them to screen like so:

Update: we used to have three epochs as a result of the above construction, but somewhere along the line we changed an implementation detail in nelpy so that when contiguous epochs are added together, they are automatically merged. This is sometimes very useful, and sometimes quite frustrating. It's unclear which approach is superior, but the automatically merged case is better for a majority of underlying epoch array processing functions. If this is an issue, you can always subtract epsilon from the right endpoint of your epochs to prevent them from being viewed as contiguous. This is a hacky-but-effective fix, with very low chances of having any unintended consequences.

It used to be that we had three resultant epochs, two in the rest state, and one while running on the track. But since all the epochs are contiguous, we now only have a single epoch. Consequently, to demonstrate what I wanted to demonstrate (before the change in nelpy), I'll artificially reduce the lengths of the epochs such that they are no longer contiguous:

which tells us that we have 3 epochs, with a total duration of approximately 30 minutes.

HOWEVER! We expected to only have two epochs! One 15 minute run, followed by a 15 minute sleep box session. So what happened? Well, we can have a closer look at the epochs, by printing out their durations:

We see that the first epoch is only about 25 seconds long, and this corresponds to when we moved the animal onto the track, so we're happy to ignore that particular epoch, and focus on the other two.

Note also that nelpy returned a nice human readable duration for each epoch. Under the hood, the durations are stored (calculated, more accurately) in seconds, and we could have computed it ourselves like so:

We see that the durations that we computed above agree with the ones we printed out earlier, but 15:59:115 minutes is arguably easier to parse than 959.11506667 seconds.

In fact, we don't even really need to call print() on many of the objects, and we can just invoke them to see a representation of the objects. Let's try that on our sleep_box object:

We see that (as we already know) there are two epochs during which we think the animal was in its sleep box, namely the first 25 seconds, and then the last 16 minutes of the recording session.

We can similarly look at the complement of the sleep box epochs, which should then be (naturally) when the animal was on the linear track:

3.b) Create 2D trajectory object¶

Next we will use our lists x1 and y1 to build a nelpy AnalogSignalArray. This is pretty easy: we pass in the position coordinates as a 2 x n_samples array, along with the corresponding timestamps, and we specify that we want this object to be defined on the epoch where the animal was NOT in its sleep box, i.e., when the animal was on the track.

Note also that we defined the trajectory object on a slightly smaller epoch, namely we shrunk the epoch by 20 seconds from both directions. We did this in order to be sure that the first and last parts of the trajectory were not those where we were moving the animal onto and off of the track.

We can plot the position/trajectory object like so:

Recall that the linear track epoch ranged from approcimately 4420 to 5380 seconds (We can always get these by simply inspecting pos again, like blow):

and so we can easily plot the trajectory, say, for the first 30 seconds while the animal was on the track. One way to do this, is to define an epoch, and then to restrict the position object to that newly defined epoch. Let's do this now:

This way of filtering / restricting the data is very powerful, and we will use it later on to get spikes only during times when the animal was running faster than some threshold, and so on.

Before we move onto the next phase (linearization), let's take a quick aside to make the plot above a bit nicer. First, we may want to set the aspect ratio to equal, and we may also want to add titles, labels, etc. We may even want to smooth the trajectory.

3.c) Linearize position¶

We can do the rest of our analysis in 2D, but sometimes it's nicer to linearize the position data first. Here, we simply use PCA (a dimensionality reduction technique) to find the direction of maximal variance in our position data, and we use this as our new 1D linear track axis.

NOTE: To use PCA from scikit-learn, we need to provide our data as an n_samples x n_features matrix. The internal representation of the nelpy AnalogSignalArray can actually be changed, but we can always explicitly ask to get the underlying data back in a particular format. We do this, by asking for _ydata_colsig, which means that we want each feature (or signal, or component) to be a column of the data matrix. Here you shouldn't be concerned with the favt that we asked for _ydata: this is simply an internal naming convention of the AnalogSignalArray, and has nothing to do with our x and y coordinates. More precisely, _ydata refers to the entire data matrix, in contrast to _tdata which stores the timestamp info. So _ydata contains both x1 and y1.

We also re-scale our data to range from 0 to 100 (so that we can express movement along the track as a percentage). Typically we would actually scale our data to some physical dimension or unit, such as cm, and not pixels as we have it above, or percentage, as we're transforming it into here.

We now have Xlinear as a scaled (0 to 100) position variable, as a numpy array. We want to put it back into a nelpy AnalogSignalArray container, so that we can slice and interogate it as before:

We see that pos1d now only has 1 signal, as expected, and if we plot it, we can confirm that it ranges between 0 and 100:

4. Extract spike times¶

We recorded extracellular activity in region CA1 of the hippocampus from a male Long-Evans rat, and we already performed spike detection and sorting (clusterting) in a separate program named MatClust (https://bitbucket.org/mkarlsso/matclust). MatClust is Matlab-based, and it gives us a .mat file with containing our sorted units and their associated spike times.

Luckily, we can read .mat files directly into Python, but the file format is still a little clunky, and so we need to parse it into a more intuitive format fisrt. For more information on the format of the .mat file, you can refer to https://github.com/nelpy/example-data/tree/master/linear-track.

At any rate, we extract spike times below, so that we end up with a list of lists, where the ith list is a list of the spike times associated with the ith unit. For example,

spikes = [[1, 2, 5, 10], [3, 5], [2, 9, 20]]

contains three units (three inner lists), with the spike times for the 1st unit being [1, 2, 5, 10] and those for the last unit being [2, 9, 20]. Spike times are given in seconds, so floating point numbers are used.

As before, we now put the list of spike times into a nelpy container object to make it easier to interact with the data. In particular, we put the spikes into a SpikeTrainArray. But first, let's also get out the time boundaries for which the spikes were recorded. We can approximate these boundaries by looking for the first and last recorded spikes, but the .mat file actually contains this information explicitly, so we'll get it directly from there:

Now we're ready to build our SpikeTrainArray. We simply pass in the list of spike times, the epoch during which spikes were recorded, and the sampling rate (recall that we previously specified the sampling rate, FS=30000 Hz).

Printing out our SpikeTrainArray (st) tells us that we have 31 units, sampled at 30,000 Hz.

We can ask several things from this SpikeTrainArray object, such as the number of spikes for each unit. To do that, we simply access the .n_spikes property, which returns the number of spikes associated with each unit as an array:

Of course, we can also compute these things manually. Let's take a step back, and compute the average firing rate for each unit, using not the nelpy SpikeTrainArray, but the original list of lists:

This wasn't too difficult, but nelpy can make this a little easier.

Like all other core nelpy objects, a SpikeTrainArray has a support on which it is defined, and so we can always get timing information by accessing the underlying support object (instead of using start and stop explicitly, as we did above).

Let's calculate the average firing rates using nelpy now:

Now wasn't that a little simpler? But the real beauty comes in when the SpikeTrainArray is not simply defined on one continuous epoch, but on many smaller epochs. For example, we might be interested in asking "what is the average firing rate for each unit, during periods when the animal was running?"---in that case, the firing rate calculation using the list of lists approach would become exceedingly painful, but we would be able to simply write

avg_firing_rates_during_run = st[run_epochs].n_spikes / st[run_epochs].support.duration 

which will work just the same, even though st[run].support now consists of many discontiguous epochs.

Visualize spikes with raster plots¶

Nelpy has several built-in plot types (we have already seen npl.plot2d() before), and one that is frequently useful is the rasterplot().

We notice that there is a clear difference in overall activity between the run session (4400--5400 seconds) and when the animal was in its sleep box (5550--). For example, units 4, 7, 8, 26, 27, and 29 seem to be largely inactive during run, but active during the sleep box. Let's use the rasterplot() to highlight those units:

The rasterplot() function has quite a bit of flexibility. As another example, let's plot the same as above, but let's collapse (stack) all the units onto each other:

We can even get a histogram with the number of spikes across all units by using the npl.rastercountplot() function. Note that this function is not yet fully developed, and still needs to be modified to increase flexibility.

One final raster example, where we call rastercountplot() on the list of units that seem to be more active in the sleep box than on the track:

5. Estimate tuning curves¶

Next, we will use the position and spikes to estimate the spatial tuning curves of our 31 units.

Place fields are usually estimated using only those epochs when the animal is running faster than some speed threshold, so we find these epochs first. We smooth our speed estimates with a Guassian kernel with sigma = 100 ms, and then we return those epochs where the animal was running at least 8 percent per second.

Note: Usually we would measure the speed in cm/s or similar, but recall that our position has been scaled to be a percentage; the track was about 117 cm long, so that percent/s is not too different from cm/s.

The get_run_epochs() function has two thresholds (both set to 8 above). More specifically, the first threshold (v1) is referred to as the primary threshold, and this must be reached or exceeded for an epoch to be considered a candidate. But then the secondary threshold (v2) is used to determine the epoch boundaries.

For example, if we had used v1=8, v2=5 then only epochs during which the animal reached a speed of at least 8 u/s would be considered, and the epoch boundaries would be when the animal ran ran at or slower than 5 u/s, starting at the peak speed of >8 u/s.

There are other possible arguments, too, so that we can e.g. ask that the animal must run at least 8 u/s for some minimum period of time, and so on.

Let's see for how much time the animal actually ran faster than 8 percent per second:

We see that there are 221! short epochs during the animal ran sufficiently fast, and that the total amount of time spent running at this speed or faster, is about 4 minutes and 18 seconds (out of approximately 16 minutes spent on the track).

As promised before, we can now easily (trivially) calculate the firing rates of all the units, only during times when the animal was running:

We can also plot the linearized trajectory only during those run epochs. Each new color indicates that it is a new epoch. The black on the bottom figure corresponds to epochs when the animal did not meet the 8 percent per second requirement, and so those are when the animal is considered to be at rest.

Now that we have the epochs during which the animal was running, we can use it to restrict our SpikeTrainArray to those epochs. Then we will use these spikes, along with the position data, to estimate the spatial tuning curves.

We bin our spikes into 50 ms bins, so that we can count the number of spikes in a small window of time, which we will then later associate with the particular position bin that the animal was at when those spikes occured.

We also apply a little bit of spike time smoothing.

We now estimate the tuning curves using the binned spikes during run (bst_run), and we partition our linear track into n_extern=50 equal-sized spatial bins, and we specify that the linearized position ranges from extmin=0 to extmax=100. We also smooth the estimated tuning curves with a 0.2 cm smoothing kernel (this should be imperceptable, and a larger value for sigma should really be used, but we can always apply the smoothing later, so it's not a real problem).

Finally, we reorder the units by their peak firing positions on the track, simply to make visualization a little nicer:

That's it! Let's see what we've got:

We see that some units have spatially localized tuning curves, and some don't. Also, some units have large firing rates, and some don't. To see the shape of all the tuning curves better, we can normalize the peak firing rates like so:

This normalized view allows us to see the shapes of the tuning curves more clearly, and we can see that units 4, 24, and 7 are essentially flat, meaning they have NO spatial localization for their firing rates. This is partly due to there just being almost no spikes during run for those units though.

At any rate, let's filter out the units a little further, to try and figure out which ones might actually be place cells.

Restrict subset of cells (units) to use for subsequent decoding and/or analysis¶

Here we may request to use

1. unimodal cells only,
2. pyramidal cells only,
3. active cells only,
4. any combination of the above, and other criteria

Here we only impose a minimum peak firing rate of 1.5 Hz, and we reject putative interneurons by imposing a maximum average firing rate of 5 Hz averaged over the entire track.

6. Evaluate decoding performance¶

Next, having defined a subset of 12 place cells, we may reasonably want to ask how well these cells can represent the animal's location? We can evaluate the expected decoding performance by using a Bayesian decoder, and by evaluating the decoding accuracy on a test set.

Here, we use 5-fold cross-validation and the tuning curves of our 12 place cells to evaluate the performance, and the results are summarized in the figure below:

We can evaluate cumhist at any value between p=[0,1], which will return the error at that percentile. For example, to get the median decoding error, we simply evaluate cumhist(0.5):

This tells us that, with probability 0.9, we can decode the animal's position to within 20 percent of the true position. Recall that in our particular example, percent and cm are almost the same. The inset shows a zoomed-in view, from where we can see that with probability 0.7, we have an error of less than 8 percent.

It's not the best decoding accuracy imaginable, but overall, not bad!

Some additional views¶

Next, we take a look at where on the track spikes from a particular unit occurred:

Conclusion¶

This demo is still a work in progress, but hopefully it highlights a few ways in which using nelpy might make it easier to combine and analyze ephys data with different sampling rates, and different epochs where they are defined. There's a lot more that we can do with nelpy, so feel free to dig a little deeper to explore!

Feedback about this notebook would be appreciated, and can be sent to era3@rice.edu.

Supplementary analysis: hidden Markov models (HMMs)¶

Preliminaries¶

Nelpy makes use of hmmlearn (https://github.com/hmmlearn/hmmlearn), but modified to include a Poisson emissions model. The modified version can be found at https://github.com/eackermann/hmmlearn . In addition to the Poisson emissions model, this modified branch also fixed up an issue in calculating the log probabilities, as well as consistency with random seeds. The random seed issue has since been fixed in the main repository, and it is possible that the log probability issue has also been fixed, and if not, I (Etienne) should make a pull request to incorporate the fixes into the main branch. At any rate, for now, nelpy makes use of the fork at https://github.com/eackermann/hmmlearn .

Installing hmmlearn should be trivial on Linux, fairly easy on MacOS, and doable on Windows. To get it to compile on Windows, all you really need is to download and install the free version of Microsoft Visual Studio (with the Python extensions), and then hmmlearn should compile with no issues. After installation, you can remove Visual Studio again. However, if you don't want to go through the trouble of compiling it yourself, there is also a 64-bit Windows binary in the nelpy repository.

Overview¶

Here we give a quick demo of how we may use nelpy to learn hidden Markov models. Typically, HMMs are used in the absence of any [observable] behavioral correlate, but here we demonstrate the approach on the run data, where we do actually have access to the position data. Nevertheless, we don't use the position data when training the HMM, nor do we need it to decode neural activity to the state space. But it does make it nice for a demo, since we can use the position data to interpret what those hidden states might represent.

Briefly, we will use a SpikeTrainArray during bouts of running activity to learn our HMM. The HMM makes use of vectors of spike counts, so we need to bin our SpikeTrainArray into an appropriate bin size. Since we are learning a model on run data, we want to choose a bin size that's small enough to capture the animal's behavioral dynamics, but large enough to capture sufficiently many spikes (lot of empty bins will do us no good). Here, we will use 250 ms bins.

We will also apply a little bit of spike smoothing (300 ms standard devitation Gaussian kernel), since this has been shown to improve decoding accuracy (in general, not just for HMMs; the idea is roughly that there is some inherent noise in the signals, and so modeling that uncertainty with spike smoothing actually improves our model robustness).

Thereafter, we train our HMM on a subset of the BinnedSpikeTrainArray (the train set), and we can evaluate our model on the test set. We demonstrate how to decode to the state space, and also how to decode back to an actual behavioral correlate (in this cas position, if we augment our HMM with an additional mapping that we've learned between the state space and the external behavioral correlates.

There are much more that we can do with HMMs, and this section is by no means meant to be comprehensive, nor definitive, but hopefully it demonstrates how to get started, and how to do basic inference.

Estimate virtual tuning curves¶

Decoding to the state space:¶